Skip to content

Commit

Permalink
unify
Browse files Browse the repository at this point in the history
  • Loading branch information
lazarusA committed Dec 28, 2024
1 parent db76f94 commit f460816
Show file tree
Hide file tree
Showing 3 changed files with 99 additions and 49 deletions.
129 changes: 92 additions & 37 deletions docs/src/tutorials/methods/cellarea.jl → docs/src/manual/cellarea.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
#=
# `cellarea` tutorial
## `cellarea`

```@meta
CollapsedDocStrings=true
Expand All @@ -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 <!-- 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.
=#

````@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!

Expand All @@ -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.
=#
masking and cropping the raster to the geometry, and then computing the statistic.
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
=#


2 changes: 1 addition & 1 deletion docs/src/tutorials/gbif_wflow.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down

0 comments on commit f460816

Please sign in to comment.