diff --git a/docs/src/tutorials/methods/cellarea.jl b/docs/src/manual/cellarea.md similarity index 51% rename from docs/src/tutorials/methods/cellarea.jl rename to docs/src/manual/cellarea.md index 2427afccd..47571ea63 100644 --- a/docs/src/tutorials/methods/cellarea.jl +++ b/docs/src/manual/cellarea.md @@ -1,5 +1,4 @@ -#= -# `cellarea` tutorial +## `cellarea` ```@meta CollapsedDocStrings=true @@ -25,80 +24,137 @@ You can specify whether you want to compute the area in the plane of your projec Let's construct a raster and see what this looks like! We'll keep it in memory. The spherical method relies on the [Proj.jl](https://github.com/JuliaGeo/Proj.jl) package to perform coordinate transformation, so that has to be loaded explicitly. -=# +````@example cella using Rasters, Proj +```` -# To construct a raster, we'll need to specify the x and y dimensions. These are called "lookups" in Rasters.jl. +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. +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. + +````@example cella 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`. +```` + +I have chosen the y-range here specifically so we can show the difference between spherical and planar `cellarea`. + +````@example cella 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: +```` + +We can just call `cellarea` on this raster, which returns cell areas in meters, on Earth, assuming it's a sphere: + +````@example cella cellarea(ras) -# and if we plot it, you can see the difference in cell area as we go from the equator to the poles: +```` + +and if we plot it, you can see the difference in cell area as we go from the equator to the poles: + +````@example cella 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`: +```` + +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`: + +````@example cella 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)! +```` + +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. -=# +````@example cella 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: +```` + +Let's plot the precipitation on the world map, and highlight Chile: + +````@example cella 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. +```` + +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. + +````@example cella cropped_precip = crop(precip; to = chile) -# Now, we mask the data such that any data outside the geometry is set to `missing`. +```` + +Now, we mask the data such that any data outside the geometry is set to `missing`. +````@example cella 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. +```` + +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. -# 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. +````@example cella 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: +```` + +Now we can compute the average precipitation per square meter. First, we compute total precipitation per grid cell: + +````@example cella precip_per_area = masked_precip .* masked_areas -# We can sum this to get the total precipitation per square meter across Chile: +```` + +We can sum this to get the total precipitation per square meter across Chile: + +````@example cella 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). +```` + +We can also sum the areas to get the total area of Chile (in this raster, at least). + +````@example cella total_area = sum(skipmissing(masked_areas)) -# And we can convert that to an average by dividing by the total area: +```` + +And we can convert that to an average by dividing by the total area: + +````@example cella 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: +```` + +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: + +````@example cella 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. +```` + +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! @@ -107,5 +163,4 @@ bad_avg_precip = bad_total_precip / length(collect(skipmissing(masked_precip))) 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. -=# \ No newline at end of file + masking and cropping the raster to the geometry, and then computing the statistic. \ No newline at end of file diff --git a/docs/src/tutorials/methods/resample_warp.jl b/docs/src/manual/resample_warp.md similarity index 98% rename from docs/src/tutorials/methods/resample_warp.jl rename to docs/src/manual/resample_warp.md index f424894ab..375c45623 100644 --- a/docs/src/tutorials/methods/resample_warp.jl +++ b/docs/src/manual/resample_warp.md @@ -1,15 +1,16 @@ -#= -# `resample` tutorial - warping a raster +## `resample` warping a raster ```@meta CollapsedDocStrings=true ``` -=# + +````@example resample using Rasters, ArchGDAL using RasterDataSources using NaNStatistics using CairoMakie -#= +```` + ## What is resampling? **[`resample`](@ref)** "re-samples" the @@ -50,12 +51,6 @@ and is generally pretty robust. However, it has the following limitations: resample ``` - - - - - - Topics: - What is resampling? - When to resample vs reproject @@ -68,5 +63,5 @@ Topics: - Show some different projections and ways of constructing them - Show how to use `size` and `res` to change the resolution of a raster - Show how to use `warp` to reproject a raster -=# + diff --git a/docs/src/tutorials/gbif_wflow.md b/docs/src/tutorials/gbif_wflow.md index 5bd6fcedc..1fe352007 100644 --- a/docs/src/tutorials/gbif_wflow.md +++ b/docs/src/tutorials/gbif_wflow.md @@ -2,7 +2,7 @@ This example shows a full Species distribution modelling workflow, from loading data, to cleaning it, to fitting an ensemble and generating predictions. -It uses GBIF and WorldClim data, which are common datasets in ecology. +It uses GBIF and WorldClim data, which are common datasets in ecology. 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 Rasters, ArchGDAL, RasterDataSources and GBIF The GBIF2 library is used to download occurrence data, RasterDataSources to conveniently access Bioclim data. ArchGDAL is necessary to load in the Bioclim data.