Create a custom Kernel density function with 3 methods

[anjelinejeline]

Hello! I am trying to create my own Kernel density function to estimate the density of a set of spatial points using the KernelDensity.jl. I would like to define 3 methods:

1. Estimate bidimensional kernel density for spatial points process with all necessary settings specified by the user
2. Estimate biodimensional kernel density for spatial points process specifying a shapefile.table from which the extent and crs are taken from
3. Estimate biodimensional kernel density for spatial points process specifying a raster which will be used for the resolution, extent and crs of the output.

using KernelDensity
using Shapefile
using GeoFormatTypes, Rasters, ArchGDAL
using GeoInterface, Extents

# Methods

# 1. Estimate bidimensional kernel density for spatial points process with all necessary settings specified by the user
# data= array of two dimensions where the first dimension is the x coordinate and the second dimension is the y coordinate
# kerne_dist=the distributional family from Distributions.jl to use as the kernel (default = Normal)
# bw= bandwidth to be passed to the kernel density 
# res= resolution of the output raster 
# ext= Extents.Extent, NamedTuple of tuples holding the lower and upper bounds for each dimension of a object. See Extents.jl
# crs= CRS GeoFormatTypes.jl

rasterKDE(data::Array{Real,2},  kernel_dist::AbstractString, bw::Real, res::Real, ext::Extents.Extent, crs::GeoFormatTypes.CRS)=
rasterKDE(data, kernel_dist, bw, res, ext, crs)

# 2. Estimate biodimensional kernel density for spatial points process specifying a shapefile.table 
# extent and crs are taken from the shapefile shp
rasterKDE(data::Array{Real,2},kernel_dist::AbstractString, bw::Real, shp::Shapefile.Table ,res::Real)=
rasterKDE(data, kernel_dist,bw,shp, res)

# 3. Estimate biodimensional kernel density for spatial points process specifying a raster 
# resolution, extent and crs are taken from the raster ras

rasterKDE(data::Array{Real,2},kernel_dist::AbstractString, bw::Real, ras::AbstractRaster)=
rasterKDE(data, kernel_dist,bw,ras)

# Define function for method 1 
function rasterKDE(data::Array{Real,2}, kernel_dist::AbstractString, bw::Real, res::Real, ext::Extents.Extent, crs::GeoFormatTypes.CRS)

    lengthx=round(Int,(ext.X[2] + abs(ext.X[1])) / res)
    lengthy=round(Int,(ext.Y[2] + abs(ext.Y[1])) / res)

    rx = KernelDensity.kde_range((ext.X[1], ext.X[2]), lengthx)
    ry = KernelDensity.kde_range((ext.Y[1], ext.Y[2]), lengthy)

    # Create a bivariate kernel density estimator
    # If the distribution is not specified use the default (normal distribution)
    if isnothing(kerne_dist)
        Bkde = KernelDensity.kde(data,(rx,ry), bandwidth=(bw,bw))
    else
        Bkde = KernelDensity.kde(data,(rx,ry), bandwidth=(bw,bw), kernel=kernel_dist)
    end
    # Create a raster file
    ras=Rasters.Raster(Bkde.density, (X(rx), Y(reverse(ry))), crs=crs)
    return ras
end

 # Define function for method 2 
 function rasterKDE(data::Array{Real,2},kernel_dist::AbstractString, bw::Real, shp::Shapefile.Table ,res::Real)

    # Retrieve the shapefile extension 
    ext=GeoInterface.bbox(shp)
    crs=GeoInterface.crs(shp)

    lengthx=round(Int,(ext.X[2] + abs(ext.X[1])) / res)
    lengthy=round(Int,(ext.Y[2] + abs(ext.Y[1])) / res)

    rx = KernelDensity.kde_range((ext.X[1], ext.X[2]), lengthx)
    ry = KernelDensity.kde_range((ext.Y[1], ext.Y[2]), lengthy)

    # Create a bivariate kernel density estimator
    # If the distribution is not specified use the default (normal distribution)
    if isnothing(kerne_dist)
        Bkde = KernelDensity.kde(data,(rx,ry), bandwidth=(bw,bw))
    else
        Bkde = KernelDensity.kde(data,(rx,ry), bandwidth=(bw,bw), kernel=kernel_dist)
    end
    # Create a raster file
    ras=Rasters.Raster(Bkde.density, (X(rx), Y(reverse(ry))), crs=crs)
    return ras
   end 

   # Define function for method 3 
function rasterKDE(data::Array{Real,2},kernel_dist::AbstractString, bw::Real, ras::AbstractRaster)

    # Retrieve the shapefile extension 
    ext=GeoInterface.bbox(ras)  # Retrieve the extent
    crs=Rasters.crs(ras)  # Retrive the crs
    res=Rasters.res(ras)  # Retrieve the resolution 

    lengthx=round(Int,(ext.X[2] + abs(ext.X[1])) / res)
    lengthy=round(Int,(ext.Y[2] + abs(ext.Y[1])) / res)

    rx = KernelDensity.kde_range((ext.X[1], ext.X[2]), lengthx)
    ry = KernelDensity.kde_range((ext.Y[1], ext.Y[2]), lengthy)

    # Create a bivariate kernel density estimator
    # If the distribution is not specified use the default (normal distribution)
    if isnothing(kerne_dist)
        Bkde = KernelDensity.kde(data,(rx,ry), bandwidth=(bw,bw))
    else
        Bkde = KernelDensity.kde(data,(rx,ry), bandwidth=(bw,bw), kernel=kernel_dist)
    end
    # Create a raster file
    ras=Rasters.Raster(Bkde.density, (X(rx), Y(reverse(ry))), crs=crs)
    return ras
   end 

A couple of comments about my code:

• method 3 won't work because for some reasons the function Rasters.crs and Rasters.res do not work with my raster (you can find it here https://drive.google.com/file/d/1R1HXv4Uh3NAO4tqPxr6krRySLX8lP7VU/view?usp=sharing ) See the difference in R (where the resolution and crs are correctly retrieved) and Julia (what am I doing wrong ??)

R

Julia

• method1: I tried out the code with the files attached

• but it does not work … I believe that it may be due to a very naive error and some misunderstanding I have about the julia language… Indeed it may be due to the type of data I passed as arguments when I defined the function… the error I got is related to the method … when I defined the different arguments I thought the following:

• data::Array{Real,2} should work with every type of 2D Array/Matrix both float and int type

• kernel_dist::AbstractStringt is a string as per the KernelDensity.jl and it should be set to the default normal distribution if it is not specified

• bw::Real to be flexible and work with both integer and float types

• res::Real to be flexible and work with both integer and float types

• ext::Extents.Extent

• crs::GeoFormatTypes.CRS to be flexible and work with all types of geformat types (given that we cannot retrieve the epsg always)

# Try this out 
# Method 1 

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

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

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

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

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

ext = Extent(X = (-17601618, 17601617), Y = (-9018991, 8751339))
crs=GeoInterface.crs(World)

hope=rasterKDE(data=coordinates, bw=30000, res=1000, ext=ext, crs=crs)

# This throw the error MethodError: no method matching rasterKDE(; data::Matrix{Float64}, bw::Int64, res::Int64, ext::Extent{(:X, :Y),  Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, crs::ESRIWellKnownText{GeoFormatTypes.CRS})

• method 2: I have not tried method 2 as this will throw a similar error of method 1

Can you help me understand what I am doing wrong? I am really new to Julia and I would like to understand the reasons of these errors,

Thanks

Angela df.csv WorldShape.zip

[rafaqz]

Thanks for moving this to Github!

Some things that makes issues easier:

1. no screenshots, always use code.
2. specify that its julia code using ```juliia at the start of the code bloc, then it gets highlighted ( I editied it already).
3. Always include complete call stacks with your error messages - the whole thing from top to bottom. That's where the information about the problem is. Also make errors ```julia
4. Minimum Working Example / MWE is a key concept - try to show your problem with as little code as possible so its fast to read. (I/devs in general have like 15 minutes for this)
5. Try to share smaller files as examples whenever you can, 500mb is a lot
6. Download the files inline in julia using Downloads.download so there is s script that just works I can cut and paste and run. (Again- think about what you want me to be doing in the available 15mins = mostly I'm downloading and unzipping here rather than looking at errors)

And some observations so far: There is no Rasters.res function! that's not a thing. Maybe it should be.

You can do:

res = map(step, dims(raster))

And Rasters.crs does work on your file:

julia> rast = Raster("/home/raf/Downloads/raster.tif")^C

julia> Rasters.crs(rast)
WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\"unknown\",GEOGCS[\"GCS_unknown\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",
\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Mollweide\"],PARAMETER[\"central_meridian\",0],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\"
,1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]")

The other things I need to see the call stack for.

[anjelinejeline]

Hi @rafaqz sorry but I am not able to make the Download function work ..probably because of some constraints to external files I have ( running the code on the VM of my Institution) How can I provide you with the call stack?

Is that what you are looking for?

MethodError: no method matching rasterKDE(; data::Matrix{Float64}, bw::Int64, res::Int64, ext::Extent{(:X, :Y), Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, crs::ESRIWellKnownText{GeoFormatTypes.CRS})

Closest candidates are:
  rasterKDE(::Matrix{Real}, ::AbstractString, ::Real, ::Real, ::Extent, ::GeoFormatTypes.CRS) got unsupported keyword arguments "data", "bw", "res", "ext", "crs"
   @ Main In[3]:11
  rasterKDE(::Matrix{Real}, ::AbstractString, ::Real, ::Shapefile.Table, ::Real) got unsupported keyword arguments "data", "bw", "res", "ext", "crs"
   @ Main In[12]:16
  rasterKDE(::Matrix{Real}, ::AbstractString, ::Real, ::AbstractRaster) got unsupported keyword arguments "data", "bw", "res", "ext", "crs"
   @ Main In[12]:22
  ...

Stacktrace:
 [1] top-level scope
   @ In[15]:1

I am also attaching the jupyter notebook where I ran the code, here you have the output of each run Rasters.KDE.zip

With regard to the raster "ras" I tried to import it on a different file but I used the function Rasters.read instead of Rasters.Raster Can that be the reason for which it was not recognized as raster?

[rafaqz]

?

using Downloads
Downloads.download(url, filename)

By the "call stack" I mean the stacktrace... the full text off the specific errors that gave you trouble. It lists every function call that led to the error, and from that we can kind of see what happened.

Just past it here - like everything from the function you called to the end of the error output.

Rasters.read doesn't do what you think. Thats just the base read function reading a file. Raster is always the constructor function for a Raster object.

[anjelinejeline]

Hi @rafaqz I added the stacktraces() at the end of the script

julia> include("/storage/panelan/prova.jl")
ERROR: LoadError: MethodError: no method matching rasterKDE(; data::Matrix{Float64}, bw::Int64, res::Int64, ext::Extent{(:X, :Y), Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, crs::ESRIWellKnownText{GeoFormatTypes.CRS})

Closest candidates are:
  rasterKDE(::Matrix{Real}, ::AbstractString, ::Real, ::Shapefile.Table, ::Real) got unsupported keyword arguments "data", "bw", "res", "ext", "crs"
   @ Main /storage/panelan/prova.jl:54
  rasterKDE(::Matrix{Real}, ::AbstractString, ::Real, ::AbstractRaster) got unsupported keyword arguments "data", "bw", "res", "ext", "crs"
   @ Main /storage/panelan/prova.jl:80
  rasterKDE(::Matrix{Real}, ::AbstractString, ::Real, ::Real, ::Extent, ::GeoFormatTypes.CRS) got unsupported keyword arguments "data", "bw", "res", "ext", "crs"
   @ Main /storage/panelan/prova.jl:33

Stacktrace:
 [1] top-level scope
   @ /storage/panelan/prova.jl:128
 [2] include(fname::String)
   @ Base.MainInclude ./client.jl:478
 [3] top-level scope
   @ REPL[8]:1
in expression starting at /storage/panelan/prova.jl:128

prova.zip

I am sorry but I am really not able to use the download function as I can't access my google drive from this virtual machine and I cannot work it locally at the moment

[rafaqz]

Try to read the message, it tells you the problem. Your methods dont match your inputs. If you are writing your own functions you need to learn these outputs

On Sun, 24 Dec 2023, 19:56 anjelinejeline, @.***> wrote:

Hi @rafaqz https://github.com/rafaqz I added the stacktraces() at the end of the script

julia> include("/storage/panelan/prova.jl") ERROR: LoadError: MethodError: no method matching rasterKDE(; data::Matrix{Float64}, bw::Int64, res::Int64, ext::Extent{(:X, :Y), Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, crs::ESRIWellKnownText{GeoFormatTypes.CRS})

Closest candidates are: rasterKDE(::Matrix{Real}, ::AbstractString, ::Real, ::Shapefile.Table, ::Real) got unsupported keyword arguments "data", "bw", "res", "ext", "crs" @ Main /storage/panelan/prova.jl:54 rasterKDE(::Matrix{Real}, ::AbstractString, ::Real, ::AbstractRaster) got unsupported keyword arguments "data", "bw", "res", "ext", "crs" @ Main /storage/panelan/prova.jl:80 rasterKDE(::Matrix{Real}, ::AbstractString, ::Real, ::Real, ::Extent, ::GeoFormatTypes.CRS) got unsupported keyword arguments "data", "bw", "res", "ext", "crs" @ Main /storage/panelan/prova.jl:33

Stacktrace: [1] top-level scope @ /storage/panelan/prova.jl:128 [2] include(fname::String) @ Base.MainInclude ./client.jl:478 [3] top-level scope @ REPL[8]:1 in expression starting at /storage/panelan/prova.jl:128

prova.zip https://github.com/rafaqz/Rasters.jl/files/13762622/prova.zip

— Reply to this email directly, view it on GitHub https://github.com/rafaqz/Rasters.jl/issues/586#issuecomment-1868576460, or unsubscribe https://github.com/notifications/unsubscribe-auth/AATKU6KAAUNTEMOIEEBDXETYLB3FTAVCNFSM6AAAAABBBJ46MOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNRYGU3TMNBWGA . You are receiving this because you were mentioned.Message ID: @.***>

[anjelinejeline]

As I suspected but I do not understand the reason .. when I defined the different arguments I thought the following:

• data::Array{Real,2} should work with every type of 2D Array/Matrix both float and int type
• kernel_dist::AbstractStringt is a string as per the KernelDensity.jl and it should be set to the default normal distribution if it is not specified (
• bw::Real to be flexible and work with both integer and float types
• res::Real to be flexible and work with both integer and float types
• ext::Extents.Extent
• crs::GeoFormatTypes.CRS to be flexible and work with all types of geformat types (given that we cannot retrieve the epsg always)

which of this argument I misinterpreted and is causing the issue?

[rafaqz]

Array{Real,2} means it is exactly typed Real ! That means an Array that can hold any mix of values, like Int, Float64 mixed together.

You probably want Array{<:Real,2} so that anything <: Real is accepted.

GeoFormatTypes.CRS is also not what you want. You want GeoFormat or Union{CoordinateReferenceSystemFormat,MixedFormat}. CRS marks something like WellKnownText{CRS} when we know its CRS not Geometry.

Can I suggest you start systematically checking your object types with typeof(obj) before writing dispatch types - guessing types from a package will never work - you need to actually know a type is the exact abstract supertype to dispatch on. supertype(typeof(obj)) shows the abstract type an object inherits from.

Like supertype(EPSG)

julia> using GeoFormatTypes

julia> supertype(EPSG)
CoordinateReferenceSystemFormat

julia> supertype(CoordinateReferenceSystemFormat)
GeoFormat

julia> supertype(WellKnownText)
AbstractWellKnownText

julia> supertype(AbstractWellKnownText)
MixedFormat

julia> supertype(MixedFormat)
GeoFormat

[rafaqz]

You can also check inheritance like this:

julia> Array{Int,2} <: Array{Real,2}
false

julia> Array{Int,2} <: Array{<:Real,2}
true

Here you can see which Array type will work in your case.

[anjelinejeline]

Hi @rafaqz First of all Merry Christmas :) Thanks this was helpful , I just try to look into the supertype of every object in order to use them in the function, I have also deleted the argument specifying the distribution type to keep it simple

using KernelDensity
using Shapefile
using GeoFormatTypes, Rasters, ArchGDAL
using GeoInterface, Extents
using CSV, DataFrames

# Try this out 
# Method 1 

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

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

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

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

#  Upload the World Shapefile 
World=Shapefile.Table("WorldShape/World.shp")

crs=GeoInterface.crs(World)
supertype(typeof(crs))

AbstractWellKnownText{GeoFormatTypes.CRS}

ext = Extent(X = (-17601618, 17601617), Y = (-9018991, 8751339))
supertype(typeof(ext))

Any

supertype(typeof(coordinates))

DenseMatrix{Float64} (alias for DenseArray{Float64, 2})

function rasterKDE(data::DenseMatrix{Float64} ,bw::Signed, res::Signed, ext::Any, crs::AbstractWellKnownText)
    lengthx=round(Int,(ext.X[2] + abs(ext.X[1])) / res)
    lengthy=round(Int,(ext.Y[2] + abs(ext.Y[1])) / res)

    rx = KernelDensity.kde_range((ext.X[1], ext.X[2]), lengthx)
    ry = KernelDensity.kde_range((ext.Y[1], ext.Y[2]), lengthy)

    # Create a bivariate kernel density estimator
    # If the distribution is not specified use the default (normal distribution)
        Bkde = KernelDensity.kde(data,(rx,ry), bandwidth=(bw,bw))
    # Create a raster file
    ras=Rasters.Raster(Bkde.density, (X(rx), Y(reverse(ry))), crs=crs)
    return ras
end
rasterKDE (generic function with 1 method)

but I still got the same error

hope=rasterKDE(data=coordinates, bw=bandwidth, res=res, ext=ext, crs=crs)
MethodError: no method matching rasterKDE(; data::Matrix{Float64}, bw::Int64, res::Int64, ext::Extent{(:X, :Y), Tuple{Tuple{Int64, Int64}, Tuple{Int64, Int64}}}, crs::ESRIWellKnownText{GeoFormatTypes.CRS})

Closest candidates are:
  rasterKDE(::DenseMatrix{Float64}, ::Signed, ::Signed, ::Any, ::AbstractWellKnownText) got unsupported keyword arguments "data", "bw", "res", "ext", "crs"
   @ Main In[8]:1

Stacktrace:
 [1] top-level scope
   @ In[10]:1

[anjelinejeline]

It seems that the issue is due to the unsupported keyword arguments "data", "bw", "res", "ext", "crs" as It was a very naive error I made as I was calling the keyarguments instead of the positions

hope=rasterKDE(data=coordinates, bw=bandwidth, res=res, ext=ext, crs=crs)

while this actually works

hope=rasterKDE(coordinates, 30000, 1000, ext, crs)