[WIP] More docs

Project.toml

@@ -102,7 +102,8 @@ RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "Proj", "RasterDataSources", "SafeTestsets", "Shapefile", "Statistics", "Test", "ZarrDatasets"]
+test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "Proj", "RasterDataSources", "SafeTestsets", "Shapefile", "Statistics", "Test", "Unitful", "ZarrDatasets"]

docs/Project.toml

@@ -10,11 +10,13 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365"
 GBIF2 = "dedd4f52-e074-43bf-924d-d6bce14ad628"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
+GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36"
 Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"

docs/make.jl

@@ -1,7 +1,7 @@
 using Documenter, Rasters, Plots, Logging, Statistics, Dates, 
-    RasterDataSources, ArchGDAL, NCDatasets, CoordinateTransformations
+    RasterDataSources, ArchGDAL, Proj, NCDatasets, CoordinateTransformations
 import Makie, CairoMakie
-using DocumenterVitepress
+using DocumenterVitepress, Documenter, Literate
 using Rasters.LookupArrays, Rasters.Dimensions
 
 # Don't output huge svgs for Makie plots
@@ -18,6 +18,20 @@ function flush_info_and_warnings()
 end
 flush_info_and_warnings()
 
+# Convert the tutorials from .jl to .md
+# Note that this converts every file in the `docs/src/tutorials/` directory,
+# so be careful when adding new tutorials.
+
+# One could add a filtering statement when looping over the files to filter out certain names.
+
+tutorials_dir = joinpath(@__DIR__, "src", "tutorials")
+for (root, dirs, files) in walkdir(tutorials_dir)
+    for file in files
+        if splitext(file)[2] == ".jl"
+            Literate.markdown(joinpath(root, file), root; flavor = Literate.DocumenterFlavor())
+        end
+    end
+end
 
 Logging.disable_logging(Logging.Warn)
 
@@ -44,8 +58,6 @@ makedocs(
         devurl = "dev";
     ),
     draft = false,
-    source = "src",
-    build = "build",
     warnonly=true,
 )
 

docs/src/array_operations.md

@@ -1,3 +1,5 @@
+# Array operations
+
 Most base methods work as for regular julia `Array`s, such as `reverse` and
 rotations like `rotl90`. Base, statistics and linear algebra methods like `mean`
 that take a `dims` argument can also use the dimension name. 

docs/src/gbif_wflow.md

@@ -1,6 +1,8 @@
-Load occurrences for the Mountain Pygmy Possum using GBIF.jl
+# GBIF workflow
 
-## Load GBIF
+In this example, we'll load occurrences for the Mountain Pygmy Possum species using [GBIF2.jl](https://github.com/rafaqz/GBIF2.jl), an interface to the [Global Biodiversity Information Facility](https://www.gbif.org/), and extract environmental variables using BioClim data from [RasterDataSources.jl](https://github.com/EcoJulia/RasterDataSources.jl).
+
+## Load GBIF species data
 
 ````@example gbif
 using Rasters, GBIF2
@@ -38,10 +40,15 @@ foreach(i -> scatter!(p, coords; subplot=i, kw...), 1:4)
 p
 ````
 
-Then extract predictor variables and write to CSV.
+Then extract predictor variables.
 ````@example gbif
-using CSV
 predictors = collect(extract(se_aus, coords))
+````
+
+These are recognized as a table format in Julia, so we can write them to file using CSV.jl:
+
+````@example gbif
+using CSV
 CSV.write("burramys_parvus_predictors.csv", predictors)
 ````
 

docs/src/tutorials/methods/cellarea.jl

@@ -0,0 +1,111 @@
+#=
+# `cellarea` tutorial
+
+```@meta
+CollapsedDocStrings=true
+```
+
+```@docs; canonical=false
+cellarea
+```
+
+`cellarea` computes the area of each cell in a raster.  
+This is useful for a number of reasons - if you have a variable like 
+population per cell, or elevation ([spatially extensive variables](https://r-spatial.org/book/05-Attributes.html#sec-extensiveintensive)),
+you'll want to account for the fact that different cells have different areas.
+
+`cellarea` returns a Raster with the same x and y dimensions as the input, 
+where each value in the raster encodes the area of the cell (in meters by default).
+
+You can specify whether you want to compute the area in the plane of your projection 
+(`Planar()`), on a sphere of some radius (`Spherical(; radius=...)`).
+<!-- or on an ellipsoid 
+(`Geodetic()`), using the first argument.-->
+
+Let's construct a raster and see what this looks like!  We'll keep it in memory.
+
+The spherical <!-- and geodetic --> method relies on the [Proj.jl](https://github.com/JuliaGeo/Proj.jl) package to perform coordinate transformation, so that has to be loaded explicitly.
+=#
+
+using Rasters, Proj
+
+# To construct a raster, we'll need to specify the x and y dimensions.  These are called "lookups" in Rasters.jl.
+using Rasters.Lookups
+# We can now construct the x and y lookups.  Here we'll use a start-at-one, step-by-five grid.
+# Note that we're specifying that the "sampling", i.e., what the coordinates actually mean, 
+# is `Intervals(Start())`, meaning that the coordinates are the starting point of each interval.
+#
+# This is in contrast to `Points()` sampling, where each index in the raster represents the value at a sampling point;
+# here, each index represents a grid cell, which is defined by the coordinate being at the start.
+x = X(1:5:30; sampling = Intervals(Start()), crs = EPSG(4326))
+y = Y(50:5:80; sampling = Intervals(Start()), crs = EPSG(4326))
+# I have chosen the y-range here specifically so we can show the difference between spherical and planar `cellarea`.
+ras = Raster(ones(x, y); crs = EPSG(4326))
+# We can just call `cellarea` on this raster, which returns cell areas in meters, on Earth, assuming it's a sphere:
+cellarea(ras)
+# and if we plot it, you can see the difference in cell area as we go from the equator to the poles:
+using Makie, CairoMakie# , GeoMakie
+heatmap(ras; axis = (; aspect = DataAspect()))
+# We can also try this using the planar method, which simply computes the area of the rectangle using `area = x_side_length * y_side_length`:
+cellarea(ras, Planar())
+# Note that this is of course wildly inaccurate for a geographic dataset - but if you're working in a projected coordinate system, like polar stereographic or Mercator, this can be very useful (and a _lot_ faster)!
+
+#=
+## Usage example
+Let's get the rainfall over Chile, and compute the average rainfall per meter squared across the country for the month of June.
+
+We'll get the precipitation data across the globe from [WorldClim](https://www.worldclim.org/data/index.html), via [RasterDataSources.jl](https://github.com/EcoJulia/RasterDataSources.jl), and use the `month` keyword argument to get the June data.
+
+Then, we can get the geometry of Chile from [NaturalEarth.jl](https://github.com/JuliaGeo/NaturalEarth.jl), and use `Rasters.mask` to get the data just for Chile.
+=#
+
+using RasterDataSources, NaturalEarth
+
+precip = Raster(WorldClim{Climate}, :prec; month = 6)
+#
+all_countries = naturalearth("admin_0_countries", 10)
+chile = all_countries.geometry[findfirst(==("Chile"), all_countries.NAME)]
+# Let's plot the precipitation on the world map, and highlight Chile:
+f, a, p = heatmap(precip; colorrange = Makie.zscale(replace_missing(precip, NaN)), axis = (; aspect = DataAspect()))
+p2 = poly!(a, chile; color = (:red, 0.5))
+f
+# You can see Chile highlighted in red, in the bottom left quadrant.
+#
+# First, let's make sure that we only have the data that we care about, and crop and mask the raster so it only has values in Chile.
+# We can crop by the geometry, which really just generates a view into the raster that is bounded by the geometry's bounding box.
+cropped_precip = crop(precip; to = chile)
+# Now, we mask the data such that any data outside the geometry is set to `missing`.
+masked_precip = mask(cropped_precip; with = chile)
+heatmap(masked_precip)
+# This is a lot of missing data, but that's mainly because the Chile geometry we have encompasses the Easter Islands as well, in the middle of the Pacific.
+
+# Now, let's compute the average precipitation per square meter across Chile.
+# First, we need to get the area of each cell in square meters.  We'll use the spherical method, since we're working with a geographic coordinate system.  This is the default.
+areas = cellarea(masked_precip)
+masked_areas = mask(areas; with = chile)
+heatmap(masked_areas; axis = (; title = "Cell area in square meters"))
+# Now we can compute the average precipitation per square meter.  First, we compute total precipitation per grid cell:
+precip_per_area = masked_precip .* masked_areas
+# We can sum this to get the total precipitation per square meter across Chile:
+total_precip = sum(skipmissing(precip_per_area))
+# We can also sum the areas to get the total area of Chile (in this raster, at least).
+total_area = sum(skipmissing(masked_areas))
+# And we can convert that to an average by dividing by the total area:
+avg_precip = total_precip / total_area
+# So on average, Chile gets about 100mm of rain per square meter in June.
+# Let's see what happens if we don't account for cell areas:
+bad_total_precip = sum(skipmissing(masked_precip))
+bad_avg_precip = bad_total_precip / length(collect(skipmissing(masked_precip)))
+# This is misestimated!  This is why it's important to account for cell areas when computing averages.
+
+#=
+!!! note
+    If you made it this far, congratulations!  
+
+    It's interesting to note that we've replicated the workflow of `zonal` here.  
+    `zonal` is a more general function that can be used to compute any function over geometries, 
+    and it has multithreading built in.  
+
+    But fundamentally, this is all that `zonal` is doing under the hood - 
+    masking and cropping the raster to the geometry, and then computing the statistic.
+=#

src/methods/reproject.jl

@@ -2,15 +2,16 @@
 """
     reproject(obj; crs)
 
-Reproject the lookups of `obj` to a different crs. 
+Reproject the lookups (axes) of `obj` to a different crs. 
 
 This is a lossless operation for the raster data, as only the 
 lookup values change. This is only possible when the axes of source
 and destination projections are aligned: the change is usually from
 a [`Regular`](@ref) and an [`Irregular`](@ref) lookup spans.
 
 For converting between projections that are rotated, 
-skewed or warped in any way, use [`resample`](@ref).
+skewed or warped in any way, *or* if you want to re-sample the 
+_data_, use [`resample`](@ref).
 
 Dimensions without an `AbstractProjected` lookup (such as a `Ti` dimension)
 are silently returned without modification.
@@ -41,8 +42,9 @@ end
 """
     reproject(source::GeoFormat, target::GeoFormat, dim::Dimension, val)
 
-`reproject` uses ArchGDAL.reproject, but implemented for a reprojecting
-a value array of values, a single dimension at a time.
+`reproject` uses Proj.jl's `Transformation` interface, 
+but implemented for reprojecting a lookup / axis array, 
+a single dimension at a time.
 """
 function reproject(source, target, dim, val)
     if source == target

test/cellarea.jl

@@ -1,6 +1,7 @@
 using Rasters, DimensionalData, Rasters.Lookups, Proj
 using Test
 using DimensionalData: @dim, YDim
+import Unitful
 include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
 
 @testset "cellarea" begin
@@ -70,4 +71,29 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
     # test missing crs throws an error
     nocrsdimz = setcrs(dimz, nothing)
     @test_throws ArgumentError cellarea(nocrsdimz)
+
+    @testset "Unitful cellarea" begin
+         # Case 1: cellarea with unitful manifold
+         # This is the simplest case
+         unitful_manifold = Spherical(; radius = Spherical().radius * Unitful.u"m")
+         unitful_cellarea = cellarea(unitful_manifold, dimz_3857)
+         @test unitful_cellarea isa Raster{<:Unitful.Quantity}
+         @test Unitful.ustrip.(parent(unitful_cellarea)) == cellarea(Spherical(; radius = unitful_manifold.radius |> Unitful.ustrip), dimz_3857)
+         # Case 2: unitful dimensions
+         ux_3857 = rebuild(x_3857; val = rebuild(val(x_3857); data = val(x_3857) .* Unitful.u"m", span = Regular(val(x_3857).span.step * Unitful.u"m")))
+         uy_3857 = rebuild(y_3857; val = rebuild(val(y_3857); data = val(y_3857) .* Unitful.u"m", span = Regular(val(y_3857).span.step * Unitful.u"m")))
+         unitful_dimz_3857 = (ux_3857, uy_3857)
+         @test cellarea(Planar(), unitful_dimz_3857) isa Raster{<:Unitful.Quantity}
+         @test Unitful.ustrip.(parent(cellarea(Planar(), unitful_dimz_3857))) == parent(cellarea(Planar(), dimz_3857))
+         # Unitful lookups shouldn't work without a unitful manifold
+         @test_throws Unitful.DimensionError cellarea(unitful_dimz_3857)
+         # The tests below fail because Unitful apparently can't convert to Float64...
+         # (see https://github.com/PainterQubits/Unitful.jl/issues/742)
+         # But we also have to re-convert back to the original unit type, which cellarea currently
+         # doesn't do.
+
+         # unitful_cellarea = cellarea(unitful_manifold, unitful_dimz_3857)
+         # @test unitful_cellarea isa Raster{<:Unitful.Quantity}
+         # @test Unitful.ustrip.(unitful_cellarea) == cellarea(unitful_manifold, dimz_3857)
+    end
 end