Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Create a custom Kernel density function with 3 methods #586

Closed
anjelinejeline opened this issue Dec 24, 2023 · 10 comments
Closed

Create a custom Kernel density function with 3 methods #586

anjelinejeline opened this issue Dec 24, 2023 · 10 comments

Comments

@anjelinejeline
Copy link

anjelinejeline commented Dec 24, 2023 •
edited by rafaqz
Loading

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:

R
R

Julia
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
Copy link
Owner

rafaqz commented Dec 24, 2023 •
edited
Loading

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
Copy link
Author

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
Copy link
Owner

rafaqz commented Dec 24, 2023 •
edited
Loading

?

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
Copy link
Author

anjelinejeline commented Dec 24, 2023 •
edited
Loading

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
Copy link
Owner

rafaqz commented Dec 24, 2023 via email

@anjelinejeline
Copy link
Author

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
Copy link
Owner

rafaqz commented Dec 24, 2023 •
edited
Loading

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
Copy link
Owner

rafaqz commented Dec 24, 2023

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
Copy link
Author

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
Copy link
Author

anjelinejeline commented Dec 25, 2023 •
edited
Loading

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)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@rafaqz @anjelinejeline