# Spatial Kernel in Julia

**URL:** <https://discourse.julialang.org/t/spatial-kernel-in-julia/107707>\
**Category:** Geo\
**Tags:** question, package\
**Created:** [December 16, 2023, 11:08am UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707 "2023-12-16T11:08:38Z")\
**Posts on this page:** 9\
**Page:** 2

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [December 20, 2023, 9:05pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/22 "2023-12-20T21:05:58Z")

</div>

`length(Bkde.x)` must be equal to `size(Bkde.density, 2)` or `size(Bkde.density, 2) + 1`. Same for y. Difference is grid registration, a must know feature that is not GMT specific but valid for all grids. See [4. Standardized command line options — GMT 6.5.0 documentation](https://docs.generic-mapping-tools.org/dev/reference/options.html#grid-registration-the-r-option)

---

<div class="post-metadata">

**Author:** ![anjelinejeline](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/anjelinejeline/32/205552_2.png) [@anjelinejeline](https://discourse.julialang.org/u/anjelinejeline)\
**Post date:** [December 20, 2023, 9:12pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/23 "2023-12-20T21:12:35Z")

</div>

Like this?

```julia
G = mat2grid(Float32.(Bkde.density), x=(Bkde.density, 2), y=(Bkde.density, 2), proj4=“+proj=moll +lon_0=0 +x_0=0 +y_0=0”)

```

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [December 20, 2023, 9:37pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/24 "2023-12-20T21:37:25Z")

</div>

Changing the resolution to 1km with Rasters.jl is a single command…

```julia
rast_1k = resample(rast; res=1000, crs=your_crs_in_meters)

```

I’m not sure why you would use GMT.jl for that when you already had a `Raster`.

(But yes just going directly to 1000 is better than resampling - But you would just wrap the output of KernelDensity as a `Raster`?)

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [December 20, 2023, 10:17pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/25 "2023-12-20T22:17:34Z")

</div>

> [@Raf](#):
>
> I’m not sure why you would use [GMT.jl](https://juliahub.com/ui/Packages/GMT) for that when you already had a `Raster`.

My fault, but after a couple of days and the problem not solved I decided to offer an alternative.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [December 20, 2023, 10:27pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/26 "2023-12-20T22:27:27Z")

</div>

Yeah, I’ve never used KernelDensity.jl so can’t really help with that part.

But writing a raster is the easy bit, probably just one line after you work out getting a matrix from the Interpolations.jl `pdf`.

Like:

```julia
pdf_matix = ... # whatever you need to do with KernelDensity.jl / Interpolations.jl
write("some.tif", Raster(pdf_matrix, (X(xs), Y(ys)); crs))

```

We could probably add a KernelDensity.jl extension to Rasters.jl so these things are easier, and it just takes any Table/DataFrame and outputs a kernel density `Raster`.

---

<div class="post-metadata">

**Author:** ![anjelinejeline](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/anjelinejeline/32/205552_2.png) [@anjelinejeline](https://discourse.julialang.org/u/anjelinejeline)\
**Post date:** [December 21, 2023, 7:03am UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/27 "2023-12-21T07:03:13Z")

</div>

Hi @Raf  
Thank you  
I eventually ended up with rasterizing the output of the pdf function (see code below)  
I still have two questions  
Looking at the line you just wrote you specifiy Raster(pdf\_matrix, (X(xs), **Y(ys)**) instead of **Y(reverse(ys))** …  
I also did not reverse the y coordinate and I do not know why this was suggested by @Fliks can anyone explain me?

Also, after masking the raster and writing it into the disk I noticed that the missing values are transformed into -Inf

This is the warning message I got  
**Warning: `missing` cant be written with gdal, missinval for `Float64` of `-Inf` used instead**

I woud like to substitute them with NA but I can pass only strings to replace\_missing

> **[API Reference - Rasters](https://rafaqz.github.io/Rasters.jl/dev/reference/#Rasters.replace_missing-Tuple%7BAny%7D)**
>
> Rasters

So ideally I should replace them with “NA”

The thing is that I would like to import them in R afterwards where missing values are coded with NA …

what do you suggest?

Thank you very much for all suppport

```julia
# Extract x and y columns
x_values = df.x
y_values = df.y

# Extract x and y coordinates
coordinates = hcat(x_values, y_values)

bw=30000::Int64

# Defining the bounding of the World ESRI:54009 Mollweide 
# and resolution of 1kmx1km (1000m x 1000m)

res=1000
xmin=-17601618
xmax=17601617
ymin=-9018991
ymax=8751339

lengthx=round(Int,(xmax + abs(xmin)) / res)
lengthy=round(Int,(ymax + abs(ymin)) / res)

rx = KernelDensity.kde_range((xmin,xmax), lengthx)
ry = KernelDensity.kde_range((ymin, ymax), lengthy)

# Create a bivariate kernel density estimator
Bkde = kde(coordinates,(rx,ry), bandwidth=(bw,bw))

# Create a raster file
ras=Raster(Bkde.density, (X(rx), Y((ry))), crs=EPSG(54009))

# Upload the World Shapefile 
World=Shapefile.Table("World.shp") |> DataFrame

plot(World.geometry)

# Crop the raster using the World Shapefile
ras=Rasters.crop(ras;to=World.geometry)

# Mask the raster using the World Shapefile 

ras_masked=Rasters.mask(ras,with=World.geometry)

Plots.plot(ras_masked)

# Define the coords reference system 
crs = ArchGDAL.toWKT(ArchGDAL.importPROJ4("+proj=moll +lon_0=0 +x_0=0 +y_0=0"))
ras_masked=setcrs(ras_masked,crs)

mean(skipmissing(ras_masked))

Rasters.write("SpatialKernelall.tif", ras_masked)

```

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [December 21, 2023, 12:34pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/28 "2023-12-21T12:34:36Z")

</div>

The standard in gdal/most of the spatial world is for Y to be reversed.

Rasters used to do this for you automatically, but it turned out to be annoying in some edge cases of skew and rotation and it means round-trip read/write returned a different array. So now you need to reverse it yourself if you want it. e.g. with `reverse(raster, Y)` afterwards, or just make the Y axis backwards from the start.

For missing values, Julia doesn’t have the concept of `NA` and GDAL/tif doesn’t either. It has either sentinel values (the most common) or a separate bitmask layer for in some formats.

Rasters only uses a sentinel value: `missingval` which you can set manually in `replace_missing` or just let it choose `typemin(eltype(rast))`. But R should read that fine and convert them to `Na`, its part of the generic gdal interface that R packages also use.

`replace_missing(rast, newmissingval)` should work? I don’t understand what ou mean with:

> I would like to substitute them with NA but I can pass only strings to replace\_missing

There is no `NA` in julia. Do you mean you want to use `NaN` ? That is only possible with `Float64`.

Also, you don’t need to do this conversion, Rasters/GDAL will do it automatically for you in `write`:

```julia
crs = ArchGDAL.toWKT(ArchGDAL.importPROJ4("+proj=moll +lon_0=0 +x_0=0 +y_0=0"))
ras_masked=setcrs(ras_masked,crs)

```

You can just do `crs = ProjString("+proj=moll +lon_0=0 +x_0=0 +y_0=0")`.

---

<div class="post-metadata">

**Author:** ![anjelinejeline](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/anjelinejeline/32/205552_2.png) [@anjelinejeline](https://discourse.julialang.org/u/anjelinejeline)\
**Post date:** [December 21, 2023, 1:00pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/29 "2023-12-21T13:00:18Z")

</div>

Hi @Raf you are right and R automatically consider them as Na  
Thank you so much for your support  
PS I agree, an extension of KernelDensity to Rasters may make things easier 🙂

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [December 21, 2023, 1:07pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/30 "2023-12-21T13:07:17Z")

</div>

Since you have learned how it works now, feel free to add it in a PR!

Basically it will just be your code in a `kerneldensity(data; to, res, size, crs, missingval)` function like most other Rasters methods. `data` could be any Table.jl compatible table or GeoInterface.jl compatible `PointTrait` or `MultiPointTrait` or vector of them.

We would add it in an extension based on KernelDensity.jl so we don’t add the dependencies for everyone. There are other extensions in Rasters.jl you can copy, (like for gdal, with e.g. `resample`).

If you get the bones of it in place and working I can help get it finished.

[Previous page](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707.md?page=1)