# 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:** 20\
**Page:** 1

<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 16, 2023, 11:08am UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/1 "2023-12-16T11:08:38Z")

</div>

Hi,  
I am really new to Julia,  
I would like to compute a spatial kernel with a resolution of 1000x1000 metres and then save the results as a raster file (“+proj=moll +lon\_0=0 +x\_0=0 +y\_0=0”).

I gave it a try with the KernelDensity package and then trying to interpolate the values but I got stuck as I do not understand how to properly use the pdf function.

Can anyone help?

Thank you for your help

This is my code till I got stuck

```julia
using CSV, DataFrames, KernelDensity, Proj, ArchGDAL, Shapefile , Colors, Distributions, Interpolations

# Load the CSV file
df = CSV.read("df.csv", DataFrame)

# Read the world shapefile
world_shapefile = Shapefile.read("WorldShape/World.shp")

# Discard the first column
df = select(df, Not(1))

# Filter by year 
df = filter(row -> 1975 <= row.outbreak_start <= 2020, df)

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

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

# Create a bivariate kernel density estimator
Bkde = kde(coordinates)

```

Here is where I face issues

```julia
# Set resolution
resolution = 1000

# Generate grid coordinates
x_grid = range(minimum(x_values), maximum(x_values), length=resolution)
y_grid = range(minimum(y_values), maximum(y_values), length=resolution)

# Evaluate KDE on the grid
kde_values = pdf(Bkde, hcat(vec(x_grid), vec(y_grid)))

```

---

<div class="post-metadata">

**Author:** ![mlkrock](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mlkrock/32/31211_2.png) [@mlkrock](https://discourse.julialang.org/u/mlkrock)\
**Post date:** [December 18, 2023, 6:15pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/2 "2023-12-18T18:15:00Z")

</div>

Hello @anjelinejeline I think this is a simple fix, just need to change your last line to  
`kde_values = pdf(Bkde, x_grid, y_grid)`

---

<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 18, 2023, 10:21pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/3 "2023-12-18T22:21:38Z")

</div>

Thank you @mlkrock it works !!  
Do you know how I can turn it into a raster Geotiff and assign the right crs ?

Thank you

---

<div class="post-metadata">

**Author:** ![Fliks](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fliks/32/2494_2.png) [@Fliks](https://discourse.julialang.org/u/Fliks)\
**Post date:** [December 18, 2023, 11:13pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/4 "2023-12-18T23:13:19Z")

</div>

You could use the Rasters package to handle the geotiff part. You can have to see whether it is necessary to reverse the y coordinates that would depend on the values in it.

```julia
using GeoFormatTypes, Rasters, ArchGDAL
julia> ras = Raster(kde_values, (X(x_grid), Y(reverse(y_grid))), crs=EPSG(4326))
1000×1000 Raster{Float64,2} with dimensions: 
  X Projected{Float64} 1.0:0.009009009009009009:10.0 ForwardOrdered Regular Points crs: EPSG,
  Y Projected{Float64} 30.0:-0.009009009009009009:21.0 ReverseOrdered Regular Points crs: EPSG
extent: Extent(X = (1.0, 10.0), Y = (21.0, 30.0))
missingval: missing
crs: EPSG:4326

julia> write(tempname()*".tif", ras)
┌ Warning: `missing` cant be written with gdal, missinval for `Float64` of `-Inf` used instead
└ @ Rasters ~/.julia/packages/Rasters/PnKS7/src/utils.jl:29
"/tmp/jl_CJkPkl8F4T.tif"

```

---

<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, 2:52pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/6 "2023-12-20T14:52:27Z")

</div>

Hi @Fliks  
Thank you  
I was able to create a kernel and interpolate it on a grid of the World but the resolution of my output does not match the resolution I initially set for the kernel.

Can you help me get a raster of 1kmx1km resolution?

Thank you

```julia
bw=30000::Int64
# Create a bivariate kernel density estimator

Bkde = kde(coordinates, bandwidth=(bw,bw))

# Set resolution
resolution = 1000

# Generate grid coordinates considering the bounding of the World ESRI:54009 Mollweide
x_grid = range(-17601618 , 17601617, length=resolution)
y_grid = range(-9018991, 8751339, length=resolution)

# Interpolate the values on the World grid 
kde_values = pdf(Bkde, x_grid, y_grid)

# Create a raster file
ras=Raster(kde_values, (X(x_grid), Y(reverse(y_grid))), crs=EPSG(54009))

```

---

<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, 4:00pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/7 "2023-12-20T16:00:46Z")

</div>

> [@anjelinejeline](#):
>
> Can you help me get a raster of 1kmx1km resolution?

Resolution should be determined by data, not by _wish_ but ofc one can do that.  
Hard to advise more without a reproducible example.

---

<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, 5:41pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/8 "2023-12-20T17:41:25Z")

</div>

Hi @joa-quim joaquin  
My objective is to create a kernel density of a set of spatial points.  
The coordinates of my points are in ESRI:54009 Mollweide  
The bandwidth I want to use is 3000 metres  
I want to interpolate the density on a grid of 1x1km resolution covering the entire world

What I have done so far

I am happy to share the df object but I cannot upload it here as csv are not supported

```julia
using CSV, DataFrames, KernelDensity, GeoFormatTypes, Rasters, ArchGDAL, Distributions, Interpolations
using Plots, Shapefile

# Load the CSV file
df = CSV.read("df.csv", DataFrame)

# Discard the first column
df = select(df, Not(1))

# Filter by year 
df = filter(row -> 1960 <= row.outbreak_start <= 2020, df)

# 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
# Create a bivariate kernel density estimator

Bkde = kde(coordinates, bandwidth=(bw,bw))

# Set resolution
resolution = 1000

# Generate grid coordinates considering the bounding of the World ESRI:54009 Mollweide
x_grid = range(-17601618 , 17601617, length=resolution)
y_grid = range(-9018991, 8751339, length=resolution)

# Interpolate the values on the World grid 
kde_values = pdf(Bkde, x_grid, y_grid)

# Create a raster file
ras=Raster(kde_values, (X(x_grid), Y(reverse(y_grid))), crs=EPSG(54009))

plot(ras)

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

plot(World.geometry)

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

# Mask the raster using the World Shapefile 

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

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)

```

---

<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, 5:53pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/9 "2023-12-20T17:53:58Z")

</div>

> [@anjelinejeline](#):
>
> want to interpolate the density on a grid of 1x1km resolution covering the entire world

Mind you that this, in a Float32 array is ~3.5 GB

> [@anjelinejeline](#):
>
> ```julia
> # Load the CSV file
> df = CSV.read("df.csv", DataFrame)
> 
> ```

This is not reproducible.

---

<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, 6:22pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/10 "2023-12-20T18:22:40Z")

</div>

[df\_file.jl](https://discourse.julialang.org/uploads/short-url/lfhqLy3ieMEwwA4qMLe5O8O48H1.jl) (30.1 KB)  
I have used the JLD2 package ([Saving data files in Julia](https://shuvomoy.github.io/blogs/posts/Saving_data_files_julia/)) to share the df with you.  
I am happy to share the data as csv but I cannot upload this format here

---

<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, 7:26pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/11 "2023-12-20T19:26:16Z")

</div>

Sorry, can’t use that data. It errors with

```julia
df = filter(row -> 1960 <= row.outbreak_start <= 2020, df)
ERROR: ArgumentError: column name :outbreak_start not found in the data frame

```

But let me show an example on how that works with GMT.jl and an example from the KernelDensity test suit.  
You also seem to be doing some clipping inside continents, which GMT can do as well without any further data (the [coast](https://www.generic-mapping-tools.org/GMTjl_doc/documentation/modules/coast/) module can do land/ocean clipping

```julia
using GMT, KernelDensity

# Use test example in the KernelDensity package
r = KernelDensity.kde_range((-2.0,2.0), 128);
k12 = kde([1.0 1.0], (r,r), bandwidth=(1,1));

# Create a GMT grid. I', expanding the Mollweide because something is not working with "epsg=54009"
G = mat2grid(Float32.(k12.density), x=k12.x, y=k12.y, proj4="moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs");

# Save as a GeoTIFF file
gmtwrite("k12.tiff", G)

# Show it
viz("k12.tiff")

```

 ![GMTjl_j](https://global.discourse-cdn.com/julialang/original/3X/f/4/f4b58647c24f11d1c655b519c0b5161534dbcac7.png)

---

<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, 7:31pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/12 "2023-12-20T19:31:07Z")

</div>

> [@joa-quim](#):
>
> `KernelDensity.kde_range((`

Thanks joaquim could you please help me understand the code?  
where did you specify the resolution?  
where should the coordinates be inserted?  
where the extent of the world should be specified?

---

<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, 7:41pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/13 "2023-12-20T19:41:51Z")

</div>

I’m not a _kde_ user but I assumed that the k12.x & k12.y contains the coordinates, which were set with the `range(-2,2,128)`. Don’t do that with `rx = range(-180,180, 43200)` (360 \* 60 \* 2 = 43200) and same for lat if you don’t have a some 7 GB to spear with Float64 array.

Note, my example is on how to create a GeoTIFF from the kde result. How to create it as you want is another issue.

---

<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, 7:44pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/14 "2023-12-20T19:44:57Z")

</div>

HI joa-quim  
I was also able to create a kernel with the code I posted  
My issue is the resolution … I want a geotiff of 1kmx1km  
Can you help me with that?  
Where did you specify the resolution of your raster in the code?

---

<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, 7:56pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/15 "2023-12-20T19:56:25Z")

</div>

But I answered that above. If you insist, try

```julia
rx = KernelDensity.kde_range((-180.0,180.0), 43200);
ry = KernelDensity.kde_range((-90.0,90.0), 21600);

```

---

<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, 8:05pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/16 "2023-12-20T20:05:08Z")

</div>

Sorry, your coordinates are Mollweid meters so should not use -180,180 but that in Mollweid meters.

---

<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, 8:13pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/17 "2023-12-20T20:13:26Z")

</div>

Sorry it’s still not clear to me …  
what does the kde\_range define ?  
Sorry but I can’t find a proper documentation here [Readme · KernelDensity.jl](https://docs.juliahub.com/KernelDensity/4QyGx/0.6.3/)

My coordinates are in metres the bounding of the Mollweide are  
x: -17601618 , 17601617  
y: -9018991, 8751339

Should I write something like this?

```julia
# Load the CSV file
df = CSV.read("df.csv", DataFrame)

# Discard the first column
df = select(df, Not(1))

# Filter by year 
df = filter(row -> 1960 <= row.outbreak_start <= 2020, df)

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

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

# Bandwidth 
bw=30000::Int64

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

rx = KernelDensity.kde_range((-17601618 , 17601617), 1000);
ry = KernelDensity.kde_range(( -9018991, 8751339), 1000);

# Create a bivariate kernel density estimator

Bkde = kde(coordinates,(rx,ry), bandwidth=(bw,bw))

```

---

<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, 8:19pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/18 "2023-12-20T20:19:45Z")

</div>

```julia
rx = KernelDensity.kde_range((-17601618 , 17601617), 1000);

```

No. length must be

julia\> round(Int,(17601617 + 17601618) / 1000)  
35203

so

```julia
rx = KernelDensity.kde_range((-17601618 , 17601617), 35202);

```

Do you see now why I say its **big**

---

<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, 8:30pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/19 "2023-12-20T20:30:33Z")

</div>

Thank you I will try right away …  
I know that it’s big that’s why I am trying to move from R to julia 😅 …  
Can you explain me how I can access the documentation of the kde\_range function within the KernelDensity package?  
I am struggling to understand how to properly find guidance on the use of the functions within each julia package

---

<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, 8:35pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/20 "2023-12-20T20:35:38Z")

</div>

> [@anjelinejeline](#):
>
> I know that it’s big that’s why I am trying to move from R to julia 😅 …

A further problem is that calculations return a Float64 array, which is probably not needed. Maybe someone knows how to make it in Float32.

> [@anjelinejeline](#):
>
> Can you explain me how I can access the documentation of the kde\_range function within the KernelDensity package?

Sorry, mentioned above. I’m no KernelDensity user. Just read a bit of the code and guessed the rest.

---

<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:00pm UTC](https://discourse.julialang.org/t/spatial-kernel-in-julia/107707/21 "2023-12-20T21:00:04Z")

</div>

> [@anjelinejeline](#):
>
> `using CSV, DataFrames, KernelDensity,`

Thanks anyway  
With regard to the package GMT the function mat2grid does not work in my case

G = mat2grid(Float32.(Bkde.density), x=Bkde.x, y=Bkde.y, proj4=“+proj=moll +lon\_0=0 +x\_0=0 +y\_0=0”)

 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/c/ec833a6e29c4ac9727cdacd5cfb8b777b1e8d8cf.png)

How can I fix it?

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