deploy: 279675f

.buildinfo

@@ -0,0 +1,4 @@
+# Sphinx build info version 1
+# This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done.
+config: 9f285b88e02085586166ffcb4a9be7e1
+tags: 645f666f9bcd5a90fca523b33c5a78b7

_downloads/2724a26e8b24c2e2bf3dfc50a70eaa7a/gwl_time_series_plots.R

@@ -0,0 +1,172 @@
+# gwl_time_series_plots.R
+#
+# Copyright (C) 2021 Santander Meteorology Group (http://meteo.unican.es)
+#
+# This work is licensed under a Creative Commons Attribution 4.0 International
+# License (CC BY 4.0 - http://creativecommons.org/licenses/by/4.0)
+
+#' @title Global Warming Level time timings displayed as multi-model time series plots
+#' @description For a given CMIP project, experiment and filtering window width, produces
+#' a plot of smoothed delta time series (1850-2100, w.r.t. the PI period 1850-1900)
+#' for each ensemble member and draws its corresponding (central-year) warming level 
+#' 
+#' @author J. Bedia
+
+
+library(magrittr)
+library(RColorBrewer)
+
+#
+# Assume Atlas repo home as base directory
+#
+
+source("warming-levels/scripts/getGWL.R")
+
+#
+# Parameter settings (current values reproduce Fig. in the README file)
+#
+
+cmip <- "CMIP6"
+gwl <- 2 # Possible values c(1.5, 2 ,3, 4)
+exp <- "ssp370" # This is a CMIP-dependent parameter. See a few lines below.
+window <- 20 # window width for centered moving average.
+             # Default to 20 (as used in the IPCC AR6 Atlas products)
+
+cmip <- match.arg(cmip, choices = c("CMIP5", "CMIP6"))
+
+#
+# CMIP-dependent parameters
+#
+
+exp <- if (cmip == "CMIP6") {
+  match.arg(exp, choices =  c("ssp126", "ssp245", "ssp370", "ssp585"))
+} else {
+  match.arg(exp, choices = c("rcp26", "rcp45", "rcp85"))
+}
+
+
+#
+# Load filenames to process
+#
+
+datadir <- sprintf("./datasets-aggregated-regionally/data/%s/%s_tas_landsea", cmip, cmip)
+filelist <- list.files(datadir)
+allfiles <- sprintf("%s/%s", datadir, filelist)
+aux <- grep("historical", filelist, value = TRUE)
+modelruns <- gsub(sprintf("%s_|_historical|\\.csv", cmip), "", aux)
+
+# Remove run IDs to simplify the plot legend later:
+gcm.names <- gsub("_r.*", "", modelruns) 
+
+#
+# Main loop
+#
+
+# First, create an empty list populated with the ensemble members.
+# Models without the target rcp/ssp are discarded
+
+l <- list()
+counter <- 0L
+
+for (i in 1:length(modelruns)) {
+
+  modelfiles <- gsub("_", "_.*", modelruns[i], fixed = TRUE) %>%
+    grep(allfiles, value = TRUE)
+  hist <- grep("historical", modelfiles, value = TRUE) %>% world.annual.mean()
+  rcp <- grep(exp, modelfiles, value = TRUE)
+
+  if (length(rcp > 0L)) {
+
+    counter <- counter + 1L
+    message("[", Sys.time(), "] Processing ", modelruns[i])
+    rcp %<>% world.annual.mean() 
+    tas <- append(hist, rcp)
+
+    # Fill with NAs if needed to get a continuous annual series 1850-2100
+    aux <- rep(NA, length(1850:2100))
+    names(aux) <- 1850:2100
+    aux[which(names(aux) %in% names(tas))] <- tas
+
+    # Delta change w.r.t. the PI period (1850-1900)
+    baseline <- hist[1:51] %>% mean(na.rm = TRUE)
+    aux <- aux - baseline
+
+    # Compute GWL central year and store as attribute
+    attr(aux, "GWL") <- getGWL(aux, window = window, GWL = gwl)[1]
+    l[[gcm.names[i]]] <- aux
+
+  }
+}
+
+
+#
+# The final dataset is stored in tabular form as a data.frame
+#
+
+mat <- do.call("cbind.data.frame", l)
+
+#
+# Apply moving average to each ensemble member
+#
+
+filtered.mat <- apply(mat, MARGIN = 2L, FUN = "filter",
+                      filter = rep(1/window, window), sides = 2)
+
+#
+# Plotting
+#
+
+# Y-axis ranges 
+
+ylims <- range(filtered.mat, na.rm = TRUE)
+ylim <- c(ylims[1] - .5, ylims[2] + .5)
+
+
+# Choosing random distinct colors for each ensemble member
+# NOTE: color palette may not be fully reproducible due to the random color selection in this step
+# https://stackoverflow.com/questions/15282580/how-to-generate-a-number-of-most-distinctive-colors-in-r
+
+n <- ncol(filtered.mat)
+qual_col_pals <- brewer.pal.info[brewer.pal.info$category == 'qual',]
+col_vector <- unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))
+set.seed(2)
+model.colors <- sample(col_vector, n)
+
+
+# Empty plot
+
+plot(1850:2100, filtered.mat[,1], ty = 'n',
+     ylim = ylim, las = 1, ylab = "Delta change w.r.t. PI period 1850-1900 (degC)", xlab = "year")
+grid()
+abline(h = gwl, col = "grey60")
+
+# Add ensemble members 
+
+for (i in 1:ncol(filtered.mat)) {
+  lines(1850:2100, filtered.mat[ ,i], col = model.colors[i])
+  abline(v = attr(l[[i]], "GWL"), col = model.colors[i], lty = 2)
+}
+
+# Add legend
+
+legend(x = "topleft", 
+       legend = names(mat),
+       lty = 1, col = model.colors,
+       cex = .7, ncol = 2, bg = "grey90")  
+
+# Add title
+
+title(paste(cmip, exp, "- GWL +", gwl, "degC"))
+mtext(paste(window, "- year window width"), line = .25)
+
+# Add range as a text label
+
+rng <- sapply(l, "attr", "GWL") %>% range(na.rm = TRUE)
+txt <- if (any(is.infinite(rng))) {
+  paste0("+", gwl , " degC GWL not reached\n by any ensemble member") 
+} else {
+  paste("GWL range:", paste(rng, collapse = "-"))
+}
+text(1975, -0.5, txt, col = "blue")
+
+

_downloads/47ba7ef97df83c510bd24c4283b0c84c/AR6-WGI-hatching.R

@@ -0,0 +1,61 @@
+# AR6.WGI_hatching.R
+#
+# Copyright (C) 2021 Santander Meteorology Group (http://meteo.unican.es)
+#
+# This work is licensed under a Creative Commons Attribution 4.0 International
+# License (CC BY 4.0 - http://creativecommons.org/licenses/by/4.0)
+
+#' @title Hatching wrapper function applying atomic hatching functions
+#' @description Function for computing two different uncertainty
+#'   measures (simple and advanced) for Atlas Product Reproducibility. 
+#' @author M. Iturbide
+
+
+AR6.WGI.hatching <- function(delta, historical.ref = NULL, method, relative.delta = NULL, map.hatching.args = list(), ...){
+  plot.args <- list(...)
+  if (method == "advanced" & is.null(historical.ref)) stop("Provide historical.ref to use the advanced method")
+  if (is.null(map.hatching.args[["threshold"]])) map.hatching.args[["threshold"]] <- 0.5
+  if (is.null(map.hatching.args[["angle"]])) map.hatching.args[["angle"]] <- "-45"
+  if (is.null(map.hatching.args[["lwd"]])) map.hatching.args[["lwd"]] <- 0.6
+  if (is.null(map.hatching.args[["density"]])) map.hatching.args[["density"]] <- 4
+  if (is.null(map.hatching.args[["upscaling.aggr.fun"]])) map.hatching.args[["upscaling.aggr.fun"]] <- list(FUN = mean)
+  if (is.null(map.hatching.args[["condition"]])) map.hatching.args[["condition"]] <- "LT"
+  if (is.null(plot.args[["backdrop.theme"]])) plot.args[["backdrop.theme"]] <- "coastline"
+
+  simple <- suppressMessages(aggregateGrid(delta, aggr.mem = list(FUN = agreement, th = 80)))  
+
+  if (method == "simple") {
+
+    map.hatching.args[["clim"]] <- suppressMessages(climatology(simple))
+    plot.args[["sp.layout"]] <- list(do.call("map.hatching", map.hatching.args))
+
+  } else if (method == "advanced") {
+
+    sign <- suppressMessages(signal(h = historical.ref, d = delta))
+
+    advanced1 <- suppressMessages(aggregateGrid(sign, aggr.mem = list(FUN = signal.ens1, th = 66)))
+
+    advanced2.aux <- suppressMessages(aggregateGrid(sign, aggr.mem = list(FUN = signal.ens2, th = 66)))
+    advanced2.aux <- gridArithmetics(advanced2.aux, simple, operator = "+") 
+    advanced2 <- binaryGrid(advanced2.aux, condition = "GT", threshold = 0)
+
+    map.hatching.args[["clim"]] <-  suppressMessages(climatology(advanced1))
+    advanced1.hatch <- do.call("map.hatching", map.hatching.args)
+    map.hatching.args[["clim"]] <-  suppressMessages(climatology(advanced2))
+    advanced2.hatch <- do.call("map.hatching", map.hatching.args)
+    map.hatching.args[["angle"]] <-  as.character(as.numeric(map.hatching.args[["angle"]]) * -1)
+    advanced2.hatch.bis <- do.call("map.hatching", map.hatching.args)
+    plot.args[["sp.layout"]] <- list( advanced1.hatch, advanced2.hatch, advanced2.hatch.bis)
+
+  } else {
+    stop("Wrong method. Select 'simple' or 'advanced'")
+  }
+
+  if (!is.null(relative.delta)) {
+    plot.args[["grid"]] <- relative.delta
+  } else {
+    plot.args[["grid"]] <- suppressMessages(aggregateGrid(delta, aggr.mem = list(FUN = mean, na.rm = TRUE)))
+  }
+  suppressWarnings(do.call("spatialPlot", plot.args))
+}
+

_downloads/61f8c9492c42eed6fb443b73ac702769/01_index_calculation_SPI.R

@@ -0,0 +1,156 @@
+# 01_index_calculation_SPI.R
+#
+# Copyright (C) 2021 Santander Meteorology Group (http://meteo.unican.es)
+#
+# This work is licensed under a Creative Commons Attribution 4.0 International
+# License (CC BY 4.0 - http://creativecommons.org/licenses/by/4.0)
+
+#' @title Generate SPI6-12 NetCDF4 files
+#' @description Generate SPI6-12 (see table below) NetCDF4 files from CMIP5
+#'   and CMIP6 Model Outputs for Atlas Product Reproducibility
+#' @author M. Iturbide
+
+# This script builds on the climate4R framework 
+# https://github.com/SantanderMetGroup/climate4R
+
+library(devtools)
+# Data loading and writting libraries:
+library(loadeR)
+library(loadeR.2nc)
+# Index calculation libraries:
+library(transformeR)
+library(drought4R)
+
+# Function for latitudinal chunking 
+source_url("https://github.com/SantanderMetGroup/climate4R/blob/devel/R/climate4R.chunk.R?raw=TRUE")
+
+# DATA ACCESS ------------------------------------------------------------------------------------
+# Data is accessed remotely from the Santander Climate Data Service
+# Obtain a user and password at http://meteo.unican.es/udg-wiki
+loginUDG(username = "yourUser", password = "yourPassword")
+
+
+# USER PARAMETER SETTING: SELECT INDEX, SCENARIO, DATASETS, NUMBER OF CHUNKS AND OUTPUT DIRECTORY--------------
+
+
+# SELECT THE ATLAS INDEX
+# Atlas indices are: 
+# | code   | description                           | units      | script 
+# | ------ | ------------------------------------  | ---------- | --------
+# | pr     | Precipitation                         | mm         | script1_index_calculation.R
+# | pr99   | 99 percentile of all days             | mm         | script1_index_calculation.R
+# | tas    | Temperature                           | degC       | script1_index_calculation.R
+# | tasmin | mean min. temp                        | degC       | script1_index_calculation.R
+# | tasmax | mean max. temp                        | degC       | script1_index_calculation.R
+# | TNn    | minimum of min.temp                   | degC       | script1_index_calculation.R
+# | TXx    | maximum of max.temp                   | degC       | script1_index_calculation.R
+# | TX35   | max. temp avobe 35degC                | days       | script1_index_calculation.R
+# | TX40   | max. temp avobe 40degC                | days       | script1_index_calculation.R
+# | TX35bc | bias corrected max. temp avobe 35degC | days       | script1_index_calculation_bias_correction.R
+# | TX40bc | bias corrected max. temp avobe 40degC | days       | script1_index_calculation_bias_correction.R
+# | SPI6   | SPI-6                                 |            | script1_index_calculation_SPI.R
+# | SPI12  | SPI-12                                |            | script1_index_calculation_SPI.R
+# | Rx1day | maximum 1-day precipitation           | mm         | script1_index_calculation.R
+# | Rx5day | maximum 5-day precipitation           | mm         | script1_index_calculation.R
+# | DS*    | dry spell, consecutive dry days CDD   | days       | script1_index_calculation.R
+# | CDD*   | cooling degree days                   | degreedays | script1_index_calculation.R
+# | HDD    | heating degree days                   | degreedays | script1_index_calculation.R
+# | FD     | frost days                            | days       | script1_index_calculation.R
+# | T21.5  | mean temperature above 21.5degC       | days       | script1_index_calculation.R
+# Indices marked with * are not ready yet
+# Indices SPI6, SPI12 are calculated with the sript "script1_index_calculation_SPI.R"
+# Indices TX35bc, TX35bc are calculated with the script "script1_index_calculation_bias_correction.R"
+
+# Select one of SPI6 and SPI12, e.g.:
+AtlasIndex <- "SPI6"
+
+# Select scenario, e.g.:
+scenario <- "rcp85"  # possible choices are c("rcp26", "rcp45", "rcp85", "ssp126", "ssp585")
+
+# Select datasets, for the historical (datasets1) and rcp scenarios (datasets2), e.g.:
+datasets1 <- UDG.datasets("CMIP5.*historical")[["CMIP5_AR5_1run"]]
+datasets2 <-  UDG.datasets(paste0("CMIP5.*", scenario))[["CMIP5_AR5_1run"]]
+
+# # TO USE THE SAME SET OF MODELS USED IN THE ATLAS uncomment this paragraph. 
+# # Use the *.csv files of the inventories 
+# # (https://github.com/SantanderMetGroup/ATLAS/tree/devel/AtlasHub-inventory).
+# # This is the root of the ATLAS repo content:
+# root <- "https://raw.githubusercontent.com/SantanderMetGroup/ATLAS/devel/"
+# # read the first column of the desired *.csv file (e.g. CMIP5_Atlas_20191211.csv) as follows:
+# file <- "AtlasHub-inventory/CMIP5_Atlas_20191211.csv"
+# datasets <- read.csv(paste0(root, file))[,1]
+# datasets1 <- grep(datasets, pattern = "historical", value = TRUE)
+# datasets2 <- grep(datasets, pattern = scenario, value = TRUE)
+
+# Select number of chunks
+# Note: chunking sequentially splits the task into manageable data chunks to avoid memory problems
+# Chunking operates by spliting the data into a predefined number latitudinal slices (n=2 in this example).
+# Further details: https://github.com/SantanderMetGroup/climate4R/tree/master/R 
+n.chunks <- 2
+
+# Output directory where the generated results will be saved, e.g. the current one:
+out.dir <- getwd()
+
+
+# PARAMETER DEFINITION BASED ON OBJECT `AtlasIndex` and COMPUTE INDEX -------------------------------------------------
+
+# Match common datasets among scenarios
+datasets1.aux <- gsub("_historical", "", datasets1)
+datasets2.aux <- gsub(paste0("_", scenario), "", datasets2)
+
+datasets1 <- datasets1[which(datasets1.aux %in% datasets2.aux)]
+datasets2 <- datasets2[which(datasets2.aux %in% datasets1.aux)]
+
+# The climate4R function to be applied:
+spifun <- function(h, f, scale, ref.start, ref.end) {
+  hf <- bindGrid(h, f, dimension = "time")
+  speiGrid(hf, scale = scale, ref.start = ref.start, ref.end = ref.end)
+}
+
+# PARAMETER DEFINITION
+switch(AtlasIndex, 
+       SPI6 = {
+         scale <- 6
+       },
+       SPI12 = {
+         scale <- 12
+       })
+
+# COMPUTE INDEX 
+
+lapply(1:length(datasets1), function(i) {
+  di <- dataInventory(datasets1[i])
+  di2 <- dataInventory(datasets2[i])
+  if (any(names(di) %in% "pr") & any(names(di2) %in% "pr")) {
+    ch <- strsplit(di[["pr"]]$Dimensions$time$Date_range, "-")[[1]][c(1,4)]
+    years <- as.numeric(ch[1]):as.numeric(ch[2])
+    years <- years[which(years < 2101)]
+    years <- years[which(years > 1949)]
+
+    # prepare output directory chain
+    fol <- paste0(out.dir, "/", AtlasIndex,"/raw/", strsplit(datasets2[i], "_")[[1]][2])
+    if (!dir.exists(fol)) dir.create(fol) 
+    fol <- paste0(fol, "/", strsplit(datasets2[i], "_")[[1]][4])
+    if (!dir.exists(fol)) dir.create(fol)
+
+    # calculate index
+    index <- climate4R.chunk(n.chunks = n.chunks,
+                             C4R.FUN.args = list(FUN = "spifun", h = list(dataset = datasets1[i], var = "pr", years = 1971:2005), 
+                                                 f = list(dataset = datasets2[i], var = "pr", years = years), 
+                                                 scale = scale,
+                                                 ref.start = c(1971, 1), ref.end = c(2010, 12)),
+                             loadGridData.args = list(aggr.m = "sum"))
+    index[["Variable"]][["varName"]] <- AtlasIndex
+    index <- redim(index, drop = TRUE)
+    message(".........", i)
+
+    # WRITE .nc FILES (write each year in separate files)
+
+    lapply(years, function(x){
+      print(paste0("------", i, "-----------------------------------------------"))
+      index.x <- subsetGrid(index, years = x)
+      grid2nc(index.x, NetCDFOutFile = paste0(fol, "/", datasets2[i], "_" , AtlasIndex, "_", x, ".nc4"))
+    })
+  }
+})
+