From 25a771c36ca3b3c51e2f9e73a91c31ba1f4369df Mon Sep 17 00:00:00 2001 From: ipcc-wgi <89838559+ipcc-wgi@users.noreply.github.com> Date: Thu, 17 Nov 2022 11:54:07 +0100 Subject: [PATCH] Add files via upload --- Figure_TS24_column3.R | 102 +++++++++++++++++++++++++++++++++++++++ Figure_TS24_columns1-2.R | 88 +++++++++++++++++++++++++++++++++ 2 files changed, 190 insertions(+) create mode 100644 Figure_TS24_column3.R create mode 100644 Figure_TS24_columns1-2.R diff --git a/Figure_TS24_column3.R b/Figure_TS24_column3.R new file mode 100644 index 0000000..bd783b7 --- /dev/null +++ b/Figure_TS24_column3.R @@ -0,0 +1,102 @@ +# THIS SCRIPT REPRODUCES COLUMN 3 OF FIGURE TS.24 (IPCC-WGI AR6). +# THIS SCRIPT USES LIBRARIES FROM THE climate4R FRAMEWORK (Iturbide et al. 2019). Go to https://github.com/SantanderMetGroup/climate4R for installation and other information. +# THIS SCRIPT USES MATERIAL FROM THE IPCC-WGI/Atlas GitHub repository (DOI:10.5281/zenodo.5171760; clone or download: https://github.com/IPCC-WG1/Atlas). +# AUTHOR: Jorge Baño-Medina + +options(java.parameters = "-Xmx8g") + +# Load climate4R libraries ------------------------------------------------------------------ +library(loadeR) +library(transformeR) +library(visualizeR) +library(geoprocessoR) + +# Load other libraries ---------------------------------------------------------------------- +library(magrittr) +library(sp) +library(rgdal) +library(smoothr) +library(RColorBrewer) +library(gridExtra) + +# Function to compute model agreement ------------------------------------------------------- +consensus <- function(x, th = 80) { + mp <- mean(x, na.rm = TRUE) + if (is.na(mp)) { + 1 + } else { + if (mp > 0) { + as.numeric(sum(as.numeric(x > 0), na.rm = TRUE) > as.integer(length(x) * th / 100)) + } else if (mp < 0) { + as.numeric(sum(as.numeric(x < 0), na.rm = TRUE) > as.integer(length(x) * th / 100)) + } else if (mp == 0){ + 1 + } + } +} + +### Choose Data parameters ------------------------------------------------------------------------------ +scenario <- "rcp26" # either rcp26 or rcp85 +period <- "2041-2060" # either 2041-2060 or 2081-2100 +local_path = "final_data/" # path to netCDF data + +# set working directory +setwd(local_path) + +### Loading vectorial Spatial Objects for map plotting ------------------------------------------------------------------------------ +# Coast +coast <- readOGR("WORLD_coastline.shp") # available in https://github.com/IPCC-WG1/Atlas/tree/main/notebooks/auxiliary-material +proj4string(coast) <- CRS("+proj=longlat +datum=WGS84 +no_defs") +coast.rob <- spTransform(coast, CRSobj = CRS("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")) +# Contour +contour <- SpatialPolygons(list(Polygons(list(Polygon(matrix(c(-180, -180, 180, 180, -180, -90, 90, 90, -90, -90), ncol = 2))), ID = "A"))) +proj4string(contour) <- proj4string(coast) +contour <- densify(contour, max_distance = 0.44) +contour.rob <- spTransform(contour, CRS("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")) +# IPCC reference regions +regs <- get(load("./IPCC-WGI-reference-regions-v4_R.rda")) # available in https://github.com/IPCC-WG1/Atlas/blob/main/reference-regions/IPCC-WGI-reference-regions-v4_R.rda +regs <- as(regs, "SpatialPolygons") +regs.rob <- spTransform(regs, CRS("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")) + +### Sea-Land Mask (available in https://github.com/IPCC-WG1/Atlas/blob/main/reference-grids/land_sea_mask_05degree.nc4) ------------------------------------------------------------------------------ +mask <- loadGridData("land_sea_mask_05degree.nc4", var = "sftlf") %>% binaryGrid(condition = "LT",threshold = 1, values = c(1,NA)) + +### Loading tx35.nc files ------------------------------------------------------------------------------ +grid <- loadGridData(paste0("./CORDEX_TX35_", scenario, "_", period, "_rel.to-1995-2014.nc"), var = "tx35isimip") + +### Hatching ------------------------------------------------------------------------------ +mask2 <- aggregateGrid(grid, aggr.mem = list(FUN = "mean", na.rm = TRUE)) %>% gridArithmetics(mask) %>% binaryGrid(condition = "GT",threshold = -9999, values = c(NA,1)) # to include land-points with NA values, e.g., Arctic or Antarctica; in addition to sea points which are NA +uncer <- aggregateGrid(grid, aggr.mem = list(FUN = consensus, th = 80)) %>% gridArithmetics(mask2) +attr(uncer$xyCoords, "projection") <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" +uncer <- warpGrid(climatology(uncer), new.CRS = CRS("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")) +l <- c(map.hatching(clim = climatology(uncer), + threshold = 0.5, + condition = "GE", + density = 8, + angle = "45", coverage.percent = 100, + upscaling.aggr.fun = list(FUN = "mean", na.rm = TRUE) +), "which" = 1, lwd = 0.5) + +### Compute the ensemble mean and proyect grid to robin ------------------------------------------------------------------------------ +grid %<>% aggregateGrid(aggr.mem = list(FUN = "mean", na.rm = TRUE)) +grid$Data[which(is.na(grid$Data), arr.ind = TRUE)] <- NA +grid %<>% gridArithmetics(mask) +attr(grid$xyCoords, "projection") <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" +grid <- warpGrid(climatology(grid), new.CRS = CRS("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")) + +### We display the spatial maps ------------------------------------------------------------------------------ +spatialPlot(grid, + col.regions = c("#808080",rev(brewer.pal(n = 11, "RdBu"))[5:11]) %>% colorRampPalette(), + at = c(-32,-31,-30,-5,-1,1,5,15,30,45,60,75,100,150,200,250), + set.min = -31,set.max = 250, + main = paste("TX35 for",period,paste0("(",scenario,")"),"rel. to 1995-2014"), + par.settings = list(axis.line = list(col = 'transparent')), + sp.layout = list( + l, + list(regs.rob, first = FALSE, lwd = 0.6), + list(coast.rob, col = "gray50", first = FALSE, lwd = 0.6), + list(contour.rob, col = "black", first = FALSE, lwd = 0.7) + )) + + + diff --git a/Figure_TS24_columns1-2.R b/Figure_TS24_columns1-2.R new file mode 100644 index 0000000..e8acde0 --- /dev/null +++ b/Figure_TS24_columns1-2.R @@ -0,0 +1,88 @@ +# THIS SCRIPT REPRODUCES COLUMNS 1 AND 2 OF FIGURE TS.24 (IPCC-WGI AR6). +# THIS SCRIPT USES LIBRARIES FROM THE climate4R FRAMEWORK (Iturbide et al. 2019). Go to https://github.com/SantanderMetGroup/climate4R for installation and other information. +# THIS SCRIPT USES MATERIAL FROM THE IPCC-WGI/Atlas GitHub repository (DOI:10.5281/zenodo.5171760; clone or download: https://github.com/IPCC-WG1/Atlas). +# AUTHOR: Maialen Iturbide Martínez de Albéniz + +# Load climate4R libraries for data loading, changing spatial projection and creating map figures ---------------------------------------------------------------- +library(loadeR) +library(geoprocessoR) +library(visualizeR) + +# Load other libraries related to the management of spatial data and graphics visualization +# (use functions install.packages("name_of_the_package") or devtools::install_github("cran/name_of_the_package") to install them) ---------------------------------- +library(RColorBrewer) +library(rgdal) +library(smoothr) +library(gridExtra) + +# Directory containing the data (changed by the user) ------------------------------------------------------------------------------------------------------------ +data.dir <- "final_data/" + +# Auxiliary data for map drawing (clone or download: https://github.com/IPCC-WG1/Atlas) -------------------------------------------------------------------------- +regions <- get(load("MyDirectory/Atlas/reference-regions/IPCC-WGI-reference-regions-v4_R.rda")) +coast <- coast <- readOGR("MyDirectory/Atlas/notebooks/auxiliary-materialWORLD_coastline.shp") +world.contour <- SpatialPolygons(list(Polygons(list(Polygon(matrix(c(-180, -180, 180, 180, -180, -90, 90, 90, -90, -90), ncol = 2))), ID = "A"))) +proj4string(world.contour) <- proj4string(coast) +world.contour <- densify(world.contour, max_distance = 0.44) + +# Create file names objects +change.files <- c("CMIP5_TX35_rcp26_2041-2060_rel.to-1995-2014.nc", "CMIP5_TX35_rcp26_2081-2100_rel.to-1995-2014.nc", + "CMIP5_TX35_rcp85_2041-2060_rel.to-1995-2014.nc", "CMIP5_TX35_rcp85_2081-2100_rel.to-1995-2014.nc", + "CMIP6_TX35_ssp126_2041-2060_rel.to-1995-2014.nc", "CMIP6_TX35_ssp126_2081-2100_rel.to-1995-2014.nc", + "CMIP6_TX35_ssp585_2041-2060_rel.to-1995-2014.nc", "CMIP6_TX35_ssp585_2081-2100_rel.to-1995-2014.nc") +uncertainty.files <- c("CMIP5_TX35_rcp26_2041-2060_rel.to-1995-2014_model_agreement.nc", "CMIP5_TX35_rcp26_2081-2100_rel.to-1995-2014_model_agreement.nc", + "CMIP5_TX35_rcp85_2041-2060_rel.to-1995-2014_model_agreement.nc", "CMIP5_TX35_rcp85_2081-2100_rel.to-1995-2014_model_agreement.nc", + "CMIP6_TX35_ssp126_2041-2060_rel.to-1995-2014_model_agreement.nc", "CMIP6_TX35_ssp126_2081-2100_rel.to-1995-2014_model_agreement.nc", + "CMIP6_TX35_ssp585_2041-2060_rel.to-1995-2014_model_agreement.nc", "CMIP6_TX35_ssp585_2081-2100_rel.to-1995-2014_model_agreement.nc") + +# Load the data -------------------------------------------------------------------------------------------------------------------------------------------------- +change.grids <- lapply(paste0(data.dir, "/", change.files), loadGridData, var = "TX35") +uncertainty.grids <- lapply(paste0(data.dir, "/", uncertainty.files), loadGridData, var = "TX35_uncertainty") + +# Change projection to Robinson ---------------------------------------------------------------------------------------------------------------------------------- +change.grids.rob <- lapply(change.grids, function(x) warpGrid(x, new.CRS = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")) +uncertainty.grids.rob <- lapply(uncertainty.grids, function(x) warpGrid(x, new.CRS = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")) +regions.rob <- spTransform(regions, CRSobj = CRS("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")) +coast.rob <- spTransform(coast, CRSobj = CRS("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")) +world.contour.rob <- spTransform(world.contour, CRS("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")) + +hatching.density <- c(rep(3, 4), rep(6, 4)) +uncertainty.sp <- lapply(1:length(uncertainty.grids.rob), function(x) map.hatching(climatology(uncertainty.grids.rob[[x]]), threshold = 0.8, angle = "45", + condition = "LT", density = hatching.density[x], lwd = 0.8, col = "grey40", + upscaling.aggr.fun = list(FUN = mean))) + + +# Set graphic parameters ----------------------------------------------------------------------------------------------------------------------------------------- +legend.bins <- c(-30, -5, -1, 1, 5, 15, 30, 45, 60, 75, 100, 150, 200, 250) +coltheme = "RdBu" +revc <- TRUE +color.palette <- colorRampPalette(brewer.pal(coltheme, n = 10))(22) +titles <- c("(a) TX35 for 2041-2060 (RCP2.6) rel. to 1995-2014", "(b) TX35 for 2081-2100 (RCP2.6) rel. to 1995-2014", + "(c) TX35 for 2041-2060 (RCP8.5) rel. to 1995-2014", "(d) TX35 for 2081-2100 (RCP8.5) rel. to 1995-2014", + "(e) TX35 for 2041-2060 (SSP1-2.6) rel. to 1995-2014", "(f) TX35 for 2081-2100 (SSP1-2.6) rel. to 1995-2014", + "(g) TX35 for 2041-2060 (SSP5-8.5) rel. to 1995-2014", "(h) TX35 for 2081-2100 (SSP5-8.5) rel. to 1995-2014") + +# Plot maps ------------------------------------------------------------------------------------------------------------------------------------------------------ +p <- lapply(1:length(change.grids.rob), function(i) { + spatialPlot(change.grids.rob[[i]], backdrop.theme = "coastline", + col.regions = c(rev(color.palette[c(13, 15)]), "white", rev(color.palette[1:11])), + at = legend.bins, set.max = max(legend.bins), set.min = min(legend.bins), + main = list(titles[i], cex = 0.8), + colorkey = list(width = 1), + sp.layout = list(list(regions.rob, lwd = 0.9, first = F), + list(coast.rob, col = "gray50", lwd = 0.7, first = FALSE), + list(world.contour.rob, first = FALSE), + uncertainty.sp[[i]]), + par.settings = list(axis.line = list(col = 'transparent'))) +}) + + +# Export figure --------------------------------------------------------------------------------------------------------------------------------------------------- +pdf("FigureTS24-columns1-2.pdf", width = 8, height = 10) +do.call("grid.arrange", c(p[c(1,5,2,6,3,7,4,8)], ncol = 2, as.table = T)) +dev.off() + + + + +