Projected climate change signals and uncertainty under global warming levels#
This notebook is a reproducibility example of the IPCC-WGI AR6 Interactive Atlas products. This work is licensed under a Creative Commons Attribution 4.0 International License.
M. Iturbide (Santander Meteorology Group. Instituto de Física de Cantabria, CSIC-UC, Santander, Spain).
This notebook is a reproducibility example for the IPCC-WGI AR6 Interactive Atlas products.
This notebook works with the data available in this Hub. In particular, the IPCC-WGI AR6 Interactive Atlas Dataset, originally published at DIGITAL.CSIC for the long-term archival, and also available through the Copernicus Data Store (CDS).
Open the Getting_started.ipynb for a description of the available data and for a quick **introduction to the climate4R framework, the bundle of R packages used in this notebook.
Contents in this notebook#
Libraries
The Global Warming Level analysis dimension
Data loading for the different GWLs
Data loading for the historical reference
Uncertainty calculation and representation
Before we start, or at any stage of the notebook, we can customize the plotting area within this notebook as follows:
library(repr)
# Change plot size
options(repr.plot.width=18, repr.plot.height=8)
1. Libraries#
The core climate4R libraries that allow for data loading, transformation (e.g. spatio-temporal aggregations), and visualization are loadeR, transformeR and visualizeR. In this example notebook we will also use the geoprocessoR package to perform spatial operations.
We will also use other useful libraries; the plotting library lattice, RColorBrewer for selecting color palettes, and the magrittr library for piping operations (through %>%).
library(loadeR)
library(transformeR)
library(visualizeR)
library(magrittr)
library(RColorBrewer)
library(lattice)
Loading required package: rJava
Loading required package: loadeR.java
Java version 20x amd64 by N/A detected
NetCDF Java Library v4.6.0-SNAPSHOT (23 Apr 2015) loaded and ready
Loading required package: climate4R.UDG
climate4R.UDG version 0.2.6 (2023-06-26) is loaded
Please use 'citation("climate4R.UDG")' to cite this package.
loadeR version 1.8.2 (2024-06-04) is loaded
Please use 'citation("loadeR")' to cite this package.
_______ ____ ___________________ __ ________
/ ___/ / / / |/ / __ /_ __/ __/ / / / / __ /
/ / / / / / /|_/ / /_/ / / / / __/ / /_/ / /_/_/
/ /__/ /__/ / / / / __ / / / / /__ /___ / / \ \
\___/____/_/_/ /_/_/ /_/ /_/ \___/ /_/\/ \_\
github.com/SantanderMetGroup/climate4R
transformeR version 2.2.2 (2023-10-26) is loaded
WARNING: Your current version of transformeR (v2.2.2) is not up-to-date
Get the latest stable version (2.2.3) using <devtools::install_github('SantanderMetGroup/transformeR')>
Please see 'citation("transformeR")' to cite this package.
visualizeR version 1.6.4 (2023-10-26) is loaded
Please see 'citation("visualizeR")' to cite this package.
2. The Global Warming Level analysis dimension#
Instead of calculating climate change anomalies over a fixed period, we will compute the anomaly based on a specific level of global warming (GWL). To achieve this, we need the information on the time windows during which global surface temperature reaches various levels of warming. These time windows vary across Global Climate Models (GCMs).
This data is available at the IPCC-WGI/Atlas GitHub repository. Specifically, we require the information in the CMIP6_Atlas_WarmingLevels.csv file, as we will be working with climate projections from CMIP6.
We can read the file remotely as follows:
url <- "https://github.com/SantanderMetGroup/ATLAS/raw/refs/heads/main/warming-levels/CMIP6_Atlas_WarmingLevels.csv"
gwls <- read.csv(url)
str(gwls)
'data.frame': 35 obs. of 17 variables:
$ model_run : chr "ACCESS-CM2_r1i1p1f1" "ACCESS-ESM1-5_r1i1p1f1" "AWI-CM-1-1-MR_r1i1p1f1" "BCC-CSM2-MR_r1i1p1f1" ...
$ X1.5_ssp126: int 2027 2030 2022 2041 NA 2013 2026 2021 2023 2027 ...
$ X2_ssp126 : int 2042 2072 2050 NA NA 2025 2042 2038 2037 2058 ...
$ X3_ssp126 : int NA NA NA NA NA NA NA NA NA NA ...
$ X4_ssp126 : int NA NA NA NA NA NA NA NA NA NA ...
$ X1.5_ssp245: int 2028 2029 2020 2035 2055 2012 2028 2024 2024 2030 ...
$ X2_ssp245 : int 2040 2045 2039 2057 2088 2023 2042 2039 2038 2048 ...
$ X3_ssp245 : int 2070 NA NA NA NA 2049 2078 2075 2066 2084 ...
$ X4_ssp245 : int NA NA NA NA NA 2083 NA NA NA NA ...
$ X1.5_ssp370: int 2027 2033 2021 2032 2046 2012 2030 2027 2025 2032 ...
$ X2_ssp370 : int 2039 2048 2037 2046 2065 2023 2043 2041 2039 2045 ...
$ X3_ssp370 : int 2062 2069 2064 2074 NA 2043 2068 2063 2063 2066 ...
$ X4_ssp370 : int 2082 NA NA NA NA 2059 2087 2085 2087 2083 ...
$ X1.5_ssp585: int 2025 2027 2019 2031 2047 2011 2023 2020 2021 2028 ...
$ X2_ssp585 : int 2038 2039 2036 2043 2064 2022 2034 2033 2033 2040 ...
$ X3_ssp585 : int 2055 2060 2059 2065 2088 2040 2055 2053 2052 2058 ...
$ X4_ssp585 : int 2071 2078 2079 NA NA 2054 2071 2068 2069 2072 ...
In this example, we will focus on the +3ºC GWL. We will use the ssp585 scenario, though any other scenario could be used, as the anomalies do not vary significantly across scenarios when analyzed through the GWL framework.
gwl3 <- gwls[,c("model_run", "X3_ssp585")] %>% print
model_run X3_ssp585
1 ACCESS-CM2_r1i1p1f1 2055
2 ACCESS-ESM1-5_r1i1p1f1 2060
3 AWI-CM-1-1-MR_r1i1p1f1 2059
4 BCC-CSM2-MR_r1i1p1f1 2065
5 CAMS-CSM1-0_r2i1p1f1 2088
6 CanESM5_r1i1p1f1 2040
7 CESM2_r4i1p1f1 2055
8 CESM2-WACCM_r1i1p1f1 2053
9 CMCC-CM2-SR5_r1i1p1f1 2052
10 CNRM-CM6-1_r1i1p1f2 2058
11 CNRM-CM6-1-HR_r1i1p1f2 2051
12 CNRM-ESM2-1_r1i1p1f2 2064
13 EC-Earth3_r1i1p1f1 2057
14 EC-Earth3-Veg_r1i1p1f1 2050
15 EC-Earth3-Veg-LR_r1i1p1f1 9999
16 FGOALS-g3_r1i1p1f1 2072
17 GFDL-CM4_r1i1p1f1 2059
18 GFDL-ESM4_r1i1p1f1 2075
19 HadGEM3-GC31-LL_r1i1p1f3 2047
20 IITM-ESM_r1i1p1f1 2075
21 INM-CM4-8_r1i1p1f1 2069
22 INM-CM5-0_r1i1p1f1 2074
23 IPSL-CM6A-LR_r1i1p1f1 2050
24 KACE-1-0-G_r2i1p1f1 2042
25 KIOST-ESM_r1i1p1f1 2063
26 MIROC-ES2L_r1i1p1f2 2070
27 MIROC6_r1i1p1f1 2076
28 MPI-ESM1-2-HR_r1i1p1f1 2073
29 MPI-ESM1-2-LR_r1i1p1f1 2071
30 MRI-ESM2-0_r1i1p1f1 2064
31 NESM3_r1i1p1f1 2054
32 NorESM2-LM_r1i1p1f1 2077
33 NorESM2-MM_r1i1p1f1 2076
34 TaiESM1_r1i1p1f1 2052
35 UKESM1-0-LL_r1i1p1f2 2046
The object gwl3 is a data frame that contains the central year within a 20-year time window during which the +3ºC GWL is reached for each model. For more information, refer to the IPCC-WGI/Atlas repository.
3. Data loading for the different GWLs#
To load the data, we first need to determine the path of the file containing the specific data we require. This path can point to either a NetCDF file or a NcML catalog file. The available paths are listed in the data_inventory.csv file. We can read this information with function read.csv and apply str for a quick overview of the content:
data.paths <- read.csv("../../../data_inventory.csv")
str(data.paths)
'data.frame': 1726 obs. of 6 variables:
$ location : chr "https://hub.climate4r.ifca.es/thredds/dodsC/ipcc/ar6/atlas/ia-monthly/CORDEX-ANT/historical/pr_CORDEX-ANT_histo"| __truncated__ "https://hub.climate4r.ifca.es/thredds/dodsC/ipcc/ar6/atlas/ia-monthly/CORDEX-ANT/historical/tn_CORDEX-ANT_histo"| __truncated__ "https://hub.climate4r.ifca.es/thredds/dodsC/ipcc/ar6/atlas/ia-monthly/CORDEX-ANT/historical/rx1day_CORDEX-ANT_h"| __truncated__ "https://hub.climate4r.ifca.es/thredds/dodsC/ipcc/ar6/atlas/ia-monthly/CORDEX-ANT/historical/tx_CORDEX-ANT_histo"| __truncated__ ...
$ type : chr "opendap" "opendap" "opendap" "opendap" ...
$ variable : chr "pr" "tn" "rx1day" "tx" ...
$ project : chr "CORDEX-ANT" "CORDEX-ANT" "CORDEX-ANT" "CORDEX-ANT" ...
$ experiment: chr "historical" "historical" "historical" "historical" ...
$ frequency : chr "mon" "mon" "mon" "mon" ...
locationrefers to the pathtyperefers to the access mode, local (netcdf) or remote (opendap).variablereferst to the climatic variable or index (i.e. tas, pr, tn, rx1day, etc.)projectrefers to the project coordinating the datasets (i.e. CORDEX-EUR, CMIP5, CMIP6, etc.)experimentreferst to the scenario (i.e. historical, rcp26, ssp126, rcp85, etc.)
We can easily apply filters to narrow down to the desired file:
filtered.path <- subset(data.paths, project == "CMIP6" & variable == "pr" & experiment == "ssp585")
filtered.path
| location | type | variable | project | experiment | frequency | |
|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | |
| 77 | https://hub.climate4r.ifca.es/thredds/dodsC/ipcc/ar6/atlas/ia-monthly/CMIP6/ssp585/pr_CMIP6_ssp585_mon_201501-210012.nc | opendap | pr | CMIP6 | ssp585 | mon |
| 940 | /home/jovyan/shared/data/IPCC-WGI_AR6_Interactive_Atlas_Dataset/CMIP6/ssp585/pr_CMIP6_ssp585_mon_201501-210012.nc | netcdf | pr | CMIP6 | ssp585 | mon |
As shown in the resulting table, we have two versions (type) of the dataset. The first version is accessed via OPeNDAP through THREDDS, while the second version is located in the shared directory of this Hub and can be accessed directly. In this example, we will choose the second option.
nc.file <- subset(filtered.path, type == "opendap")[["location"]]
Before data loading do an inventory of the data using the dataInventory function:
nc.file
di <- dataInventory(nc.file)
str(di)
[2025-01-09 12:34:05.166797] Doing inventory ...
[2025-01-09 12:34:05.205669] Opening dataset...
[2025-01-09 12:34:06.143759] The dataset was successfuly opened
[2025-01-09 12:34:06.235806] Retrieving info for 'pr' (0 vars remaining)
[2025-01-09 12:34:06.357078] Done.
List of 1
$ pr:List of 7
..$ Description: chr "Monthly mean of daily accumulated precipitation"
..$ DataType : chr "float"
..$ Shape : int [1:4] 33 1032 180 360
..$ Units : chr "mm"
..$ DataSizeMb : num 8827
..$ Version : logi NA
..$ Dimensions :List of 4
.. ..$ member:List of 3
.. .. ..$ Type : chr "Ensemble"
.. .. ..$ Units : chr ""
.. .. ..$ Values: chr [1:33] "CSIRO-ARCCSS_ACCESS-CM2_r1i1p1f1" "CSIRO_ACCESS-ESM1-5_r1i1p1f1" "BCC_BCC-CSM2-MR_r1i1p1f1" "CAMS_CAMS-CSM1-0_r2i1p1f1" ...
.. ..$ time :List of 4
.. .. ..$ Type : chr "Time"
.. .. ..$ TimeStep : chr "30.436 days"
.. .. ..$ Units : chr "days since 1850-01-01 00:00:00"
.. .. ..$ Date_range: chr "2015-01-01T00:00:00Z - 2100-12-01T00:00:00Z"
.. ..$ lat :List of 5
.. .. ..$ Type : chr "Lat"
.. .. ..$ Units : chr "degrees_north"
.. .. ..$ Values : num [1:180] -89.5 -88.5 -87.5 -86.5 -85.5 -84.5 -83.5 -82.5 -81.5 -80.5 ...
.. .. ..$ Shape : int 180
.. .. ..$ Coordinates: chr "lat"
.. ..$ lon :List of 5
.. .. ..$ Type : chr "Lon"
.. .. ..$ Units : chr "degrees_east"
.. .. ..$ Values : num [1:360] -180 -178 -178 -176 -176 ...
.. .. ..$ Shape : int 360
.. .. ..$ Coordinates: chr "lon"
The result includes information about the units, variables, dimensions, and more. Since we are working with warming levels and need to load a different period of years for each model, we will first extract the information about the members from the inventory. This information is essential to match the GWL array obtained earlier. In other words, we need to create an index object that provides the model position in the data for each row in the gwl3 object. We can achieve this using the grep function in a loop.
ind <- lapply(gwl3$model_run, grep, x = di$pr$Dimensions$member$Values)
This index (ind) will enable us to load the data for each model separately, requesting a different 20-year time period each time. Before loading the data, we can set common parameters such as the geographical domain, the target season (boreal winter in this example), and the reference historical baseline (in this case, the pre-industrial period) for simplicity.
ref.period <- 1850:1900
lonLim <- c(-11, 35)
latLim <- c(35, 72)
season <- c(12, 1, 2)
We will now apply the loadGridData function within a lapply loop to obtain a list of climate4R grids, where each slot represents a grid corresponding to a particular model. For a detailed explanation of what a climate4R grid is, please open the Getting_started.ipynb notebook.
The GWL period for each model is passed to the loadGridData function using the expression (gwl3$X3_ssp585[x]-9):(gwl3$X3_ssp585[x]+10), which generates a 20-year vector starting from the central year minus 9 years and ending at the central year plus 10 years (the central years are defined in the gwl3 object created earlier).
We also introduce a tryCatch expression to prevent the loop from being interrupted in cases where models do not reach the global warming level (GWL) within the current century. Additionally, we pipe the suppressMessages function to minimize lengthy output logs. This loading process can take a few minutes and varies depending on the region and period requested.
Note that if fixed periods are considered for the analysis, the process becomes much simpler because a loop is unnecessary, as the period to be loaded is common to all models. Please open the Getting_started.ipynb notebook for a simple example using fixed periods for the analysis.
cmip6.ssp585.l <- lapply(1:length(ind), function(x) {
loadGridData(dataset = nc.file,
var = "pr",
lonLim = lonLim,
latLim = latLim,
years = (gwl3$X3_ssp585[x]-9):(gwl3$X3_ssp585[x]+10),
season = season,
members = ind[[x]]) %>%
tryCatch(error = function(x) message("no data for this member and variable"))
}) %>% suppressMessages
We can eliminate the “no data” slots in the list as follows:
nodata.members <- lapply(cmip6.ssp585.l, is.null) %>% unlist %>% which
cmip6.ssp585.l <- cmip6.ssp585.l[-nodata.members]
In the next code chunk, we perform an annual aggregation of the variable. Additionally, we generate a time-series plot to illustrate the concept of GWLs, where each model spans a different time period. The temporalPlot function calculates the spatial mean by default
cmip6.ssp585.l.y <- lapply(cmip6.ssp585.l, aggregateGrid, aggr.y = list(FUN = "mean", na.rm = TRUE)) %>% suppressMessages
temporalPlot(cmip6.ssp585.l.y) %>% suppressMessages
To simplify the rest of the analysis, let’s combine all models into the same grid. We can achieve this by applying the bindGrid function and setting the skip.temporal.check parameter to TRUE, as each model covers a different time period.
cmip6.ssp585 <- bindGrid(cmip6.ssp585.l, dimension = "member", skip.temporal.check = TRUE)
Now we plot the map of the climatologies (mean of the period; 20-year means in this case) using the climatology and spatialPlot functions. The climatology function calculates the mean by default, though other functions can be specified to perform different aggregations over the period.
In the resulting plot, each panel represents a CMIP6 member (a model). Note that we previously set the plotting window space in Jupyter.
options(repr.plot.width=20, repr.plot.height=8)
spatialPlot(climatology(cmip6.ssp585),
backdrop.theme = "coastline",
at = seq(0, 12),
set.max = 12, set.min = 0,
layout = c(11, 3),
main = "Mean precipitation (mm/day) for +3ºC GWL",
strip = strip.custom(factor.levels = cmip6.ssp585$Members))
[2025-01-09 12:36:52.795996] - Computing climatology...
[2025-01-09 12:36:53.262786] - Done.
4. Data loading for the historical reference#
To get the climate change signal, we first need to load data from the historical scenario. To do this we repeat the filtering process to get the required NetCDF path.
nc.file.hist <- subset(data.paths, project == "CMIP6" & variable == "pr" & experiment == "historical" & type == "opendap")[["location"]] %>%
print
[1] "https://hub.climate4r.ifca.es/thredds/dodsC/ipcc/ar6/atlas/ia-monthly/CMIP6/historical/pr_CMIP6_historical_mon_185001-201412.nc"
As mentioned above, we are considering the pre-industrial period (defined earlier in object ref.period). In this case, the reference period is common to all models, therefore, the process becomes simpler; all models are loaded at once into a single grid object and we do not need to set a loop.
cmip6.hist <- loadGridData(nc.file.hist,
var = "pr",
lonLim = lonLim, latLim = latLim,
years = ref.period,
season = season) %>% suppressWarnings
[2025-01-09 12:36:59.843148] Opening dataset...
[2025-01-09 12:36:59.926658] The dataset was successfuly opened
[2025-01-09 12:36:59.931243] Defining geo-location parameters
[2025-01-09 12:36:59.996818] Defining time selection parameters
[2025-01-09 12:37:00.030582] Retrieving data subset ...
[2025-01-09 12:37:13.196329] Done
Before continuing, we must retain the same models (in the same order) in both the future and historical grids. To do so, we apply the intersectGrid function.
cmip6.common.mems <- intersectGrid(cmip6.hist, cmip6.ssp585, type = "members", which.return = 1:2)
cmip6.hist <- cmip6.common.mems[[1]] ; cmip6.ssp585 <- cmip6.common.mems[[2]]
It is always a good practice to make sure that the members in both grids are identical:
identical(cmip6.hist$Members, cmip6.ssp585$Members)
Now we can calculate the anomaly by computing the difference between both climatologies. We can easily do this using the climatology and gridArithmetics functions.
delta <- gridArithmetics(climatology(cmip6.ssp585), climatology(cmip6.hist), operator = "-")
[2025-01-09 12:37:13.564657] - Computing climatology...
[2025-01-09 12:37:14.315452] - Done.
[2025-01-09 12:37:14.40912] - Computing climatology...
[2025-01-09 12:37:15.077509] - Done.
Since the delta object results from an arithmetic operation between two climatologies it isn’t necessary to apply the climatology function to plot the corresponding maps.
spatialPlot(delta,
backdrop.theme = "coastline",
at = seq(-3, 3, 0.1),
set.max = 3, set.min = -3,
layout = c(11, 3),
main = "Mean precipitation anomaly for +3ºC GWL (rel. to preindustrial period)",
strip = strip.custom(factor.levels = delta$Members))
These maps represent the climate change signal as an absolute anomaly. We can calculate the relative anomaly as follows.
delta.rel <- gridArithmetics(delta, climatology(cmip6.hist), 100, operator = c("/", "*"))
[2025-01-09 12:37:21.048818] - Computing climatology...
[2025-01-09 12:37:22.086641] - Done.
To plot the resulting maps, we choose a specific Brewer palette, “BrBG” in this case (see the ColorBrewer website).
spatialPlot(delta.rel,
backdrop.theme = "coastline",
at = seq(-50, 50, 1),
set.max = 50, set.min = -50,
color.theme = "BrBG",
layout = c(11, 3),
main = "Mean precipitation relative anomaly for +3ºC GWL (rel. to preindustrial period)",
strip = strip.custom(factor.levels = delta$Members))
To calculate the multi-model mean, use function aggregateGrid.
ens.mean <- aggregateGrid(delta.rel, aggr.mem = list(FUN = mean))
[2025-01-09 12:37:27.821689] - Aggregating members...
[2025-01-09 12:37:27.832333] - Done.
# Change plot size
options(repr.plot.width=18, repr.plot.height=6)
spatialPlot(ens.mean,
backdrop.theme = "coastline",
at = seq(-50, 50, 1),
set.max = 50, set.min = -50,
color.theme = "BrBG",
main = "Mean precipitation relative anomaly for +3ºC GWL (rel. to preindustrial period) and the multi-model ensemble mean")
5. Uncertainty Calculation and Representation#
climate4R implements the “simple” and “advanced” methods for uncertainty characterization defined in the IPCC Sixth Assessment Report. The function to apply is computeUncertainty. Please refer to the AR6 WGI Cross-Chapter Box Atlas 1 (Gutiérrez et al., 2021) for more information.
In this example, we will use the ‘simple’ method.
uncert <- computeUncertainty(historical = aggregateGrid(cmip6.hist, aggr.y = list(FUN = "mean")),
anomaly = delta,
method = "simple")
spatialPlot(uncert)
[2025-01-09 12:37:28.164966] - Aggregating members...
[2025-01-09 12:37:28.182019] - Done.
The outcome yields a binary grid, readily convertible into a vectorial spatial object representing hatches. To achieve this transformation, utilize the map.hatching function. This function verifies a particular attribute within the grid, confirming its status as a “climatology” and ensuring that the length of the time dimension is 1. While our grid does not strictly meet the criteria of a climatology, its time dimension is singular. Therefore, we can successfully employ the function by applying the climatology function to our uncert object, as the values of the object remain the same.
cat1 <- map.hatching(clim = uncert %>% climatology,
condition = "LT",
threshold = 1,
density = 1) %>% suppressMessages
The newly created object (cat1) is a list containing a SpatialLines object from the sp R package. It also includes other plotting parameters, which can be fine-tuned within the map.hatching function. This object can be added as a layer to our map using the sp.layout parameter, placing it inside a list, as this parameter allows for the addition of multiple spatial layers.
spatialPlot(ens.mean,
backdrop.theme = "coastline",
sp.layout = list(cat1),
color.theme = "BrBG",
set.min = -50, set.max = 50,
at = seq(-50, 50, 1),
main = "CMIP6 - Total precipitation (PR) Change % - Warming 3°C SSP5-8.5 (rel. to 1850-1900)")
To export the figure run the following (change the output path as desired).
pdf("../../../CMIP6_Total-precipitation-relative-change_Warming-3degC_SSP585_rel-to-1850-1900.pdf", width = 10)
spatialPlot(ens.mean,
backdrop.theme = "coastline",
sp.layout = list(cat1),
color.theme = "BrBG",
set.min = -50, set.max = 50,
at = seq(-50, 50, 1),
main = "CMIP6 - Total precipitation (PR) Change % - Warming 3°C SSP5-8.5 (rel. to 1850-1900)")
dev.off()
Complete Workflow Code Summary#
# Library loading
library(loadeR)
library(transformeR)
library(visualizeR)
library(magrittr)
library(RColorBrewer)
# Read global warming level periods information
url <- "https://github.com/SantanderMetGroup/ATLAS/raw/refs/heads/main/warming-levels/CMIP6_Atlas_WarmingLevels.csv"
gwls <- read.csv(url)
gwl3 <- gwls[,c("model_run", "X3_ssp585")]
# Obtain the data path of future projections
data.paths <- read.csv("../../data_inventory.csv")
filtered.path <- subset(data.paths, project == "CMIP6" & variable == "pr" & experiment == "ssp585" & type == "opendap")
nc.file <- filtered.path[["location"]]
# Perform the inventory of the data within the NetCDF file
di <- dataInventory(nc.file)
# Match the list of members (CMIP models) in the global warming periods information and in the data file
ind <- lapply(gwl3$model_run, grep, x = di$pr$Dimensions$member$Values)
# LOAD FUTURE SSP585 PROJECTION DATA in a loop (each iteration loads a single member for its corresponding warming level period)
cmip6.ssp585.l <- lapply(1:length(ind), function(x) {
loadGridData(dataset = nc.file,
var = "pr",
lonLim = c(-11, 35),
latLim = c(35, 72),
years = (gwl3$X3_ssp585[x]-9):(gwl3$X3_ssp585[x]+10),
season = c(12, 1, 2),
members = ind[[x]]) %>%
tryCatch(error = function(x) message("no data for this member and variable"))
}) %>% suppressMessages
# Remove members with no data
nodata.members <- lapply(cmip6.ssp585.l, is.null) %>% unlist %>% which
cmip6.ssp585.l <- cmip6.ssp585.l[-nodata.members]
# Create a multi-member grid
cmip6.ssp585 <- bindGrid(cmip6.ssp585.l, dimension = "member", skip.temporal.check = TRUE)
# Obtain the data path for the historical experiment
nc.file.hist <- subset(data.paths, project == "CMIP6" & variable == "pr" & experiment == "historical" & type == "opendap")[["location"]]
# LOAD HISTORICAL DATA for the defined reference period
cmip6.hist <- loadGridData(nc.file.hist,
var = "pr",
lonLim = c(-11, 35),
latLim = c(35, 72),
years = 1850:1900,
season = c(12, 1, 2))
# Keep common members in the historical and future data
cmip6.common.mems <- intersectGrid(cmip6.hist, cmip6.ssp585, type = "members", which.return = 1:2)
cmip6.hist <- cmip6.common.mems[[1]] ; cmip6.ssp585 <- cmip6.common.mems[[2]]
# CALCULATE THE FUTURE RELATIVE ANOMALY
delta <- gridArithmetics(climatology(cmip6.ssp585), climatology(cmip6.hist), operator = "-")
delta.rel <- gridArithmetics(delta, climatology(cmip6.hist), 100, operator = c("/", "*"))
# Calculate the ensemble mean
ens.mean <- aggregateGrid(delta.rel, aggr.mem = list(FUN = mean))
# Uncertainty calculation
uncert <- computeUncertainty(historical = aggregateGrid(cmip6.hist, aggr.y = list(FUN = "mean")),
anomaly = delta,
method = "simple")
# Create the hatches
cat1 <- map.hatching(clim = uncert %>% climatology,
condition = "LT",
threshold = 1,
density = 1)
# FINAL MAP REPRESENTATION
spatialPlot(ens.mean,
backdrop.theme = "coastline",
sp.layout = list(cat1),
color.theme = "BrBG",
set.min = -50, set.max = 50,
at = seq(-50, 50, 1),
main = "CMIP6 - Total precipitation (PR) Change % - Warming 3°C SSP5-8.5 (rel. to 1850-1900)")
[2025-01-09 12:37:31.182602] Doing inventory ...
[2025-01-09 12:37:31.216494] Opening dataset...
[2025-01-09 12:37:31.290151] The dataset was successfuly opened
[2025-01-09 12:37:31.350713] Retrieving info for 'pr' (0 vars remaining)
[2025-01-09 12:37:31.421061] Done.
[2025-01-09 12:40:03.723499] Opening dataset...
[2025-01-09 12:40:03.794434] The dataset was successfuly opened
[2025-01-09 12:40:03.798872] Defining geo-location parameters
[2025-01-09 12:40:03.860172] Defining time selection parameters
Warning message in getTimeDomain(grid, dic, season, years, time, aggr.d, aggr.m, :
“First date in dataset: 1850-01-01. Seasonal data for the first year requested not available”
[2025-01-09 12:40:03.892769] Retrieving data subset ...
[2025-01-09 12:40:16.494596] Done
[2025-01-09 12:40:17.058449] - Computing climatology...
[2025-01-09 12:40:17.549951] - Done.
[2025-01-09 12:40:17.700109] - Computing climatology...
[2025-01-09 12:40:18.608355] - Done.
[2025-01-09 12:40:18.716086] - Computing climatology...
[2025-01-09 12:40:19.364917] - Done.
[2025-01-09 12:40:19.45488] - Aggregating members...
[2025-01-09 12:40:19.466504] - Done.
[2025-01-09 12:40:19.469941] - Aggregating members...
[2025-01-09 12:40:19.488055] - Done.
[2025-01-09 12:40:19.490852] - Computing climatology...
[2025-01-09 12:40:19.501907] - Done.
Session Info#
sessionInfo()
R version 4.3.3 (2024-02-29)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS
Matrix products: default
BLAS/LAPACK: /opt/conda/envs/climate4r/lib/libopenblasp-r0.3.27.so; LAPACK version 3.12.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=en_US.UTF-8
[9] LC_ADDRESS=en_US.UTF-8 LC_TELEPHONE=en_US.UTF-8
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=en_US.UTF-8
time zone: Etc/UTC
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lattice_0.22-6 RColorBrewer_1.1-3 magrittr_2.0.3
[4] visualizeR_1.6.4 transformeR_2.2.2 loadeR_1.8.1
[7] climate4R.UDG_0.2.6 loadeR.java_1.1.1 rJava_1.0-11
[10] repr_1.1.7
loaded via a namespace (and not attached):
[1] dotCall64_1.1-1 spam_2.10-0 latticeExtra_0.6-30
[4] vctrs_0.6.5 mapplots_1.5.2 tools_4.3.3
[7] bitops_1.0-8 generics_0.1.3 parallel_4.3.3
[10] tibble_3.2.1 proxy_0.4-27 fansi_1.0.6
[13] SpecsVerification_0.5-3 pkgconfig_2.0.3 KernSmooth_2.23-24
[16] data.table_1.15.4 uuid_1.2-1 lifecycle_1.0.4
[19] compiler_4.3.3 kohonen_3.0.12 deldir_2.0-4
[22] fields_16.2 munsell_0.5.1 akima_0.6-3.4
[25] class_7.3-22 htmltools_0.5.8.1 maps_3.4.2
[28] RCurl_1.98-1.16 pillar_1.9.0 crayon_1.5.3
[31] MASS_7.3-60.0.1 classInt_0.4-10 boot_1.3-31
[34] abind_1.4-5 vioplot_0.5.0 padr_0.6.2
[37] verification_1.42 tidyselect_1.2.1 digest_0.6.37
[40] sf_1.0-16 dplyr_1.1.4 fastmap_1.2.0
[43] grid_4.3.3 colorspace_2.1-1 cli_3.6.3
[46] base64enc_0.1-3 utf8_1.2.4 e1071_1.7-16
[49] sm_2.2-6.0 IRdisplay_1.1 withr_3.0.1
[52] scales_1.3.0 sp_2.1-4 IRkernel_1.3.2
[55] jpeg_0.1-10 interp_1.1-6 pbdZMQ_0.3-13
[58] RcppEigen_0.3.4.0.2 zoo_1.8-12 png_0.1-8
[61] pbapply_1.7-2 evaluate_1.0.0 dtw_1.23-1
[64] viridisLite_0.4.2 rlang_1.1.4 Rcpp_1.0.13
[67] DBI_1.2.3 glue_1.7.0 easyVerification_0.4.5
[70] CircStats_0.2-6 jsonlite_1.8.9 R6_2.5.1
[73] units_0.8-5