Add support to mask/reveal/hide by polygons

Project.toml

@@ -1,29 +1,36 @@
 name = "SpeciesDistributionToolkit"
 uuid = "72b53823-5c0b-4575-ad0e-8e97227ad13b"
 authors = ["Timothée Poisot <[email protected]>"]
-version = "0.1.1"
+version = "0.1.2"
 
 [deps]
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 Fauxcurrences = "a2d61402-033a-4ca9-aef4-652d70cf7c9c"
 GBIF = "ee291a33-5a6c-5552-a3c8-0f29a1181037"
+GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
 MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
 Phylopic = "c889285c-44aa-4473-b1e1-56f5d4e3ccf5"
+PolygonOps = "647866c9-e3ac-4575-94e7-e3d426903924"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SimpleSDMDatasets = "2c7d61d0-5c73-410d-85b2-d2e7fbbdcefa"
 SimpleSDMLayers = "2c645270-77db-11e9-22c3-0f302a89c64c"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe"
+ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
 
 [compat]
 Distances = "0.10"
 Fauxcurrences = "0.2"
 GBIF = "0.4"
+GeoJSON = "0.8"
 MakieCore = "0.8"
 Phylopic = "0.0"
+PolygonOps = "0.1"
 Reexport = "1.2"
 SimpleSDMDatasets = "0.2"
 SimpleSDMLayers = "1.0"
 StatsBase = "0.34"
 TestItems = "1"
+ZipFile = "0.10"
 julia = "1.8"

docs/Project.toml

@@ -5,6 +5,7 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
 Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"

docs/make.jl

@@ -40,6 +40,7 @@ makedocs(;
             "tutorials/consensus.md",
             "tutorials/layers_and_occurrences.md",
             "tutorials/bioclim.md",
+            "tutorials/polygons.md",
             "tutorials/pseudoabsences.md",
             "tutorials/fauxcurrences.md",
             "tutorials/futureclimate.md",

docs/src/tutorials/polygons.jl

@@ -0,0 +1,116 @@
+# # Working with polygons
+
+using SpeciesDistributionToolkit
+using CairoMakie
+import GeoJSON
+CairoMakie.activate!(; type = "png", px_per_unit = 3.0) #hide
+
+# In this tutorial, we will clip a layer to a polygon (in GeoJSON format), then
+# use the same polygon to filter GBIF records.
+
+# ::: warning About coordinates
+# 
+# [GeoJSON](https://geojson.org/) coordinates are expressed in WGS84. For this
+# reason, *any* polygon is assumed to be in this CRS, and all operations will be
+# done by projecting the layer coordinates to this CRS.
+# 
+# :::
+
+# We provide a very lightweight wrapper around the
+# [GADM](https://gadm.org/index.html) database, which will return data as
+# ready-to-use GeoJSON files. For example, we can get the borders of
+# Switzerland, as featured in [Damaris Zurell excellent SDM
+# tutorial](https://damariszurell.github.io/SDM-Intro/).
+
+CHE = SpeciesDistributionToolkit.gadm("CHE")
+
+# ::: details More about GADM
+# 
+# The only other function to interact with GADM is `gadmlist`, which takes either
+# a series of places, as in
+#
+SpeciesDistributionToolkit.gadmlist("FRA", "Bretagne")
+# 
+# or a level, to provide a list of places within a three-letter coded area at a
+# specific depth:
+# 
+SpeciesDistributionToolkit.gadmlist("FRA", 3)[1:3]
+# 
+# :::
+
+# The next step is to get a layer, and so we will download the data about
+# deciduous broadleaf trees from [EarthEnv](/datasets/EarthEnv#landcover):
+
+provider = RasterData(WorldClim2, Elevation)
+layer = SDMLayer(
+    provider;
+    resolution = 0.5,
+    left = 0.0,
+    right = 20.0,
+    bottom = 35.0,
+    top = 55.0,
+)
+
+# We can check that this polygon is larger than we want:
+
+heatmap(layer; colormap = :navia, axis = (; aspect = DataAspect()))
+
+# We can now mask this layer according to the polygon:
+
+mask!(layer, CHE[1].geometry)
+heatmap(layer; colormap = :navia, axis = (; aspect = DataAspect()))
+
+# This is a much larger layer than we need! For this reason, we will trim it so that the empty areas are removed:
+
+heatmap(trim(layer); colormap = :navia, axis = (; aspect = DataAspect()))
+
+# Let's now get some occurrences in the area defined by the layer boundingbox,
+# while specifying the dataset key for the [eBird Observation
+# Dataset](https://www.gbif.org/dataset/4fa7b334-ce0d-4e88-aaae-2e0c138d049e):
+
+ouzel = taxon("Turdus torquatus")
+presences = occurrences(
+    ouzel,
+    trim(layer),
+    "occurrenceStatus" => "PRESENT",
+    "limit" => 300,
+    "datasetKey" => "4fa7b334-ce0d-4e88-aaae-2e0c138d049e",
+)
+
+# ::: details Occurrences from a layer
+# 
+# The `GBIF.occurrences` method can accept a layer as its second argument, to
+# limit to the occurrence of a species within the bounding box of this layer. If
+# a layer is used as the sole argument, all occurrences in the bounding box are
+# queried.
+# 
+# :::
+
+while length(presences) < count(presences)
+    occurrences!(presences)
+end
+
+# We can plot the layer and all occurrences:
+
+heatmap(trim(layer); colormap = :navia, axis = (; aspect = DataAspect()))
+scatter!(presences; color = :orange, markersize = 4)
+current_figure() #hide
+
+# Some of these occurrences are outside of the masked region in the layer. For
+# this reason, we will use the *non-mutating* `mask` method on the GBIF records:
+
+heatmap(trim(layer); colormap = :navia, axis = (; aspect = DataAspect()))
+scatter!(mask(presences, CHE[1].geometry); color = :orange, markersize = 4)
+current_figure() #hide
+
+# ::: details A note about vectors of occurrences
+# 
+# The reason why `mask` called on a GBIF result is not mutating the result is
+# that GBIF results also store the query that was used. For this reason, it
+# makes little sense to modify this object. The non-mutating `mask` returns a
+# vector of GBIF records, which for most purposes can be used in-place of the
+# result.
+# 
+# :::
+
+#-

src/SpeciesDistributionToolkit.jl

@@ -18,6 +18,11 @@ using Reexport
 import Distances
 const _distance_function = Fauxcurrences._distancefunction
 
+import GeoJSON
+import PolygonOps
+import ZipFile
+import Downloads
+
 # Functions to get latitudes/longitudes
 include("latlon.jl")
 export latitudes, longitudes
@@ -40,4 +45,9 @@ include("pseudoabsences.jl")
 export WithinRadius, SurfaceRangeEnvelope, RandomSelection, DistanceToEvent
 export pseudoabsencemask, backgroundpoints
 
+# Functions to deal with polygons
+include("polygons/polygons.jl")
+include("polygons/gadm.jl")
+export trim
+
 end # module SpeciesDistributionToolkit

src/integrations/gbif_layers.jl

@@ -65,10 +65,46 @@ function SimpleSDMLayers.mask(layer::SDMLayer, records::GBIF.GBIFRecords)
     return out
 end
 
-function SimpleSDMLayers.mask(layer::SDMLayer, records::GBIF.GBIFRecords, ::Type{T}) where {T <: Number}
+function SimpleSDMLayers.mask(
+    layer::SDMLayer,
+    records::GBIF.GBIFRecords,
+    ::Type{T},
+) where {T <: Number}
     out = zeros(layer, T)
     for record in records
         out[record] += one(T)
     end
     return out
 end
+
+function _bbox_from_layer(layer::SDMLayer)
+    EL = eastings(layer)
+    NL = northings(layer)
+    prj = SimpleSDMLayers.Proj.Transformation(layer.crs, "+proj=longlat +datum=WGS84 +no_defs"; always_xy = true)
+
+    b1 = [prj(EL[1], n) for n in NL]
+    b2 = [prj(EL[end], n) for n in NL]
+    b3 = [prj(e, NL[1]) for e in EL]
+    b4 = [prj(e, NL[end]) for e in EL]
+    bands = vcat(b1, b2, b3, b4)
+
+    nx = extrema(first.(bands))
+    ny = extrema(last.(bands))
+    return (; left = nx[1], right = nx[2], bottom = ny[1], top = ny[2])
+end
+
+function GBIF.occurrences(t::GBIFTaxon, layer::SDMLayer, query::Pair...)
+    spatial_extent = SpeciesDistributionToolkit._bbox_from_layer(layer)
+    query = (query..., "decimalLatitude" => (spatial_extent.bottom, spatial_extent.top))
+    query = (query..., "decimalLongitude" => (spatial_extent.left, spatial_extent.right))
+    query = (query..., "hasCoordinate" => true)
+    return occurrences(t, query...)
+end
+
+function GBIF.occurrences(layer::SDMLayer, query::Pair...)
+    spatial_extent = SpeciesDistributionToolkit._bbox_from_layer(layer)
+    query = (query..., "decimalLatitude" => (spatial_extent.bottom, spatial_extent.top))
+    query = (query..., "decimalLongitude" => (spatial_extent.left, spatial_extent.right))
+    query = (query..., "hasCoordinate" => true)
+    return occurrences(query...)
+end

src/integrations/makie.jl

@@ -14,6 +14,12 @@ function sprinkle(records::GBIFRecords)
     return (lon, lat)
 end
 
+function sprinkle(records::Vector{GBIFRecord})
+    lon = Float32.(replace(longitudes.(records), missing => NaN))
+    lat = Float32.(replace(latitudes.(records), missing => NaN))
+    return (lon, lat)
+end
+
 function MakieCore.convert_arguments(::MakieCore.GridBased, layer::SDMLayer)
     return sprinkle(convert(SDMLayer{Float32}, layer))
 end
@@ -22,6 +28,8 @@ MakieCore.convert_arguments(P::MakieCore.NoConversion, layer::SDMLayer) =
     MakieCore.convert_arguments(P, values(layer))
 MakieCore.convert_arguments(P::MakieCore.PointBased, records::GBIFRecords) =
     MakieCore.convert_arguments(P, sprinkle(records)...)
+MakieCore.convert_arguments(P::MakieCore.PointBased, records::Vector{GBIFRecord}) =
+    MakieCore.convert_arguments(P, sprinkle(records)...)
 
 function MakieCore.convert_arguments(
     P::MakieCore.PointBased,

src/polygons/gadm.jl

@@ -0,0 +1,69 @@
+function _get_gadm_file(code, level)
+    # Prepare to store the data
+    gadm_root = joinpath(SimpleSDMDatasets._LAYER_PATH, "GADM", code)
+    if !ispath(gadm_root)
+        mkpath(gadm_root)
+    end
+    # Get the URL for the file
+    fname = "gadm41_$(code)_$(level).json"
+    if level > 0
+        fname *= ".zip"
+    end
+    url = "https://geodata.ucdavis.edu/gadm/gadm4.1/json/$(fname)"
+    # Download the file
+    if ~isfile(joinpath(gadm_root, fname))
+        try
+            Downloads.download(url, joinpath(gadm_root, fname))
+        catch
+            return nothing
+        end
+    end
+    # Read the file in GeoJSON
+    gadm_file = joinpath(gadm_root, fname)
+    if level > 0
+        gadm_zip = ZipFile.Reader(gadm_file)
+        for f in gadm_zip.files
+            if f.name == replace(fname, ".zip" => "")
+                out = open(joinpath(gadm_root, f.name), "w")
+                write(out, read(f, String))
+                close(out)
+            end
+        end
+    end
+    # Return
+    return GeoJSON.read(replace(gadm_file, ".zip" => ""))
+end
+
+gadm(code::String) = _get_gadm_file(code, 0)
+
+function gadm(code::String, places::String...)
+    level = length(places)
+    avail = _get_gadm_file(code, level)
+    position = only(
+        reduce(
+            intersect,
+            [
+                findall(isequal(replace(places[i], " " => "")), getproperty(avail, Symbol("NAME_$(i)")))
+                for
+                i in 1:level
+            ],
+        ),
+    )
+    return avail.geometry[position]
+end
+
+gadmlist(code::String) = getproperty(_get_gadm_file(code, 1), :NAME_1)
+gadmlist(code::String, level::Integer) =
+    getproperty(_get_gadm_file(code, level), Symbol("NAME_$(level)"))
+function gadmlist(code::String, places::String...)
+    level = length(places) + 1
+    avail = _get_gadm_file(code, level)
+    position = reduce(
+        intersect,
+        [
+            findall(isequal(replace(places[i], " " => "")), getproperty(avail, Symbol("NAME_$(i)"))) for
+            i in 1:(level - 1)
+        ],
+    )
+    return getproperty(avail, Symbol("NAME_$(level)"))[position]
+end