From c50a22fcb22f19202018eed7c4c22ebdc854ae56 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= Date: Wed, 18 Sep 2024 15:58:31 -0400 Subject: [PATCH] WIP --- docs/src/vignettes/simplesdmlayers.jl | 39 +++++++------------ src/extensions/SpeciesDistributionToolkit.jl | 40 +++++++------------- test/units/R1_simplesdmlayers.jl | 19 +++------- 3 files changed, 33 insertions(+), 65 deletions(-) diff --git a/docs/src/vignettes/simplesdmlayers.jl b/docs/src/vignettes/simplesdmlayers.jl index 62f307e..5dc41f3 100644 --- a/docs/src/vignettes/simplesdmlayers.jl +++ b/docs/src/vignettes/simplesdmlayers.jl @@ -22,7 +22,6 @@ using SpatialBoundaries using SpeciesDistributionToolkit using CairoMakie -import Plots # Required for radial histograms # First we can start by defining the extent of the Southwestern islands of # Hawaii, which can be used to restrict the extraction of the various landcover @@ -32,7 +31,7 @@ import Plots # Required for radial histograms hawaii = (left = -160.2, right = -154.5, bottom = 18.6, top = 22.5) dataprovider = RasterData(EarthEnv, LandCover) landcover_classes = SimpleSDMDatasets.layers(dataprovider) -landcover = [SimpleSDMPredictor(dataprovider; layer=class, full=true, hawaii...) for class in landcover_classes] +landcover = [SDMLayer(dataprovider; layer=class, full=true, hawaii...) for class in landcover_classes] # We can remove all the areas that contain 100% water from the landcover data as # our question of interest is restricted to the terrestrial realm. We do this by @@ -41,7 +40,8 @@ landcover = [SimpleSDMPredictor(dataprovider; layer=class, full=true, hawaii...) ow_index = findfirst(isequal("Open Water"), landcover_classes) not_water = landcover[ow_index] .!== 0x64 -lc = [mask(not_water, layer) for layer in landcover] +nodata!(not_water, false) +[mask!(layer, not_water) for layer in landcover] # As layers one through four of the EarthEnv data are concerned with data on # woody cover (*i.e.* "Evergreen/Deciduous Needleleaf Trees", "Evergreen @@ -51,7 +51,7 @@ lc = [mask(not_water, layer) for layer in landcover] # woody cover: classes_with_trees = findall(contains.(landcover_classes, "Trees")) -tree_lc = convert(Float32, reduce(+, lc[classes_with_trees])) +tree_lc = convert(SDMLayer{Float64}, reduce(+, landcover[classes_with_trees])) heatmap(tree_lc; colormap=:linear_kbgyw_5_98_c62_n256) # Although we have previously summed the four landcover layers for the actual @@ -61,7 +61,7 @@ heatmap(tree_lc; colormap=:linear_kbgyw_5_98_c62_n256) # woody cover layers. This will give as a vector containing four `LatticeWomble` # objects (since the input data was in the form of a matrix). -wombled_layers = wombling.(lc[classes_with_trees]); +wombled_layers = wombling.(landcover[classes_with_trees]); # As we are interested in overall woody cover for Southwestern islands we can # take the `wombled_layers` vector and use them with the `mean` function to get @@ -69,29 +69,21 @@ wombled_layers = wombling.(lc[classes_with_trees]); # cover. This will 'flatten' the four wombled layers into a single # `LatticeWomble` object. -wombled_mean = mean(wombled_layers); - -# From the `wombled_mean` object we can 'extract' the layers for both the mean -# rate and direction of change. For ease of plotting we will also convert these -# layers to `SimpleSDMPredictor` type objects. It is also possible to call these -# matrices directly from the `wombled_mean` object using `wombled_mean.m` or -# `wombled_mean.θ` for the rate and direction respectively. - -rate, direction = SimpleSDMPredictor(wombled_mean) +rate = mosaic(mean, (x -> x.rate).(wombled_layers)) +direction = mosaic(a -> atan(mean(sin.(deg2rad.(a))), mean(cos.(deg2rad.(a)))), (x -> x.direction).(wombled_layers)) # Lastly we can identify candidate boundaries using the `boundaries`. Here we # will use a thresholding value (t) of 0.1 and save these candidate boundary # cells as `b`. Note that we are now working with a `SimpleSDMResponse` object # and this is simply for ease of plotting. -b = similar(rate) -b.grid[boundaries(wombled_mean, 0.1; ignorezero = true)] .= 1.0 +(quantize(rate).>=0.1) |> heatmap # We will overlay the identified boundaries (in green) over the rate of change # (in levels of grey): heatmap(rate, colormap=[:grey95, :grey5]) -heatmap!(b, colormap=[:transparent, :green]) +#heatmap!(b, colormap=[:transparent, :green]) current_figure() # For this example we will plot the direction of change as radial plots to get @@ -117,14 +109,11 @@ direction_candidate = mask(b, direction) # collect the cells and convert them from degrees to radians. Then we can start # by plotting the direction of change of *all* cells. -Plots.stephist( - deg2rad.(values(direction_all)); - proj=:polar, - lab="", - c=:teal, - nbins = 36, - yshowaxis=false, - normalize = false) +f = Figure() +ax = PolarAxis(f[1,1]) +h = fit(Histogram, values(direction); nbins=50) +stairs!(ax, h.edges[1], h.weights[[axes(h.weights,1)...,1]]) +current_figure() # Followed by plotting the direction of change only for cells that are # considered as candidate boundary cells. diff --git a/src/extensions/SpeciesDistributionToolkit.jl b/src/extensions/SpeciesDistributionToolkit.jl index 12c867f..6df0ab6 100644 --- a/src/extensions/SpeciesDistributionToolkit.jl +++ b/src/extensions/SpeciesDistributionToolkit.jl @@ -1,37 +1,23 @@ """ - wombling(layer::T; convert_to::Type=Float64) where {T <: SimpleSDMLayer} + wombling(layer::SDMLayer; convert_to::Type=Float64) Performs a lattice wombling on a `SimpleSDMLayer`. """ -function wombling(layer::T; convert_to::Type = Float64) where {T <: SimpleSDMLayer} - try - global nan = convert(convert_to, NaN) - catch - throw(ArgumentError("The type given as `convert_to` must have a `NaN` value.")) - end - +function wombling(layer::SDMLayer) + # Get the values for x and y - y = collect(SimpleSDMLayers.longitudes(layer)) - x = collect(SimpleSDMLayers.latitudes(layer)) + y = collect(SimpleSDMLayers.eastings(layer)) + x = collect(SimpleSDMLayers.northings(layer)) # Get the grid - z = convert(Matrix{Union{Nothing, convert_to}}, layer.grid) - replace!(z, nothing => nan) - return wombling(x, y, convert(Matrix{convert_to}, z)) -end + z = convert(Matrix{Float64}, layer.grid) + z[findall(!, layer.indices)] .= NaN -# We want to make sure that the layers are returned without NaN values, and -# adding this method makes the code easier to write -Base.isnan(::Nothing) = false - -function SimpleSDMLayers.SimpleSDMPredictor(W::T) where {T <: LatticeWomble} - rate = SimpleSDMLayers.SimpleSDMPredictor(W.m, extrema(W.y)..., extrema(W.x)...) - direction = SimpleSDMLayers.SimpleSDMPredictor(W.θ, extrema(W.y)..., extrema(W.x)...) - rate.grid[findall(isnan, rate.grid)] .= nothing - direction.grid[findall(isnan, direction.grid)] .= nothing - return (rate, direction) -end + # Womble + W = wombling(Float64.(x), Float64.(y), z) -function SimpleSDMLayers.SimpleSDMResponse(W::T) where {T <: LatticeWomble} - return convert.(SimpleSDMLayers.SimpleSDMResponse, SimpleSDMLayers.SimpleSDMPredictor(W)) + # Prepare to return + rate = SimpleSDMLayers.SDMLayer(W.m, (!isnan).(W.m), extrema(W.y), extrema(W.x), layer.crs) + direction = SimpleSDMLayers.SDMLayer(W.θ, (!isnan).(W.m), extrema(W.y), extrema(W.x), layer.crs) + return (; rate, direction) end diff --git a/test/units/R1_simplesdmlayers.jl b/test/units/R1_simplesdmlayers.jl index 4debfe9..92aa665 100644 --- a/test/units/R1_simplesdmlayers.jl +++ b/test/units/R1_simplesdmlayers.jl @@ -4,9 +4,10 @@ using SpatialBoundaries using SpeciesDistributionToolkit using Test -precipitation = SimpleSDMPredictor( - RasterData(WorldClim2, BioClim), - layer=12; +precipitation = SDMLayer( + RasterData(WorldClim2, BioClim); + resolution=0.5, + layer="BIO12", left = -80.0, right = -56.0, bottom = 44.0, @@ -15,15 +16,7 @@ precipitation = SimpleSDMPredictor( W = wombling(precipitation) -@test isa(W, LatticeWomble) - -Lt, Ld = SimpleSDMPredictor(W) -@test isa(Lt, SimpleSDMPredictor) - -Rt, Rd = SimpleSDMResponse(W) -@test isa(Rt, SimpleSDMResponse) - -@test !any(map(isnan, Rt.grid)) -@test !any(map(isnan, Rd.grid)) +@test isa(W.rate, SDMLayer) +@test isa(W.direction, SDMLayer) end