Calculation and hatching of the uncertainty in projected delta changes.

# Atlas notebooks

***

> This notebook reproduces and extends parts of the figures and products of the AR6-WGI Atlas. It is part of a notebook collection available at https://github.com/IPCC-WG1/Atlas for reproducibility and reusability purposes. This work is licensed under a [Creative Commons Attribution 4.0 International License](http://creativecommons.org/licenses/by/4.0).
>
> ![Creative Commons License >](https://i.creativecommons.org/l/by/4.0/88x31.png)

Calculation and hatching of the uncertainty in projected delta changes.#

08/07/2021

M. Iturbide (Santander Meteorology Group. Institute of Physics of Cantabria, CSIC-UC, Santander, Spain).

This notebook is an example of calculation and visualization of the AR6-WGI uncertainty methods (simple and advanced) for projected delta changes. Please refer to the AR6 WGI Cross-Chapter Box Atlas 1 for more information.

Example data available at auxiliary-material is used. The base hatching R-functions are defined in datasets-interactive-atlas/hatching-functions/hatching-functions.R.

Two main sections in this notebook provide two ways of creating hatched maps:

(A) An easy one-step way to create maps including uncertainty information
(B) An explicit step-by-step way to illustrate the inner workings of including uncertainty information

Load libraries and functions#

We will use the following climate4R libraries:

loadeR to load data ([The R-based climate4R Open Framework for Reproducible Climate Data Access and Post-Processing - ScienceDirect])
visualizeR for data visualization ([The R-based climate4R Open Framework for Reproducible Climate Data Access and Post-Processing - ScienceDirect])

We will also need:

sp to work with spatial data ([Pebesma and Bivand])

library(loadeR)
library(visualizeR)
library(sp)

The base hatching R-functions are defined in datasets-interactive-atlas/hatching-functions/hatching-functions.R. Additionally, a wrapper function applying the previous is also available in datasets-interactive-atlas/hatching-functions/AR6-WGI-hatching.R. To load the functions in the working environment use the source R base function as follows.

source("../datasets-interactive-atlas/hatching-functions/hatching-functions.R")
source("../datasets-interactive-atlas/hatching-functions/AR6-WGI-hatching.R")

The scripts are available for download:

(A) Easy one-step way to create maps including uncertainty information#

This one-step way uses the AR6.WGI.hatching function. It requires one or two data objects; the delta change and, if the advanced method is chosen, also the historical reference.

The data used in this section is available under the ‘’’sh notebooks/auxiliary-material ‘’’ folder. There are two NetCDFs of CMIP6 JJA (June-July-August) data for South America, containing the SSP5-8.5 precipitation delta change for the long term (2081-2100 w.r.t baseline 1995-2014) and annual JJA precipitation data for the historical pre-industrial reference 1850-1900, respectively.

You can use the dataInventory function from loadeR to check basic data information, e.g.:

di <- dataInventory("auxiliary-material/CMIP6_SAM_historical_JJA_pr_1850-1900_annual.nc")
str(di)

Note that there are 33 CMIP6 models (“member” dimension).

These data can be loaded into variables using loadGridData.

ref <- loadGridData("auxiliary-material/CMIP6_SAM_historical_JJA_pr_1850-1900_annual.nc", var = "pr")
delta <- loadGridData("auxiliary-material/CMIP6_SAM_ssp585_JJA_pr_2081-2100_delta.nc", var = "pr")

[2023-05-07 17:01:47] Defining geo-location parameters

[2023-05-07 17:01:47] Defining time selection parameters

[2023-05-07 17:01:47] Retrieving data subset ...

[2023-05-07 17:01:49] Done

[2023-05-07 17:01:53] Defining geo-location parameters

[2023-05-07 17:01:53] Defining time selection parameters

[2023-05-07 17:01:53] Retrieving data subset ...

[2023-05-07 17:01:53] Done

We only need to apply the AR6.WGI.hatching function, which includes the method parameter to select the uncertainty method. For instance, in the next example the simple method is selected (method = simple), which is based on model agreement. Optionally, a list of additional parameters from the map.hatching function (package visualizeR) can be passed to argument map.hatching.args to control the appearance of the hatches, i.e. density, direction, color, … (type ?map.hatching for specific documentation).

Finally, additional arguments from the spatialPlot function can be also included to control graphical parameters, such as the color theme (e.g. color.theme = "BrBG"). A list of color themes to use from RColorBrewer package can be obtained by issuing RColorBrewer::display.brewer.all().

pl.simple.sam <- AR6.WGI.hatching(delta = delta, method = "simple", 
                                  map.hatching.args = list(density = 2, lwd = 1.2), 
                                  color.theme = "BrBG")

pl.simple.sam

../_images/3fab889c557d1f2ff93596a647c707ec4c03a371054031aecd9503c07d52bf75.png

To plot the relative delta, it needs to be explicitly provided. In this case, it is available in this repository as auxiliary material. Add more graphical parameters until the desired output is obtained. E.g. we can manually set the range of the color bar.

rel.delta <- loadGridData("auxiliary-material/CMIP6_SAM_ssp585_JJA_pr_2081-2100_relative_delta.nc", var = "pr")
pl.simple.sam <- AR6.WGI.hatching(delta = delta, relative.delta = rel.delta, method = "simple", 
                                  map.hatching.args = list(density = 2, lwd = 1.2), 
                                  color.theme = "BrBG", set.max = 50, set.min = -50)

[2023-05-07 17:01:58] Defining geo-location parameters

[2023-05-07 17:01:58] Defining time selection parameters

[2023-05-07 17:01:58] Retrieving data subset ...

[2023-05-07 17:01:58] Done

pl.simple.sam

../_images/9b9a51f4837125f5e27522103c2aafec0196cdca6cf90b03b4558f7077b823fa.png

We can change to the advanced method easily (method = advanced). The advanced method is based on the presence/absence of signal with respect to the interannual variability of a reference period (here pre-industrial, 1850-1900). This method has two components:

The areas with absence of signals are hatched
The areas with presence of signal but absence of model agreement are crossed.

pl.advanced.sam <- AR6.WGI.hatching(delta = delta, historical.ref = ref, relative.delta = rel.delta, method = "advanced", 
                                  map.hatching.args = list(density = 1, lwd = 0.7), 
                                  color.theme = "BrBG", set.max = 50, set.min = -50)

pl.advanced.sam

../_images/ffef216d5e91aef0f6087ae2661695760f5fe236a540a49d4729e1152992d01b.png

To export the map you can use the following code:

pdf("mymap.pdf")
pl.advanced.sam
dev.off()

(B) Explicit step-by-step illustration of the inner workings to include uncertainty information#

The data used in this section is global and it is available under the ´´´sh notebooks/auxiliary-material ´´´ folder. The data consists of two NetCDFs containing CMIP6 JJA historical and SSP5-8.5 precipitation climatologies for the 1986-2005 and 2081-2100 periods, respectively.

Use the dataInventory function from loadeR to check basic data information:

di <- dataInventory("auxiliary-material/CMIP6_historical_JJA_pr_1995-2014.nc")

[2023-05-07 17:02:02] Doing inventory ...

[2023-05-07 17:02:06] Retrieving info for 'pr' (0 vars remaining)

[2023-05-07 17:02:06] Done.

Load the data into variables using loadGridData.

hist <- loadGridData("auxiliary-material/CMIP6_historical_JJA_pr_1995-2014.nc", var = "pr")
scen <- loadGridData("auxiliary-material/CMIP6_ssp585_JJA_pr_2081-2100.nc", var = "pr")

[2023-05-07 17:02:10] Defining geo-location parameters

[2023-05-07 17:02:10] Defining time selection parameters

[2023-05-07 17:02:10] Retrieving data subset ...

[2023-05-07 17:02:10] Done

[2023-05-07 17:02:14] Defining geo-location parameters

[2023-05-07 17:02:14] Defining time selection parameters

[2023-05-07 17:02:14] Retrieving data subset ...

[2023-05-07 17:02:14] Done

Calculate delta change#

Future delta changes are the arithmetic difference between future and historical time slice climatologies. Relative deltas can also be calculated with respect to the multi-model mean.

delta <- gridArithmetics(climatology(scen), climatology(hist), operator = "-")
ensemble.mean <- function(grid) aggregateGrid(grid, aggr.mem = list(FUN = mean, na.rm = T))
delta.ens <- ensemble.mean(delta)
hist.ens <- ensemble.mean(climatology(hist))
# Relative delta
rel.delta <- gridArithmetics(delta.ens, 
                             hist.ens, 
                             100, 
                             operator = c("/", "*"))

[2023-05-07 17:02:14] - Computing climatology...

[2023-05-07 17:02:20] - Done.

[2023-05-07 17:02:20] - Computing climatology...

[2023-05-07 17:02:25] - Done.

[2023-05-07 17:02:26] - Aggregating members...

[2023-05-07 17:02:27] - Done.

[2023-05-07 17:02:27] - Computing climatology...

[2023-05-07 17:02:32] - Done.

[2023-05-07 17:02:32] - Aggregating members...

[2023-05-07 17:02:32] - Done.

Check the output with a first visualization of the delta change for each member of the ensemble.

spatialPlot(delta, backdrop.theme = "coastline", main = "delta", 
           set.max = 5, set.min = -5, color.theme = "BrBG")

../_images/0f0bfe574b1d61fcd389287ba4011fad020685f1c4f67f2aaeb0f413e51ba110.png

And this would be the relative delta change for the ensemble mean.

spatialPlot(rel.delta, backdrop.theme = "coastline", main = "relative delta", 
           set.max = 60, set.min = -60, color.theme = "BrBG")

../_images/b64d39d5947a43a732339b01c3675119edf885f2cd39acc80830a6310a37dd1f.png

The “simple” uncertainty method#

We compute now the simple uncertainty measure (agreement function) and check the result visually. Note that argument th allows to set the agreement threshold.

simple <- aggregateGrid(delta, aggr.mem = list(FUN = agreement, th = 80))
spatialPlot(simple)

This simple approach represents model agreement with two categories: 1 indicates high model agreement, where at least 80% of models agree on the sign of change; 0 indicates low model agreement, where fewer than 80% of models agree on sign of change. Note that model agreement is computed at a gridbox level and is not representative of regionally aggregated results over larger regions.

This binary mask can be used to create the spatial object of the hatches, which we will then overlay on the original plots, hatching the areas where there is not enough model agreement. Let’s obtain the hatching and check visually.

simple.hatch <- map.hatching(clim = climatology(simple), threshold = "0.5", angle = "-45",
                             condition = "LT", density = 4,  lwd = 0.6,
                             upscaling.aggr.fun = list(FUN = mean))
plot(simple.hatch[[2]])

[2023-05-07 17:02:41] - Computing climatology...

[2023-05-07 17:02:41] - Done.

../_images/f5d2d722bb6150e9152a9a5e67d7817122c53bac681d5d57e6521c597809000e.png

We can overlay the hatching on the relative delta to produce the final figure (the hatches are passed to argument sp.layout)

pl.simple <- spatialPlot(rel.delta, 
            color.theme = "BrBG", 
            at = seq(-60, 60, 5), 
            set.max = 60, set.min = -60,
            backdrop.theme = "coastline",
            main = list("Mean delta change", cex = 0.8),
            xlab = list("Period: 2081-2100, Season: JJA", cex = 0.8),
            sp.layout = list(simple.hatch),
            par.settings = list(axis.line = list(col = 'transparent')))
pl.simple

../_images/604bb322585c3dbf4dc606578c89013815e575387eb665fdd89b76c0eb0902a9.png

The “advanced” uncertainty method#

The advanced method is based on the presence/absence of signal with respect to the interannual variability of a reference period (here 1850-1900). This method has two components:

The areas with absence of signal are hatched
The areas with presence of signal but absence of model agreement are crossed.

To calculate the signal we need the data of the delta change and anual data of the historical reference. Due to file size restrictions on GitHub we do not include the global historical reference. However it would be loaded as follows:

ref <- loadGridData("auxiliary-material/CMIP6_historical_JJA_1850-1900_annual.nc", var = "pr")

Then we would calculate the signal as follows:

sign <- signal(h = ref, d = delta)

This repository includes the pre-calculated grid of the signal, load it by executing the following cell.

sign <- loadGridData("auxiliary-material/CMIP6_ssp585_JJA_pr_2081-2100_signal.nc", var = "pr_signal")

[2023-05-07 17:02:46] Defining geo-location parameters

[2023-05-07 17:02:46] Defining time selection parameters

[2023-05-07 17:02:46] Retrieving data subset ...

[2023-05-07 17:02:46] Done

To obtain the overall signal of the ensemble we use the auxiliary functions signal.ens1 and signal.ens2, which apply a threshold to the ensemble mean signal.

advanced1 <- aggregateGrid(sign, aggr.mem = list(FUN = signal.ens1, th = 66))

[2023-05-07 17:02:46] - Aggregating members...

[2023-05-07 17:02:46] - Done.

advanced2.aux <- aggregateGrid(sign, aggr.mem = list(FUN = signal.ens2, th = 66))

[2023-05-07 17:02:46] - Aggregating members...

[2023-05-07 17:02:47] - Done.

signal.ens1 returns the first component of the advanced method (where 0 = absence of signal), while signal.ens2 returns an intermediate grid (where 0 = presence of signal) to calculate the second component in an additional step.

spatialPlot(advanced1, main = "no signal = 0")
spatialPlot(advanced2.aux, main = "no signal = 1")

../_images/f97f016f5559c751e055811f14e97c70b4bb344aa155b22e8289be152bfd41c9.png

../_images/dc7b6361a33a5e4e262f1b64bf65b01e5d92905e71c8990e6417418859311470.png

The additional step consists on combining this map (advanced2.aux) with the map we obtained by the simple hatching approach (simple). As a result, we get a three-valued map.

advanced2.aux <- gridArithmetics(advanced2.aux, simple, operator = "+") 
spatialPlot(advanced2.aux)

../_images/87d188009fa2f78314cdbc93cb18dfb66ac6dd482495b9fdb0fb266bcea97431.png

The different values arise from:

no signal (1) + agreement (1) = 2
no signal (1) + disagreement (0) = 1
signal (0) + agreement (1) = 1
signal (0) + disagreement (0) = 0

The second component of the advanced method hatches the forth outcome, i.e. the 0 value, where there is signal but there is no model agreement.

advanced2 <- binaryGrid(advanced2.aux, condition = "GT", threshold = 0)
spatialPlot(advanced2)

../_images/2cb66fad91399e671254c226a891221610ed4d772d53efd39ab4e61b7a74de2e.png

We can now create the hatches for the two components (advanced1 and advanced2) of the advanced hatching method. For the second component, we opted for using crosses to obscure the map. To do this, we simply apply again the map.hatching function, but using different angles (angle = "45" and angle = "-45") each time.

my.hatching <- function(grid, angle) map.hatching(grid, threshold = 0.8, angle = angle,
                                       condition = "LT", density = 4,  lwd = 0.6,
                                       upscaling.aggr.fun = list(FUN = mean))
advanced1.hatch <- my.hatching(climatology(advanced1), "-45")
advanced2.hatch <- my.hatching(climatology(advanced2),  "45")
advanced2.hatch.bis <- my.hatching(climatology(advanced2), "-45")

[2023-05-07 17:02:48] - Computing climatology...

[2023-05-07 17:02:48] - Done.

[2023-05-07 17:02:48] - Computing climatology...

[2023-05-07 17:02:48] - Done.

[2023-05-07 17:02:49] - Computing climatology...

[2023-05-07 17:02:49] - Done.

plot(advanced1.hatch[[2]])
plot(advanced2.hatch[[2]], add = TRUE)
plot(advanced2.hatch.bis[[2]], add = TRUE)

../_images/ddf37be18163ede3b054bc5979b4cc41ee9d86f7e6f9776c43a9177043c85816.png

We can overlay the hatching on the relative delta to produce the final figure (the hatches are passed to argument sp.layout)

pl.advanced <- spatialPlot(rel.delta, 
            color.theme = "BrBG", 
            at = seq(-60, 60, 5), 
            set.max = 60, set.min = -60,
            backdrop.theme = "coastline",
            main = list("Mean delta change", cex = 0.8),
            xlab = list("Period: 2081-2100, Season: JJA", cex = 0.8),
            sp.layout = list(advanced1.hatch, advanced2.hatch, advanced2.hatch.bis),
            par.settings = list(axis.line = list(col = 'transparent')))
pl.advanced

../_images/740dca20e1a5a6ae2ac86ea94b6645e237b6f55b1341232f55d730ff0175dc7b.png

We could also add the hatching to the ensemble mean of the absolute deltas (instead of the relative ones) just replacing in the above code the first line by:

pl <- spatialPlot(ensemble.mean(delta), ...

Session Information#

sessionInfo()

Show code cell output Hide code cell output

R version 3.6.3 (2020-02-29)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Fedora Linux 38 (Workstation Edition)

Matrix products: default
BLAS/LAPACK: /home/phanaur/mambaforge/envs/tfg/lib/libopenblasp-r0.3.21.so

locale:
 [1] LC_CTYPE=es_ES.UTF-8          LC_NUMERIC=C                 
 [3] LC_TIME=es_ES.UTF-8           LC_COLLATE=es_ES.UTF-8       
 [5] LC_MONETARY=es_ES.UTF-8       LC_MESSAGES=es_ES.UTF-8      
 [7] LC_PAPER=es_ES.UTF-8          LC_NAME=es_ES.UTF-8          
 [9] LC_ADDRESS=es_ES.UTF-8        LC_TELEPHONE=es_ES.UTF-8     
[11] LC_MEASUREMENT=es_ES.UTF-8    LC_IDENTIFICATION=es_ES.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] sp_1.4-5            visualizeR_1.6.1    transformeR_2.1.3  
[4] loadeR_1.7.1        climate4R.UDG_0.2.3 loadeR.java_1.1.1  
[7] rJava_1.0-4        

loaded via a namespace (and not attached):
 [1] viridis_0.6.1           maps_3.3.0              jsonlite_1.7.2         
 [4] viridisLite_0.4.0       SpecsVerification_0.5-3 dotCall64_1.0-1        
 [7] kohonen_3.0.10          sm_2.2-5.6              latticeExtra_0.6-29    
[10] pillar_1.6.1            lattice_0.20-44         glue_1.4.2             
[13] RcppEigen_0.3.3.9.1     uuid_0.1-4              digest_0.6.27          
[16] RColorBrewer_1.1-2      colorspace_2.0-1        htmltools_0.5.1.1      
[19] Matrix_1.3-3            pkgconfig_2.0.3         padr_0.5.3             
[22] purrr_0.3.4             scales_1.1.1            jpeg_0.1-8.1           
[25] CircStats_0.2-6         dtw_1.22-3              tibble_3.1.2           
[28] proxy_0.4-25            generics_0.1.0          ggplot2_3.3.3          
[31] ellipsis_0.3.2          repr_1.1.3              verification_1.42      
[34] pbapply_1.4-3           magrittr_2.0.1          crayon_1.4.1           
[37] easyVerification_0.4.4  evaluate_0.14           fansi_0.4.2            
[40] MASS_7.3-54             tools_3.6.3             data.table_1.14.0      
[43] lifecycle_1.0.0         munsell_0.5.0           akima_0.6-2.1          
[46] compiler_3.6.3          mapplots_1.5.1          rlang_0.4.11           
[49] grid_3.6.3              RCurl_1.98-1.3          pbdZMQ_0.3-5           
[52] IRkernel_1.2            spam_2.6-0              bitops_1.0-7           
[55] tcltk_3.6.3             base64enc_0.1-3         boot_1.3-28            
[58] vioplot_0.3.6           gtable_0.3.0            abind_1.4-5            
[61] R6_2.5.0                gridExtra_2.3           zoo_1.8-9              
[64] rgdal_1.5-16            dplyr_1.0.6             utf8_1.2.1             
[67] parallel_3.6.3          IRdisplay_1.0           Rcpp_1.0.6             
[70] fields_12.3             vctrs_0.3.8             png_0.1-7              
[73] tidyselect_1.1.1       

References#

Rba(1,2): The R-based climate4R open framework for reproducible climate data access and post-processing - ScienceDirect. https://doi.org/10.1016/j.envsoft.2018.09.009.
PB: Pebesma and Bivand. Rnews 2005.