From 6d1b11256f58deb28a8f2450f4733a925d9c1a39 Mon Sep 17 00:00:00 2001 From: tiemvanderdeure Date: Wed, 15 Jan 2025 16:03:44 +0100 Subject: [PATCH 1/3] add multidimensional spatial means tutorial --- docs/src/tutorials/spatial_mean.md | 51 ++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/docs/src/tutorials/spatial_mean.md b/docs/src/tutorials/spatial_mean.md index 5c26d9a45..f914ce11a 100644 --- a/docs/src/tutorials/spatial_mean.md +++ b/docs/src/tutorials/spatial_mean.md @@ -141,3 +141,54 @@ We've also seen how to use the `cellarea` function to compute the area of each c We've seen that the spatial mean is not the same as the arithmetic mean, and that we need to account for the area of each cell when computing the average. +## Bonus: Computing spatial means across dimensions + +As a next step, we would like to know how precipitation will change in Chile until the end of the 21st century. To do this, we can use climate model outputs. This data can come from multiple climate models (GCMs) and under different socio-economic scenarios (SSPs). We'll use additional dimensions to keep track of these. + +First we define a simple function takes an SSP (socioeconomic scenario) and a GCM (climate model) as input, and downloads the appropriate climate data. + +````@example zonal +using Dates +getfutureprec(ssp, gcm) = Raster(WorldClim{Future{Climate, CMIP6, gcm, ssp}}, :prec, date = Date(2090)) +```` + +We will leverage some tools from [DimensionalData](https://github.com/rafaqz/DimensionalData.jl), which is the package that underlies Rasters.jl. Rather than having a seperate Raster for each combination of GCM and SSP, `gcm` and `ssp` will be additional dimensions, and our Raster will be 4-dimensional (X-Y-gcm-ssp). + +To do this, we first define two dimensions that correspond to the SSPs and GCMs we are interested in, then use the `@d` macro from [DimensionalData](https://github.com/rafaqz/DimensionalData.jl) to preserve these dimensions as we get the data, and then combine all Rasters into a single object using `Rasters.combine` + +````@example cellarea +SSPs = Dim{:ssp}([SSP126, SSP370]) # SSP126 is a low-emission scenario, SSP370 is a high-emission scenario +GCMs = Dim{:gcm}([GFDL_ESM4, IPSL_CM6A_LR]) # These are different general circulation (climate) models + +precip_future = (@d getfutureprec.(SSPs, GCMs)) |> RasterSeries |> Rasters.combine +```` + +Since the format of WorldClim's datasets for future climate is slightly different from the dataset for the historical period, this actually returned a 5-dimensional raster, with a `Band` dimension that represents months. Here we'll just select the 6th month, matching the selection above. We will also replace the `NaN` missing value by the more standard `missing` using [`replace_missing`](@ref). + +````@example cellarea +precip_future = precip_future[Band = 6] +precip_future = replace_missing(precip_future) +```` + +On our 4-dimensional raster, functions like `crop` and `mask`, as well as broadcasting, will still work. + +Here we repeat the procedure from above to mask out areas so we only have data for Chile, and then multiply by the cell area. + +````@example cellarea +masked_precip_future = mask(crop(precip_future; to = chile); with = chile) + +precip_litres_future = masked_precip_future .* areas +```` + +Now we calculate the average precipitation for each SSP and each GCM. Annoyingly, the future WorldClim doesn't have data for all land pixels, so we have to re-calculate the total area. + +````@example cellarea +masked_areas_future = mask(areas, with = masked_precip_future[ssp = 1, gcm = 1]) +total_area_f = sum(skipmissing(masked_areas_future)) + +avg_prec_future = map(eachslice(precip_litres_future; dims = (:ssp, :gcm))) do slice + sum(skipmissing(slice)) / total_area_f +end +```` + +Which shows us that June rainfall in Chile will be slightly lower in the future, especially under the high-emission SSP370 scenario. \ No newline at end of file From f2b3bddc9c31d623a2177324ecee4a9b9772404b Mon Sep 17 00:00:00 2001 From: tiemvanderdeure Date: Mon, 20 Jan 2025 12:17:52 +0100 Subject: [PATCH 2/3] a couple of text changes --- docs/src/tutorials/spatial_mean.md | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/src/tutorials/spatial_mean.md b/docs/src/tutorials/spatial_mean.md index f914ce11a..53582fce5 100644 --- a/docs/src/tutorials/spatial_mean.md +++ b/docs/src/tutorials/spatial_mean.md @@ -143,18 +143,18 @@ We've seen that the spatial mean is not the same as the arithmetic mean, and tha ## Bonus: Computing spatial means across dimensions -As a next step, we would like to know how precipitation will change in Chile until the end of the 21st century. To do this, we can use climate model outputs. This data can come from multiple climate models (GCMs) and under different socio-economic scenarios (SSPs). We'll use additional dimensions to keep track of these. +As a next step, we would like to know how precipitation will change in Chile until the end of the 21st century. To do this, we can use climate model outputs. This is a bit more complicated than calculating historical precipitation, because the forecast data can come from multiple climate models (GCMs), which each can be run under different socio-economic scenarios (SSPs). Here, we'll show how to use additional dimensions to keep track of this type of data. -First we define a simple function takes an SSP (socioeconomic scenario) and a GCM (climate model) as input, and downloads the appropriate climate data. +To start, we define a simple function that takes an SSP (socioeconomic scenario) and a GCM (climate model) as input, and return the appropriate climate data. ````@example zonal using Dates getfutureprec(ssp, gcm) = Raster(WorldClim{Future{Climate, CMIP6, gcm, ssp}}, :prec, date = Date(2090)) ```` -We will leverage some tools from [DimensionalData](https://github.com/rafaqz/DimensionalData.jl), which is the package that underlies Rasters.jl. Rather than having a seperate Raster for each combination of GCM and SSP, `gcm` and `ssp` will be additional dimensions, and our Raster will be 4-dimensional (X-Y-gcm-ssp). +Rather than having a seperate Raster object for each combination of GCM and SSP, we will do our analysis on a single Raster, which will have `gcm` and `ssp` as additional dimensions. In total, our Raster will have four dimensions: X, Y, gcm, and ssp. -To do this, we first define two dimensions that correspond to the SSPs and GCMs we are interested in, then use the `@d` macro from [DimensionalData](https://github.com/rafaqz/DimensionalData.jl) to preserve these dimensions as we get the data, and then combine all Rasters into a single object using `Rasters.combine` +To accomplish this, we will leverage some tools from [DimensionalData](https://github.com/rafaqz/DimensionalData.jl), which is the package that underlies Rasters.jl. We start by defining two dimensions that correspond to the SSPs and GCMs we are interested in, then use the `@d` macro from [DimensionalData](https://github.com/rafaqz/DimensionalData.jl) to preserve these dimensions as we get the data, and then combine all Rasters into a single object using `Rasters.combine`. ````@example cellarea SSPs = Dim{:ssp}([SSP126, SSP370]) # SSP126 is a low-emission scenario, SSP370 is a high-emission scenario @@ -163,7 +163,7 @@ GCMs = Dim{:gcm}([GFDL_ESM4, IPSL_CM6A_LR]) # These are different general circul precip_future = (@d getfutureprec.(SSPs, GCMs)) |> RasterSeries |> Rasters.combine ```` -Since the format of WorldClim's datasets for future climate is slightly different from the dataset for the historical period, this actually returned a 5-dimensional raster, with a `Band` dimension that represents months. Here we'll just select the 6th month, matching the selection above. We will also replace the `NaN` missing value by the more standard `missing` using [`replace_missing`](@ref). +Since the format of WorldClim's datasets for future climate is slightly different from the dataset for the historical period, this actually returned a 5-dimensional raster, with a `Band` dimension that represents months. Here we'll just select the 6th month, matching the selection above (but note that the analysis would also work for all Bands simultaneously). We will also replace the `NaN` missing value by the more standard `missing` using [`replace_missing`](@ref). ````@example cellarea precip_future = precip_future[Band = 6] From 3dbd8189c84b92a78236213c9b4ce353bc7481b3 Mon Sep 17 00:00:00 2001 From: tiemvanderdeure Date: Mon, 20 Jan 2025 12:20:37 +0100 Subject: [PATCH 3/3] text changes in first half --- docs/src/tutorials/spatial_mean.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/tutorials/spatial_mean.md b/docs/src/tutorials/spatial_mean.md index 53582fce5..a52fa3cbf 100644 --- a/docs/src/tutorials/spatial_mean.md +++ b/docs/src/tutorials/spatial_mean.md @@ -88,13 +88,13 @@ You can see here that cells are largest towards the equator, and smallest away f ## Computing the spatial mean -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 over each grid cell. (The units of this Raster will be m^2 * mm, which happens to be equal to liter.) ````@example cellarea 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 across Chile: ````@example cellarea total_precip = sum(skipmissing(precip_per_area)) @@ -106,7 +106,7 @@ We can also sum the areas to get the total area of Chile (in this raster, at lea 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 (in mm) by dividing by the total area: ````@example cellarea avg_precip = total_precip / total_area