From 86806bbe4ab8bfa785ce885250db6259168b9a6d Mon Sep 17 00:00:00 2001 From: rafaqz Date: Fri, 22 Dec 2023 23:07:34 +0100 Subject: [PATCH 1/5] add kerneldensity method --- Project.toml | 6 +- .../RastersKernelDensityExt.jl | 82 +++++++++++++++++++ src/Rasters.jl | 2 +- src/extensions.jl | 35 ++++++++ src/methods/rasterize.jl | 1 + src/methods/shared_docstrings.jl | 6 ++ src/utils.jl | 2 +- test/runtests.jl | 2 +- 8 files changed, 132 insertions(+), 4 deletions(-) create mode 100644 ext/RastersKernelDensityExt/RastersKernelDensityExt.jl diff --git a/Project.toml b/Project.toml index 1173f1e1c..34e9264cc 100644 --- a/Project.toml +++ b/Project.toml @@ -28,6 +28,7 @@ Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3" CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" +KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36" @@ -36,6 +37,7 @@ RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36" RastersArchGDALExt = "ArchGDAL" RastersCoordinateTransformationsExt = "CoordinateTransformations" RastersHDF5Ext = "HDF5" +RastersKernelDensityExt = "KernelDensity" RastersMakieExt = "Makie" RastersNCDatasetsExt = "NCDatasets" RastersRasterDataSourcesExt = "RasterDataSources" @@ -54,6 +56,7 @@ Flatten = "0.4" GeoFormatTypes = "0.4" GeoInterface = "1" HDF5 = "0.14, 0.15, 0.16" +KernelDensity = "0.6" Makie = "0.19, 0.20" Missings = "0.4, 1" NCDatasets = "0.10, 0.11, 0.12, 0.13" @@ -74,6 +77,7 @@ CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" +KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" @@ -84,4 +88,4 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeometryBasics", "HDF5", "NCDatasets", "Plots", "RasterDataSources", "SafeTestsets", "Shapefile", "Statistics", "Test"] +test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeometryBasics", "HDF5", "KernelDensity", "NCDatasets", "Plots", "RasterDataSources", "SafeTestsets", "Shapefile", "Statistics", "Test"] diff --git a/ext/RastersKernelDensityExt/RastersKernelDensityExt.jl b/ext/RastersKernelDensityExt/RastersKernelDensityExt.jl new file mode 100644 index 000000000..3085bff37 --- /dev/null +++ b/ext/RastersKernelDensityExt/RastersKernelDensityExt.jl @@ -0,0 +1,82 @@ +module RastersKernelDensityExt + +using Extents, + GeoInterface, + KernelDensity, + Rasters + +using Rasters: Tables, maybeshiftlocus, Center, lookup, _extent2dims + +const GI = GeoInterface + +const VectorTuple = Tuple{<:AbstractVector{<:Number},<:AbstractVector{<:Number}} + +function Rasters.kerneldensity(geom; geomcolumn=nothing, kw...) + if Tables.istable(typeof(geom)) + if geomcolumn isa Tuple + # Use the columns directly and skip _to_vectors + data = (Tables.getcolumn(geom, geomcolumn[1]), Tables.getcolumn(geom, geomcolumn[2])) + return _kerneldensity(data; kw...) + elseif isnothing(geomcolumn) + geomcolumn = GI.geometrycolumns(geom)[1] + end + geom = GI.getcolumn(geom, geomcolumn) + end + return _kerneldensity(_to_vectors(geom); kw...) +end +# Also accept the regular KernelDensity.jl data formats +Rasters.kerneldensity(data::VectorTuple; kw...) = _kerneldensity(data; kw...) +Rasters.kerneldensity(data::AbstractMatrix{<:Number}; kw...) = + _kerneldensity((data[:, 1], data[:, 2]); kw...) + +function _kerneldensity(data::VectorTuple; + to=nothing, size=nothing, res=nothing, crs=nothing, kw... +) + # Get the extent from `data` if there is no `to` keyword + if isnothing(to) + x, y = map(extrema, data) + to = Extents.Extent(; X=x, Y=y) + end + # Convert input keywords to Dimensions + ds = _extent2dims(to; size, res, crs) + + # KernelDensity needs midpoints + midpoint_ranges = map(parent ∘ lookup, maybeshiftlocus(Center(), ds)) + + # Create a bivariate kernel density estimator + bkde = KernelDensity.kde(data, midpoint_ranges; kw...) + + # Return a raster + return Raster(bkde.density, ds) +end + +_to_vectors(geoms) = _to_vectors(GI.trait(geoms), geoms) +function _to_vectors(::Nothing, geoms::AbstractArray) + if all(GI.geomtrait(g) isa PointTrait for g in geoms) + return _fill_vectors(geoms, length(geoms)) + elseif all(GI.geomtrait(g) isa MultiPointTrait for g in geoms) + n = sum(GI.npoint(g) for g in geoms) + points = Iterators.flatten(GI.getpoint(g) for g in geoms) + return _fill_vectors(points, n) + end +end +function _to_vectors(::GI.MultiPointTrait, geom) + n = GI.npoint(geom) + points = Iterators.flatten(GI.getpoint(g) for g in geoms) + return _fill_vectors(points, n) +end +_to_vectors(::PointTrait, geom) = [GI.x(geom), GI.y(geom)] +function _to_vectors(trait, geom) + throw(ArgumentError("$trait not supported: only `MultiPointTrait` `PointTrait` and `MultiPointTrait` are supported")) +end + +function _fill_vectors(points, n) + vec1 = Vector{Float64}(undef, n) + vec2 = Vector{Float64}(undef, n) + for (i, p) in enumerate(points) + vec1[i], vec2[i] = GI.x(p), GI.y(p) + end + return (vec1, vec2) +end + +end diff --git a/src/Rasters.jl b/src/Rasters.jl index 64bc8131c..818edd1bc 100644 --- a/src/Rasters.jl +++ b/src/Rasters.jl @@ -59,7 +59,7 @@ export missingval, boolmask, missingmask, replace_missing, replace_missing!, aggregate, aggregate!, disaggregate, disaggregate!, mask, mask!, resample, warp, zonal, crop, extend, trim, slice, combine, points, classify, classify!, mosaic, mosaic!, extract, rasterize, rasterize!, - coverage, coverage!, setcrs, setmappedcrs, smapseries, cellsize + coverage, coverage!, setcrs, setmappedcrs, smapseries, cellsize, kerneldensity export crs, mappedcrs, mappedindex, mappedbounds, projectedindex, projectedbounds export reproject, convertlookup diff --git a/src/extensions.jl b/src/extensions.jl index 4dff7411e..f725d5529 100644 --- a/src/extensions.jl +++ b/src/extensions.jl @@ -195,6 +195,41 @@ $EXPERIMENTAL """ cellsize(args...; kw...) = throw_extension_error(cellsize, "ArchGDAL", :RastersArchGDALExt, args) + +""" + kerneldensity(geoms; kw...) + +Kalculate the kernel density of the points in `geoms` using KernelDensity.jl. + +Returns a `Raster`. + +## Keywords + +- `geoms`: an iterable of GeoInterface.jl compatible points or multipoints, + or a Tables.jl compatible table with a geometry column. KernelDensity.jl + compatible 2-column `Matrix` or `Tuple` of `Vectors` are also supported. +$TO_KEYWORD +$RES_KEYWORD +$SIZE_KEYWORD +$CRS_KEYWORD +$GEOMCOLUMN_KEYWORD +- `kernel`: a distribution from Distributions.jl, default is `Normal()`. +- `bandwidth` bandwidth to be passed to the kernel density estimator. +- `weights`: weights to be passed to the kernel density estimator. + +## Example + +```julia +using Rasters, KernelDensity +points = zip(rand(Normal(), 100), rand(Normal(), 100)) + +kerneldensity(points; res=0.1, weights=rand(100)) +``` + +$EXPERIMENTAL +""" +kerneldensity(args...; kw...) = throw_extension_error(warp, "KernelDensity", :RastersKernelDensityExt, args) + # Other shared stubs function layerkeys end function smapseries end diff --git a/src/methods/rasterize.jl b/src/methods/rasterize.jl index 916274eb8..7bd81e46c 100644 --- a/src/methods/rasterize.jl +++ b/src/methods/rasterize.jl @@ -382,6 +382,7 @@ These are detected automatically from `data` where possible. $GEOM_KEYWORDS $RASTERIZE_KEYWORDS +$GEOMCOLUMN_KEYWORD $FILENAME_KEYWORD $SUFFIX_KEYWORD diff --git a/src/methods/shared_docstrings.jl b/src/methods/shared_docstrings.jl index f184534c4..fc2c30860 100644 --- a/src/methods/shared_docstrings.jl +++ b/src/methods/shared_docstrings.jl @@ -57,3 +57,9 @@ const PROGRESS_KEYWORD = """ const VERBOSE_KEYWORD = """ - `vebose`: whether to print messages about potential problems. `true` by default. """ + +const GEOMCOLUMN_KEYWORD = """ +- `geomcolumn`: `Symbol` to manually select the column the geometries are in + when `data` is a Tables.jl compatible table, or a tuple of `Symbol` for columns + of point coordinates. +""" diff --git a/src/utils.jl b/src/utils.jl index 3280fab59..197e8c7b6 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,5 +1,5 @@ filter_ext(path, ext::AbstractString) = - filter(fn -> splitext(fn)[2] == ext, readdir(path; join=true)) + filte(fn -> splitext(fn)[2] == ext, readdir(path; join=true)) filter_ext(path, exts::Union{Tuple,AbstractArray}) = filter(fn -> splitext(fn)[2] in exts, readdir(path; join=true)) filter_ext(path, ext::Nothing) = readdir(path; join=true) diff --git a/test/runtests.jl b/test/runtests.jl index eaaf8d944..6e50b771a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -23,7 +23,7 @@ end @time @safetestset "reproject" begin include("reproject.jl") end @time @safetestset "warp" begin include("warp.jl") end @time @safetestset "resample" begin include("resample.jl") end - +@time @safetestset "kernel density" begin include("kerneldensity.jl") end # Only test SMAP locally for now, also RasterDataSources because CI dowloads keep breaking @time @safetestset "ncdatasets" begin include("sources/ncdatasets.jl") end From 1b16bb169be40676048942d2e535fd638f260e9c Mon Sep 17 00:00:00 2001 From: rafaqz Date: Sun, 24 Dec 2023 01:35:30 +0100 Subject: [PATCH 2/5] better tests --- src/utils.jl | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 197e8c7b6..e4cdfde37 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -128,15 +128,14 @@ function _without_mapped_crs(f, st::AbstractRasterStack, mappedcrs::GeoFormat) return x end -function _extent2dims(to; size=nothing, res=nothing, crs=nothing, kw...) - _extent2dims(to, size, res, crs) -end -function _extent2dims(to::Extents.Extent, size::Nothing, res::Nothing, crs) - isnothing(res) && throw(ArgumentError("Pass either `size` or `res` keywords or a `Tuple` of `Dimension`s for `to`.")) -end -function _extent2dims(to::Extents.Extent, size, res, crs) - isnothing(res) || _size_and_res_error() -end +_extent2dims(to; size=nothing, res=nothing, crs=nothing, kw...) = _extent2dims(to, size, res, crs) +_extent2dims(to::DimTuple, size::Nothing, res::Nothing, crs) = setcrs(to, crs) +_extent2dims(to::DimTuple, size::Nothing, res::Nothing, crs::Nothing) = to +_extent2dims(to, size::Nothing, res::Nothing, crs) = _extent2dims(dims(to), size, res, crs) +_extent2dims(to, size, res, crs) = _extent2dims(Extents.extent(to), size, res, crs) +_extent2dims(to::Extents.Extent, size, res, crs) = _size_and_res_error() +_extent2dims(to::Union{Nothing,Extents.Extent}, size::Nothing, res::Nothing, crs) = + throw(ArgumentError("Pass either `size` or `res` keywords or a `Tuple` of `Dimension`s for `to`.")) function _extent2dims(to::Extents.Extent{K}, size::Nothing, res::Real, crs) where K tuple_res = ntuple(_ -> res, length(K)) _extent2dims(to, size, tuple_res, crs) From 02fa2cd544ab0c9dae1376233264c2a755ce1713 Mon Sep 17 00:00:00 2001 From: rafaqz Date: Sun, 24 Dec 2023 20:16:28 +0100 Subject: [PATCH 3/5] add test --- test/kerneldensity.jl | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) create mode 100644 test/kerneldensity.jl diff --git a/test/kerneldensity.jl b/test/kerneldensity.jl new file mode 100644 index 000000000..6e589f723 --- /dev/null +++ b/test/kerneldensity.jl @@ -0,0 +1,36 @@ +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) From 1ebcb972a6536d509437571debf9bab02f5d288f Mon Sep 17 00:00:00 2001 From: rafaqz Date: Sun, 24 Dec 2023 20:19:06 +0100 Subject: [PATCH 4/5] typo --- src/extensions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/extensions.jl b/src/extensions.jl index f725d5529..bee5ce9b9 100644 --- a/src/extensions.jl +++ b/src/extensions.jl @@ -199,7 +199,7 @@ cellsize(args...; kw...) = throw_extension_error(cellsize, "ArchGDAL", :RastersA """ kerneldensity(geoms; kw...) -Kalculate the kernel density of the points in `geoms` using KernelDensity.jl. +Calculate the kernel density of the points in `geoms` using KernelDensity.jl. Returns a `Raster`. From 4036f8c3d11e12b289fad8db5034afe694c480a0 Mon Sep 17 00:00:00 2001 From: rafaqz Date: Sun, 24 Dec 2023 20:49:09 +0100 Subject: [PATCH 5/5] fix typepo --- src/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index e4cdfde37..7699c6cc2 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,5 +1,5 @@ filter_ext(path, ext::AbstractString) = - filte(fn -> splitext(fn)[2] == ext, readdir(path; join=true)) + filter(fn -> splitext(fn)[2] == ext, readdir(path; join=true)) filter_ext(path, exts::Union{Tuple,AbstractArray}) = filter(fn -> splitext(fn)[2] in exts, readdir(path; join=true)) filter_ext(path, ext::Nothing) = readdir(path; join=true)