Mosaic with polygons (#266)

* feat(sdt)!: zonal mosaic

* doc: some more example

* doc: show the top 10 districts

* doc: example with more polygons

* doc: map with polygons

* doc: move around

* doc: docstrings

* doc: fix boxes

* doc: update

* doc: update

* bug: GeoJSONT

* feat: work with Feature and FeaturesCollection

* feat: pass arguments to f in mosaic

* docstring

* feat: support for keyword args

* doc: tutorial text update

* doc: docstrings

* doc: fix output

docs/make.jl

@@ -44,6 +44,7 @@ makedocs(;
             "tutorials/pseudoabsences.md",
             "tutorials/fauxcurrences.md",
             "tutorials/futureclimate.md",
+            "tutorials/zonal.md",
         ],
         "How-to..." => [
             "howto/get_gbif_data.md",
@@ -54,9 +55,11 @@ makedocs(;
             "howto/makie.md",
         ],
         "Documentation" => [
-            "Work with species occurrence data" => "manual/SpeciesDistributionToolkit/index.md",
-            "Occurrences and layers" => "manual/SpeciesDistributionToolkit/gbif.jl.md",
-            "Pseudo-absences" => "manual/SpeciesDistributionToolkit/pseudoabsences.md",
+            "manual/index.md",
+            "manual/gbif.jl.md",
+            "manual/pseudoabsences.md",
+            "manual/polygons.md",
+            "manual/gadm.md",
         ]
     ],
 )

docs/src/manual.md

@@ -44,13 +44,3 @@ This will automatically install all the sub-packages.
 This manual is split into two sections: tutorials, which are medium to long
 examples of using the full functionality of the package; and how-tos, which are
 shorter (and denser) summaries of how to achieve a specific task.
-
-In addition, each of the component packages has its own documentation.
-
-## List of components
-
-- [GBIF.jl](https://poisotlab.github.io/SpeciesDistributionToolkit.jl/GBIF/)
-- [Phylopic.jl](https://poisotlab.github.io/SpeciesDistributionToolkit.jl/Phylopic/)
-- [SimpleSDMDatasets.jl](https://poisotlab.github.io/SpeciesDistributionToolkit.jl/SimpleSDMDatasets/)
-- [SimpleSDMLayers.jl](https://poisotlab.github.io/SpeciesDistributionToolkit.jl/SimpleSDMLayers/)
-- [Fauxcurrences.jl](https://poisotlab.github.io/SpeciesDistributionToolkit.jl/Fauxcurrences/)

docs/src/manual/gadm.md

@@ -0,0 +1,22 @@
+# Access GADM polygons
+
+The package come with a *very* lightweight series of convenience functions to
+interact with [GADM](https://gadm.org/). The [`GADM.jl`
+package](https://github.com/JuliaGeo/GADM.jl) is an alternative solution to the
+same problem.
+
+All methods assume that the first argument is an alpha-3 code valid under [ISO
+3166-1](https://www.iso.org/obp/ui/#search), and the following levels are
+sub-divisions of this territory.
+
+## Accessing polygons
+
+```@docs
+SpeciesDistributionToolkit.gadm
+```
+
+## Listing polygons
+
+```@docs
+SpeciesDistributionToolkit.gadmlist
+```

docs/src/manual/index.md

@@ -0,0 +1,62 @@
+# Documentation of integration functions
+
+These manual pages are focused on the functions integrating the different
+components of `SpeciesDistributionToolkit` together.
+
+The full docstrings for the different packages are in their own documentation
+webpages.
+
+::: details `GBIF.jl`
+
+A wrapper around the GBIF API, to retrieve taxa and occurrence datasets, and
+perform filtering on these occurrence data based on flags.
+
+[Read the documentation](https://poisotlab.github.io/SpeciesDistributionToolkit.jl/GBIF/)
+
+![GitHub Release](https://img.shields.io/github/v/release/poisotlab/speciesdistributiontoolkit.jl?filter=GBIF-*&style=flat-square&label=GBIF.jl) ![Lifecycle:Stable](https://img.shields.io/badge/Lifecycle-Stable-97ca00?style=flat-square)
+
+:::
+
+::: details `SimpleSDMDatasets.jl`
+
+An efficient way to download and store environmental raster data for consumption
+by other packages.
+
+[Read the documentation](https://poisotlab.github.io/SpeciesDistributionToolkit.jl/SimpleSDMDatasets/) 
+
+![GitHub Release](https://img.shields.io/github/v/release/poisotlab/speciesdistributiontoolkit.jl?filter=SimpleSDMDatasets-*&style=flat-square&label=SimpleSDMDatasets.jl) ![Lifecycle:Stable](https://img.shields.io/badge/Lifecycle-Stable-97ca00?style=flat-square)
+
+:::
+
+
+::: details `SimpleSDMLayers.jl`
+
+A series of types and common operations on raster data.
+
+[Read the documentation](https://poisotlab.github.io/SpeciesDistributionToolkit.jl/SimpleSDMLayers/)
+
+![GitHub Release](https://img.shields.io/github/v/release/poisotlab/speciesdistributiontoolkit.jl?filter=SimpleSDMLayers-*&style=flat-square&label=SimpleSDMLayers.jl) ![Lifecycle:Maturing](https://img.shields.io/badge/Lifecycle-Maturing-007EC6?style=flat-square)
+
+:::
+
+::: details `Fauxcurrences.jl`
+
+A package to simulate realistic species occurrence data from a know series of
+occurrences, with additional statistical constraints.
+
+[Read the documentation](https://poisotlab.github.io/SpeciesDistributionToolkit.jl/Fauxcurrences/)
+
+![GitHub Release](https://img.shields.io/github/v/release/poisotlab/speciesdistributiontoolkit.jl?filter=Fauxcurrences-*&style=flat-square&label=Fauxcurrences.jl) ![Lifecycle:Maturing](https://img.shields.io/badge/Lifecycle-Maturing-007EC6?style=flat-square)
+
+:::
+
+
+::: details `Phylopic.jl`
+
+A wrapper around the Phylopic API.
+
+[Read the documentation](https://poisotlab.github.io/SpeciesDistributionToolkit.jl/Phylopic/)
+
+![GitHub Release](https://img.shields.io/github/v/release/poisotlab/speciesdistributiontoolkit.jl?filter=Phylopic-*&style=flat-square&label=Phylopic.jl) ![Lifecycle:Stable](https://img.shields.io/badge/Lifecycle-Stable-97ca00?style=flat-square)
+
+:::

docs/src/manual/polygons.md

@@ -0,0 +1,27 @@
+# Working with polygons
+
+The package assumes that the polygons are read using [the `GeoJSON.jl`
+package](https://github.com/JuliaGeo/GeoJSON.jl). As per
+[RFC7946](https://datatracker.ietf.org/doc/html/rfc7946), the coordinates in the
+polygon are assumed to be WGS84.
+
+## Masking
+
+ ```@docs
+ mask!
+ mask
+ ```
+
+## Trimming
+
+```@docs
+trim
+```
+
+## Mosaic and zonal-like operations
+
+```@docs
+mosaic
+zone
+byzone
+```

docs/src/tutorials/polygons.jl

@@ -51,13 +51,13 @@ layer = SDMLayer(
     top = 55.0,
 )
 
-# We can check that this polygon is larger than we want:
+# We can check that this polygon is larger than the area we want:
 
 heatmap(layer; colormap = :navia, axis = (; aspect = DataAspect()))
 
 # We can now mask this layer according to the polygon:
 
-mask!(layer, CHE[1].geometry)
+mask!(layer, CHE)
 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:
@@ -100,7 +100,7 @@ current_figure() #hide
 # 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)
+scatter!(mask(presences, CHE); color = :orange, markersize = 4)
 current_figure() #hide
 
 # ::: details A note about vectors of occurrences

docs/src/tutorials/zonal.jl

@@ -0,0 +1,64 @@
+# # Zonal statistics
+
+# In this tutorial, we will
+
+using SpeciesDistributionToolkit
+using Statistics
+using CairoMakie
+CairoMakie.activate!(; type = "png", px_per_unit = 3.0) #hide
+
+#-
+
+spatial_extent = (left = 165.739746, bottom = -47.587547, right = 180.812988, top = -33.649514)
+
+#-
+
+dataprovider = RasterData(CHELSA2, BioClim)
+
+# precipitation of coldest quarter
+
+layer = SDMLayer(dataprovider; layer = "BIO19", spatial_extent...)
+
+#-
+
+mask!(layer, SpeciesDistributionToolkit.gadm("NZL"))
+layer = trim(layer)
+heatmap(layer; axis=(; aspect=DataAspect()))
+
+#-
+
+districts = SpeciesDistributionToolkit.gadm("NZL", 2);
+
+# We can start looking at how these map onto the landscape:
+
+heatmap(zone(layer, districts); colormap=:hokusai, axis=(; aspect=DataAspect()))
+
+# The layer resulting from the `zone` operation has integer values, and the
+# values correspond to the polygon to which each pixel belongs. Note that the
+# pixels that are not within a polygon are turned off.
+
+z = mosaic(median, layer, districts)
+nodata!(z, 0.0)
+
+#-
+
+heatmap(z; axis=(; aspect=DataAspect()))
+
+#-
+
+districtnames = SpeciesDistributionToolkit.gadmlist("NZL", 2)
+districtnames[1:10]
+
+#- 
+
+top5 = first.(sort(byzone(median, layer, districts, districtnames); by=(x) -> x.second, rev=true)[1:5])
+
+#-
+
+f, ax, plt = heatmap(layer; axis=(; aspect=DataAspect()), colormap=[:lightgrey, :black])
+[lines!(ax, districts[i], label=districtnames[i], linewidth=3) for i in indexin(top5, districtnames)]
+axislegend(position=(0, 0.7), nbanks=1)
+hidedecorations!(ax) #hide
+hidespines!(ax) #hide
+current_figure() #hide
+

src/SpeciesDistributionToolkit.jl

@@ -47,7 +47,10 @@ export pseudoabsencemask, backgroundpoints
 
 # Functions to deal with polygons
 include("polygons/polygons.jl")
+include("polygons/zonal.jl")
+include("polygons/mosaic.jl")
 include("polygons/gadm.jl")
 export trim
+export zone, byzone
 
 end # module SpeciesDistributionToolkit

src/polygons/gadm.jl

@@ -34,16 +34,32 @@ function _get_gadm_file(code, level)
     return GeoJSON.read(replace(gadm_file, ".zip" => ""))
 end
 
+"""
+    gadm(code::String)
+
+Returns the `GeoJSON` object associated to a the alpha-3 code defined by
+[ISO](https://www.iso.org/obp/ui/#search/code/).
+"""
 gadm(code::String) = _get_gadm_file(code, 0)
 
+"""
+    gadm(code::String, places::String...)
+
+Returns all polygons nested within an arbitrary sequence of areas according to
+GADM. For example, getting all counties in Oklahoma is achieved with the
+arguments `"USA", "Oklahoma"`.
+"""
 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)")))
+                findall(
+                    isequal(replace(places[i], " " => "")),
+                    getproperty(avail, Symbol("NAME_$(i)")),
+                )
                 for
                 i in 1:level
             ],
@@ -52,18 +68,95 @@ function gadm(code::String, places::String...)
     return avail.geometry[position]
 end
 
+"""
+    gadm(code::String, level::Integer)
+
+Returns all areas within the top-level territory `code`, at the level `level`.
+For example, getting all *départements* in France is done with the arguments
+`"FRA", 2`.
+"""
+function gadm(code::String, level::Integer)
+    avail = _get_gadm_file(code, level)
+    return avail.geometry
+end
+
+"""
+    gadm(code::String, level::Integer, places::String...)
+
+Returns all areas within the `places` within a country defined by `code` at a
+specific level. For example, the *districts* in the French region of Bretagne
+are obtained with the arguments `"FRA", 3, "Bretagne"`.
+"""
+function gadm(code::String, level::Integer, places::String...)
+    level = max(length(places), level)
+    avail = _get_gadm_file(code, level)
+    position = reduce(
+        intersect,
+        [
+            findall(
+                isequal(replace(p, " " => "")),
+                getproperty(avail, Symbol("NAME_$(i)")),
+            )
+            for
+            (i, p) in enumerate(places)
+        ],
+    )
+    return avail.geometry[position]
+end
+
+"""
+    gadmlist(code::String)
+
+Returns all top-level divisions of the territory defined by its `code`.
+"""
 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
+            findall(
+                isequal(replace(places[i], " " => "")),
+                getproperty(avail, Symbol("NAME_$(i)")),
+            ) for
             i in 1:(level - 1)
         ],
     )
     return getproperty(avail, Symbol("NAME_$(level)"))[position]
+end
+
+"""
+    gadmlist(code::String, level::Integer)
+
+Returns all `level` divisions of the territory defined by its `code`, regardless
+of which higher-level divisions they belong to.
+"""
+function gadmlist(code::String, level::Integer)
+    avail = _get_gadm_file(code, level)
+    return getproperty(avail, Symbol("NAME_$(level)"))
+end
+
+"""
+    gadmlist(code::String, level::Integer, places::String...)
+
+Returns all `level` divisions of the territory defined by its `code`, that
+belong to the hierarchy defined by `places`.
+"""
+function gadmlist(code::String, level::Integer, places::String...)
+    level = max(length(places), level)
+    avail = _get_gadm_file(code, level)
+    position = reduce(
+        intersect,
+        [
+            findall(
+                isequal(replace(p, " " => "")),
+                getproperty(avail, Symbol("NAME_$(i)")),
+            )
+            for
+            (i, p) in enumerate(places)
+        ],
+    )
+    return getproperty(avail, Symbol("NAME_$(level)"))[position]
 end