add a kerneldensity method

[rafaqz]

trafficstars

@anjelinejeline if you want to try this out you can check out this branch with ] add Rasters#kerneldensity and it should just work. You can try the code in the tests here to get started.

[anjelinejeline]

@rafaqz I have installed the extension as you suggested ] add Rasters#kerneldensity from the REPL

(@v1.9) pkg> add Rasters#kerneldensity
     Cloning git-repo `https://github.com/rafaqz/Rasters.jl.git`
    Updating git-repo `https://github.com/rafaqz/Rasters.jl.git`
   Resolving package versions...
    Updating `~/.julia/environments/v1.9/Project.toml`
  [a3a2b9e3] ~ Rasters v0.10.0 ⇒ v0.10.0 `https://github.com/rafaqz/Rasters.jl.git#kerneldensity`
    Updating `~/.julia/environments/v1.9/Manifest.toml`
  [a3a2b9e3] ~ Rasters v0.10.0 ⇒ v0.10.0 `https://github.com/rafaqz/Rasters.jl.git#kerneldensity`
Precompiling project...
  3 dependencies successfully precompiled in 15 seconds. 275 already precompiled.

(@v1.9) pkg> 

but the function is not recognized

using Rasters

# Load the raster template 
ras=Rasters.Raster("raster.tif")

# 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)

# Compute the kernel 

try1=Rasters.kerneldensity(data=coordinates,to=ras, bandwidth=30000)

Error:

UndefVarError: `kerneldensity` not defined

Stacktrace:
 [1] getproperty(x::Module, f::Symbol)
   @ Base ./Base.jl:31
 [2] top-level scope
   @ In[18]:1`

[rafaqz]

Maybe restart your session? If kerneldensity is not defined then you probably don't have this branch active. (if you added that in a session with Rasters already active the new extension code wont be picked up)

Then you also need to do using KernelDensity as wel as using Rasters

The extension part means it only loads when KernelDensity.jl is also loaded.

First maybe try running the code from the tests so we are on the same page:

using Test, Rasters, KernelDensity, Extents

points = collect(zip(rand(100) .* 10, rand(100)))
# Using just res, detecting extent from points
res_rast = kerneldensity(points; res=(0.1, 0.01))
# This doesn't exactly match, someting is wrong with `_extent2dims`
@test_broken map(step, dims(res_rast)) == (0.1, 0.01)
# Using just size, detecting extent from points
size_rast = kerneldensity(points; size=250)
@test size(size_rast) == (250, 250)
# Using a defined extent with res
to = Extents.Extent(X=(0, 10), Y=(0, 5))
ext_res_rast = kerneldensity(points; to, res=0.05)
@test size(ext_rast) == (200, 100)
@test map(step, dims(ext_rast)) == (0.05, 0.05)
# Using a defined extent with size
to = Extents.Extent(X=(0, 10), Y=(0, 5))
ext_res_rast = kerneldensity(points; to, size=(200, 100))
@test size(ext_rast) == (200, 100)
@test map(step, dims(ext_rast)) == (0.05, 0.05)
# Using a defined DimArray
to = rand(X(0:0.01:10), Y(0:0.001:1))
template_rast = kerneldensity(points; to)
@test size(to) == size(template_rast)
@test map(step, dims(template_rast)) == (0.01, 0.001)

# Test plots
using GLMakie
p = Makie.heatmap(res_rast)
Makie.scatter!(points)

p = Makie.heatmap(size_rast)
Makie.scatter!(points)

p = Makie.heatmap(ext_rast)
Makie.scatter!(points)

p = Makie.heatmap(template_rast)
Makie.scatter!(points)

[anjelinejeline]

Hi @rafaqz I just tried to run the test but something is missing ( ext_rast pheraphs?)

using Test, Rasters, KernelDensity, Extents

points = collect(zip(rand(100) .* 10, rand(100)))
# Using just res, detecting extent from points
res_rast = kerneldensity(points; res=(0.1, 0.01))
# This doesn't exactly match, someting is wrong with `_extent2dims`
@test_broken map(step, dims(res_rast)) == (0.1, 0.01)
# Using just size, detecting extent from points
size_rast = kerneldensity(points; size=250)
@test size(size_rast) == (250, 250)
# Using a defined extent with res
to = Extents.Extent(X=(0, 10), Y=(0, 5))
ext_res_rast = kerneldensity(points; to, res=0.05)
@test size(ext_rast) == (200, 100)
@test map(step, dims(ext_rast)) == (0.05, 0.05)
# Using a defined extent with size
to = Extents.Extent(X=(0, 10), Y=(0, 5))
ext_res_rast = kerneldensity(points; to, size=(200, 100))
@test size(ext_rast) == (200, 100)
@test map(step, dims(ext_rast)) == (0.05, 0.05)
# Using a defined DimArray
to = rand(X(0:0.01:10), Y(0:0.001:1))
template_rast = kerneldensity(points; to)
@test size(to) == size(template_rast)
@test map(step, dims(template_rast)) == (0.01, 0.001)

Error During Test at In[2]:14
  Test threw exception
  Expression: size(ext_rast) == (200, 100)
  UndefVarError: `ext_rast` not defined
  Stacktrace:
   [1] top-level scope
     @ /software/22.04/julia/julia-1.9.3/share/julia/stdlib/v1.9/Test/src/Test.jl:478
   [2] eval
     @ ./boot.jl:370 [inlined]
   [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
     @ Base ./loading.jl:1903
   [4] softscope_include_string(m::Module, code::String, filename::String)
     @ SoftGlobalScope ~/.julia/packages/SoftGlobalScope/u4UzH/src/SoftGlobalScope.jl:65
   [5] execute_request(socket::ZMQ.Socket, msg::IJulia.Msg)
     @ IJulia ~/.julia/packages/IJulia/Vo51o/src/execute_request.jl:67
   [6] #invokelatest#2
     @ ./essentials.jl:819 [inlined]
   [7] invokelatest
     @ ./essentials.jl:816 [inlined]
   [8] eventloop(socket::ZMQ.Socket)
     @ IJulia ~/.julia/packages/IJulia/Vo51o/src/eventloop.jl:8
   [9] (::IJulia.var"#15#18")()
     @ IJulia ./task.jl:514

[rafaqz]

Oh yeah the raster name changed.

[rafaqz]

These are the new fixed tests. You can plot by running specific lines and then the plot commands at the bottom.

using Test, Rasters, KernelDensity, Extents

points = append!(collect(zip(rand(100) .* 10, rand(100))), [(0.0, 0.0), (10.0, 1.0)])

@testset "Using just res, detecting extent from points" begin
    rast = kerneldensity(points; res=(0.1, 0.01))
    @test map(step, dims(rast)) == (0.1, 0.01)
end

@testset "Using just res, detecting extent from points" begin
    rast = kerneldensity(points; size=250)
    @test size(rast) == (250, 250)
end

@testset "Using a defined extent with res" begin
    to = Extents.Extent(X=(0, 10), Y=(0, 5))
    rast = kerneldensity(points; to, res=0.05)
    @test size(rast) == (200, 100)
    @test map(step, dims(rast)) == (0.05, 0.05)
end

@testset "Using a defined extent with size" begin
    to = Extents.Extent(X=(0, 10), Y=(0, 5))
    rast = kerneldensity(points; to, size=(200, 100))
    @test size(rast) == (200, 100)
    @test map(step, dims(rast)) == (0.05, 0.05)
end

@testset "Using a preddefined DimArray" begin
    to = rand(X(0:0.01:10), Y(0:0.001:1))
    rast = kerneldensity(points; to)
    @test size(to) == size(rast)
    @test map(step, dims(rast)) == (0.01, 0.001)
end

# Test plots
# using GLMakie
# p = Makie.heatmap(rast)
# Makie.scatter!(points)

[anjelinejeline]

Hi @rafaqz I have trouble plotting the test as you suggested

 Info: Precompiling GLMakie [e9467ef8-e4e7-5192-8a1a-b1aee30e663a]
Errors encountered while load FileIO.File{FileIO.DataFormat{:PNG}, String}("/home/panelan/.julia/packages/Makie/Qvk4f/assets/icons/icon-128.png").
All errors:
===========================================
Evaluation into the closed module `ImageIO` breaks incremental compilation because the side effects will not be permanent. This is likely due to some other module mutating `ImageIO` with `eval` during precompilation - don't do this.
===========================================
ArgumentError: Package ImageMagick [6218d12a-5da1-5696-b52f-db25d2ecc6d1] is required but does not seem to be installed:
 - Run `Pkg.instantiate()` to install all recorded dependencies.

===========================================

Fatal error:
ERROR: LoadError: Evaluation into the closed module `ImageIO` breaks incremental compilation because the side effects will not be permanent. This is likely due to some other module mutating `ImageIO` with `eval` during precompilation - don't do this.
Stacktrace:
  [1] eval
    @ ./boot.jl:370 [inlined]
  [2] (::ImageIO.var"#1#2"{Symbol})()
    @ ImageIO ~/.julia/packages/ImageIO/bBdyA/src/ImageIO.jl:19
  [3] lock(f::ImageIO.var"#1#2"{Symbol}, l::ReentrantLock)
    @ Base ./lock.jl:229
  [4] checked_import
    @ ~/.julia/packages/ImageIO/bBdyA/src/ImageIO.jl:8 [inlined]
  [5] load(f::FileIO.File{FileIO.DataFormat{:PNG}, String}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ImageIO ~/.julia/packages/ImageIO/bBdyA/src/ImageIO.jl:27
  [6] load(f::FileIO.File{FileIO.DataFormat{:PNG}, String})
    @ ImageIO ~/.julia/packages/ImageIO/bBdyA/src/ImageIO.jl:26
  [7] #invokelatest#2
    @ ./essentials.jl:819 [inlined]
  [8] invokelatest
    @ ./essentials.jl:816 [inlined]
  [9] action(::Symbol, ::Vector{Union{Base.PkgId, Module}}, ::FileIO.Formatted; options::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ FileIO ~/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:219
 [10] action
    @ ~/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:196 [inlined]
 [11] action(::Symbol, ::Vector{Union{Base.PkgId, Module}}, ::Symbol, ::String; options::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ FileIO ~/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:185
 [12] action
    @ ~/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:185 [inlined]
 [13] load(::String; options::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ FileIO ~/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:113
 [14] load
    @ ~/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:109 [inlined]
 [15] _broadcast_getindex_evalf
    @ ./broadcast.jl:683 [inlined]
 [16] _broadcast_getindex
    @ ./broadcast.jl:656 [inlined]
 [17] getindex
    @ ./broadcast.jl:610 [inlined]
 [18] copy
    @ ./broadcast.jl:912 [inlined]
 [19] materialize
    @ ./broadcast.jl:873 [inlined]
 [20] icon()
    @ Makie ~/.julia/packages/Makie/Qvk4f/src/Makie.jl:313
 [21] empty_screen(debugging::Bool; reuse::Bool)
    @ GLMakie ~/.julia/packages/GLMakie/WlUOp/src/screen.jl:260
 [22] empty_screen
    @ ~/.julia/packages/GLMakie/WlUOp/src/screen.jl:222 [inlined]
 [23] singleton_screen(debugging::Bool)
    @ GLMakie ~/.julia/packages/GLMakie/WlUOp/src/screen.jl:329
 [24] macro expansion
    @ ~/.julia/packages/GLMakie/WlUOp/src/precompiles.jl:21 [inlined]
 [25] macro expansion
    @ ~/.julia/packages/PrecompileTools/kmH5L/src/workloads.jl:78 [inlined]
 [26] macro expansion
    @ ~/.julia/packages/GLMakie/WlUOp/src/precompiles.jl:18 [inlined]
 [27] macro expansion
    @ ~/.julia/packages/PrecompileTools/kmH5L/src/workloads.jl:140 [inlined]
 [28] top-level scope
    @ ~/.julia/packages/GLMakie/WlUOp/src/precompiles.jl:16
Stacktrace:
  [1] handle_error(e::ErrorException, q::Base.PkgId, bt::Vector{Union{Ptr{Nothing}, Base.InterpreterIP}})
    @ FileIO ~/.julia/packages/FileIO/BE7iZ/src/error_handling.jl:61
  [2] handle_exceptions(exceptions::Vector{Tuple{Any, Union{Base.PkgId, Module}, Vector}}, action::String)
    @ FileIO ~/.julia/packages/FileIO/BE7iZ/src/error_handling.jl:56
  [3] action(::Symbol, ::Vector{Union{Base.PkgId, Module}}, ::FileIO.Formatted; options::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ FileIO ~/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:228
  [4] action
    @ ~/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:196 [inlined]
  [5] action(::Symbol, ::Vector{Union{Base.PkgId, Module}}, ::Symbol, ::String; options::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ FileIO ~/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:185
  [6] action
    @ ~/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:185 [inlined]
  [7] load(::String; options::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ FileIO ~/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:113
  [8] load
    @ ~/.julia/packages/FileIO/BE7iZ/src/loadsave.jl:109 [inlined]
  [9] _broadcast_getindex_evalf
    @ ./broadcast.jl:683 [inlined]
 [10] _broadcast_getindex
    @ ./broadcast.jl:656 [inlined]
 [11] getindex
    @ ./broadcast.jl:610 [inlined]
 [12] copy
    @ ./broadcast.jl:912 [inlined]
 [13] materialize
    @ ./broadcast.jl:873 [inlined]
 [14] icon()
    @ Makie ~/.julia/packages/Makie/Qvk4f/src/Makie.jl:313
 [15] empty_screen(debugging::Bool; reuse::Bool)
    @ GLMakie ~/.julia/packages/GLMakie/WlUOp/src/screen.jl:260
 [16] empty_screen
    @ ~/.julia/packages/GLMakie/WlUOp/src/screen.jl:222 [inlined]
 [17] singleton_screen(debugging::Bool)
    @ GLMakie ~/.julia/packages/GLMakie/WlUOp/src/screen.jl:329
 [18] macro expansion
    @ ~/.julia/packages/GLMakie/WlUOp/src/precompiles.jl:21 [inlined]
 [19] macro expansion
    @ ~/.julia/packages/PrecompileTools/kmH5L/src/workloads.jl:78 [inlined]
 [20] macro expansion
    @ ~/.julia/packages/GLMakie/WlUOp/src/precompiles.jl:18 [inlined]
 [21] macro expansion
    @ ~/.julia/packages/PrecompileTools/kmH5L/src/workloads.jl:140 [inlined]
 [22] top-level scope
    @ ~/.julia/packages/GLMakie/WlUOp/src/precompiles.jl:16
 [23] include(mod::Module, _path::String)
    @ Base ./Base.jl:457
 [24] include(x::String)
    @ GLMakie ~/.julia/packages/GLMakie/WlUOp/src/GLMakie.jl:1
 [25] top-level scope
    @ ~/.julia/packages/GLMakie/WlUOp/src/GLMakie.jl:79
 [26] include
    @ ./Base.jl:457 [inlined]
 [27] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt128}}, source::String)
    @ Base ./loading.jl:2049
 [28] top-level scope
    @ stdin:3
in expression starting at /home/panelan/.julia/packages/GLMakie/WlUOp/src/precompiles.jl:15
in expression starting at /home/panelan/.julia/packages/GLMakie/WlUOp/src/GLMakie.jl:1
in expression starting at stdin:3
ERROR: LoadError: Failed to precompile GLMakie [e9467ef8-e4e7-5192-8a1a-b1aee30e663a] to "/home/panelan/.julia/compiled/v1.9/GLMakie/jl_6zJ265".
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, keep_loaded_modules::Bool)
    @ Base ./loading.jl:2300
  [3] compilecache
    @ ./loading.jl:2167 [inlined]
  [4] _require(pkg::Base.PkgId, env::String)
    @ Base ./loading.jl:1805
  [5] _require_prelocked(uuidkey::Base.PkgId, env::String)
    @ Base ./loading.jl:1660
  [6] macro expansion
    @ ./loading.jl:1648 [inlined]
  [7] macro expansion
    @ ./lock.jl:267 [inlined]
  [8] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:1611
  [9] include(fname::String)
    @ Base.MainInclude ./client.jl:478
 [10] top-level scope
    @ REPL[31]:1

[rafaqz]

Did you install ImageMagick like the message says? Do that and try again (there are a few possible image backends so it doesn't install automatically, unfortunately)

[anjelinejeline]

Yes, thanks this issue can be fixed if ImageMagick is installed .. but I am now getting an error related to rast which is not defined apparently

[rafaqz]

You need to run the specific line individually not in the test block just copy and paste it out

The run the plot code for whatever youre interested in

[anjelinejeline]

Hi rafa .. it works, but I cannot plot it .. this code returns a Scatter{Tuple{Vector{Point{2, Float32}}}}

using GLMakie, Makie
p = Makie.heatmap(rast)
Makie.scatter!(points)

[rafaqz]

Not sure how you are running these things? It should call display automatically if you do that in the REPL.

But you can call display manually on the output if thats what you want. display(p)

(Are you always running things as whole scripts? the best way to work with Julia is in a live REPL session)

[anjelinejeline]

Hi rafa .. it works!! thanks .. the output of the function also seems OKAY