## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----monte carlo 1, message=FALSE, warning=FALSE, eval=FALSE------------------ # library(gamlss) # library(gamlss.dist) # library(MuMIn) # library(spsUtil) # # Designing the selection function based on the AICc # select.model <- function(rain.week) { # best <- matrix(NA, 1, 2) # colnames(best) <- c("NonLinear", "Linear") # t.gam <- quiet(gamlss( # rain.week ~ 1, # family = GA, # mu.link = "log", # sigma.link = "log" # )) # t.gam.ns10 <- quiet(gamlss( # rain.week ~ time, # family = GA, # mu.link = "log", # sigma.link = "log" # )) # t.gam.ns01 <- quiet(gamlss( # rain.week ~ 1, # sigma.formula = ~time, # family = GA, # mu.link = "log", # sigma.link = "log" # )) # t.gam.ns11 <- quiet(gamlss( # rain.week ~ time, # sigma.formula = ~time, # family = GA, # mu.link = "log", # sigma.link = "log" # )) # t.gam.ns20 <- quiet(gamlss( # rain.week ~ time + I(time^2), # family = GA, # mu.link = "log", # sigma.link = "log" # )) # t.gam.ns02 <- quiet(gamlss( # rain.week ~ 1, # family = GA, # mu.link = "log", # sigma.formula = ~ time + # I(time^2), # sigma.link = "log" # )) # t.gam.ns21 <- quiet(gamlss( # rain.week ~ time + I(time^2), # sigma.formula = ~time, # family = GA, # mu.link = "log", # sigma.link = "log" # )) # t.gam.ns12 <- quiet(gamlss( # rain.week ~ time, # sigma.formula = ~ time + I(time^2), # family = GA, # mu.link = "log", # sigma.link = "log" # )) # t.gam.ns22 <- quiet(gamlss( # rain.week ~ time + I(time^2), # sigma.formula = ~ time + I(time^2), # family = GA, # mu.link = "log", # sigma.link = "log" # )) # t.gam.ns30 <- quiet(gamlss( # rain.week ~ time + # I(time^2) + # I(time^3), # family = GA, # mu.link = "log", # sigma.link = "log" # )) # t.gam.ns03 <- quiet(gamlss( # rain.week ~ 1, # family = GA, # mu.link = "log", # sigma.formula = ~ time + # I(time^2) + # I(time^3), # sigma.link = "log" # )) # t.gam.ns31 <- quiet(gamlss( # rain.week ~ time + # I(time^2) + # I(time^3), # sigma.formula = ~time, # family = GA, # mu.link = "log", # sigma.link = "log" # )) # t.gam.ns13 <- quiet(gamlss( # rain.week ~ time, # sigma.formula = ~ time + # I(time^2) + # I(time^3), # family = GA, # mu.link = "log", # sigma.link = "log" # )) # t.gam.ns32 <- quiet(gamlss( # rain.week ~ time + # I(time^2) + # I(time^3), # sigma.formula = ~ time + # I(time^2), # family = GA, # mu.link = "log", # sigma.link = "log" # )) # t.gam.ns23 <- quiet(gamlss( # rain.week ~ time + # I(time^2), # sigma.formula = ~ time + # I(time^2) + # I(time^3), # family = GA, # mu.link = "log", # sigma.link = "log" # )) # t.gam.ns33 <- quiet(gamlss( # rain.week ~ time + # I(time^2) + # I(time^3), # sigma.formula = ~ time + # I(time^2) + # I(time^3), # family = GA, # mu.link = "log", # sigma.link = "log" # )) # # Selection among all 16 models # best[1, 1] <- which.min(list( # MuMIn::AICc(t.gam, k = k), # MuMIn::AICc(t.gam.ns10, k = k), # MuMIn::AICc(t.gam.ns01, k = k), # MuMIn::AICc(t.gam.ns11, k = k), # MuMIn::AICc(t.gam.ns20, k = k), # MuMIn::AICc(t.gam.ns02, k = k), # MuMIn::AICc(t.gam.ns21, k = k), # MuMIn::AICc(t.gam.ns12, k = k), # MuMIn::AICc(t.gam.ns22, k = k), # MuMIn::AICc(t.gam.ns30, k = k), # MuMIn::AICc(t.gam.ns03, k = k), # MuMIn::AICc(t.gam.ns31, k = k), # MuMIn::AICc(t.gam.ns13, k = k), # MuMIn::AICc(t.gam.ns32, k = k), # MuMIn::AICc(t.gam.ns23, k = k), # MuMIn::AICc(t.gam.ns33, k = k) # )) # # Selection among the stationary and linear models (models 1 to 4) # best[1, 2] <- which.min(list( # MuMIn::AICc(t.gam, k = k), # MuMIn::AICc(t.gam.ns10, k = k), # MuMIn::AICc(t.gam.ns01, k = k), # MuMIn::AICc(t.gam.ns11, k = k) # )) # # return(best) # } # library(gamlss) # library(gamlss.dist) # library(MuMIn) # library(spsUtil) # # Monte Carlo Simulation with three Sample Sizes # sample_sizes <- c(30, 60, 90) # Define the sample sizes # ns <- 10000 # Number of simulations # k.seq <- seq(from = 2, to = 7, by = 0.2) # penalty <- length(k.seq) # # results_list <- list() # selected_nonlinear_list <- list() # selected_linear_list <- list() # # for (n in sample_sizes) { # time <- 1L:n # rain.week <- matrix(NA, nrow = n, ncol = ns) # # sigma0 <- 0.22 # mu0 <- 809 # slope.mu <- 0 # trend-free # slope.sigma <- 0 # trend-free # sigma <- sigma0 + slope.sigma * time # mu <- mu0 + slope.mu * time # # # Simulating precipitation data # rain.week <- replicate(ns, sapply(1:n, function(t) rGA(1, mu = mu[t], sigma = sigma[t]))) # # result <- matrix(NA, nrow = ns, ncol = 2 * penalty) # selected.nonlinear <- matrix(NA, 16, penalty) # selected.linear <- matrix(NA, 4, penalty) # # # Calculating loop for each k # for (j in seq_along(k.seq)) { # k <- k.seq[j] # result1 <- apply(rain.week, 2, select.model) # result[, (2 * j - 1):(2 * j)] <- t(result1) # # for (model in 1:16) { # selected.nonlinear[model, j] <- 100 * (sum(result[, 2 * j - 1] == model) / ns) # if (model <= 4) { # selected.linear[model, j] <- 100 * (sum(result[, 2 * j] == model) / ns) # } # } # } # rownames(selected.nonlinear) <- paste0("Model", 1:16) # rownames(selected.linear) <- paste0("Model", 1:4) # colnames(selected.nonlinear) <- as.character(k.seq) # colnames(selected.linear) <- as.character(k.seq) # # results_list[[as.character(n)]] <- result # selected_nonlinear_list[[as.character(n)]] <- selected.nonlinear # selected_linear_list[[as.character(n)]] <- selected.linear # } # # for (n in sample_sizes) { # cat("\nSample Size:", n, "\n") # cat("\nSelected Nonlinear Models:\n") # print(selected_nonlinear_list[[as.character(n)]]) # cat("\nSelected Linear Models:\n") # print(selected_linear_list[[as.character(n)]]) # } ## ----monte carlo 2, message=FALSE, warning=FALSE, eval=FALSE------------------ # ## Monte Carlo Simulation # n <- 90 # ns <- 10000 # time <- 1L:n # # rain.week <- matrix(NA, nrow = n, ncol = ns) # sigma0 <- 0.22 # mu0 <- 809 # slope.mu <- (-0.004 * mu0) # super-imposed trend # slope.sigma <- 0 # trend-free # sigma <- sigma0 + slope.sigma * time # mu <- mu0 + slope.mu * time # k.seq <- seq(from = 2, to = 7, by = 0.2) # penalty <- length(k.seq) # # result <- matrix(NA, nrow = ns, ncol = 2 * penalty) # selected.nonlinear <- matrix(NA, 16, penalty) # selected.linear <- matrix(NA, 4, penalty) # # simulating again # rain.week <- replicate(ns, sapply(1:n, function(t) rGA(1, mu = mu[t], sigma = sigma[t]))) # # for (j in seq_along(k.seq)) { # k <- k.seq[j] # result1 <- apply(rain.week, 2, select.model) # result[, (2 * j - 1):(2 * j)] <- t(result1) # # for (model in 1:16) { # selected.nonlinear[model, j] <- 100 * (sum(result[, 2 * j - 1] == model) / ns) # if (model <= 4) { # selected.linear[model, j] <- 100 * (sum(result[, 2 * j] == model) / ns) # } # } # } # # rownames(selected.nonlinear) <- paste0("Model", 1:16) # rownames(selected.linear) <- paste0("Model", 1:4) # colnames(selected.nonlinear) <- as.character(k.seq) # colnames(selected.linear) <- as.character(k.seq) # # print(selected.nonlinear) # print(selected.linear) ## ----get-oxford-rain---------------------------------------------------------- # The {curl} library is used to handle possible download issues with timeouts library(curl) h <- new_handle() handle_setopt(h, timeout = 360L) rain_file <- file.path(tempdir(), "oxford_rain.csv") curl::curl_download( url = "https://figshare.com/ndownloader/files/21950895", destfile = rain_file, mode = "wb", quiet = FALSE, handle = h ) Oxford.1815 <- read.csv(rain_file) Oxford <- Oxford.1815[which(Oxford.1815$YYY >= 1827), ] # Rainfall records, including traces (NA), started in 1827 summary(Oxford) ## ----replace-gaps------------------------------------------------------------- # create a new vector of rain to work with OxfordRain <- Oxford$Rainfall.mm.1.dpl.no.traces # set to 0 OxfordRain[is.na(OxfordRain)] <- 0 # ensure no NAs summary(OxfordRain) ## ----Case Study 1------------------------------------------------------------- library(SPIChanges) rainTS4 <- TSaggreg(daily.rain = OxfordRain, start.date = "1827-01-01", TS = 4) ### Fit All Models: Stationary, Linear and Nonlinear Oxford.Changes <- SPIChanges(rain.at.TS = rainTS4, only.linear = "No") ### Fit Only the Stationary and Linear Models Oxford.Changes.linear <- SPIChanges(rain.at.TS = rainTS4, only.linear = "Yes") ## ----selected models, fig.width=8, fig.height=5------------------------------- eixo.x.Month <- Oxford.Changes$model.selection[, 1] eixo.y.quasiWeek <- Oxford.Changes$model.selection[, 2] valores <- Oxford.Changes$model.selection[, 3] dados <- cbind(eixo.x.Month, eixo.y.quasiWeek, valores) df <- as.data.frame(dados) library(ggplot2) df$valores <- as.factor(df$valores) fig1 <- ggplot(df) + aes(x = eixo.x.Month, y = eixo.y.quasiWeek, fill = valores) + geom_tile() + scale_fill_manual(values = c( "1" = "#D3D3D3", "2" = "#ADD8E6", "3" = "#87CEEB", "4" = "#4682B4", "6" = "#FA8072", "11" = "#FF8C00" )) + theme_bw() + labs(x = "Months", y = "quasi-week", fill = NULL, title = "only.linear = 'No'", caption = ("Figure 1. Gamma-based models selected by the SPIChanges package. The plot corresponds to setting \n the `only.linear` argument of the SPIChanges() function to `No`. Monthly precipitation series recorded \n at the Radcliffe Observatory site in Oxford-UK (January 1827 to January 2020).")) + scale_x_continuous(breaks = 1:12, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")) + geom_text(aes(label = valores), color = "black", size = 3, fontface = "bold") + theme( strip.text = element_text(color = "black", face = "bold"), axis.text.x = element_text(size = 10, color = "black"), axis.text.y = element_text(size = 10, color = "black"), axis.title.x = element_text(size = 10, face = "bold", color = "black"), axis.title.y = element_text(size = 10, face = "bold", color = "black"), plot.title = element_text(face = "bold", size = 14, color = "black"), plot.caption = element_text(hjust = 0, size = 10, color = "black", lineheight = 1.0), plot.margin = margin(10, 10, 20, 10) ) fig1 ## ----year-to-year changes, fig.width=10, fig.height=5------------------------- library(dplyr) library(ggplot2) library(patchwork) library(purrr) # Define constants eixo.x.years <- 1827L:2019L # Define conditions conditions <- list( jan1 = list(Month = 1, quasiWeek = 1, col = 4), jan2 = list(Month = 1, quasiWeek = 2, col = 4), jan3 = list(Month = 1, quasiWeek = 3, col = 4), jan4 = list(Month = 1, quasiWeek = 4, col = 4), jul3 = list(Month = 7, quasiWeek = 3, col = 4), jul4 = list(Month = 7, quasiWeek = 4, col = 4), dec4mu = list(Month = 12, quasiWeek = 4, col = 4), jun4 = list(Month = 6, quasiWeek = 4, col = 5), jul1 = list(Month = 7, quasiWeek = 1, col = 5), jul2 = list(Month = 7, quasiWeek = 2, col = 5), sep3 = list(Month = 9, quasiWeek = 3, col = 5), dec2 = list(Month = 12, quasiWeek = 2, col = 5), dec4 = list(Month = 12, quasiWeek = 4, col = 5), mar2 = list(Month = 3, quasiWeek = 2, col = 5), mar3 = list(Month = 3, quasiWeek = 3, col = 5), sep1 = list(Month = 9, quasiWeek = 1, col = 5) ) # Extract data using purrr::map results <- map(names(conditions), ~ { cond <- conditions[[.x]] Oxford.Changes$Statistics %>% filter(Month == cond$Month, quasiWeek == cond$quasiWeek) %>% pull(cond$col) }) %>% set_names(names(conditions)) # Assign results to variables list2env(results, envir = .GlobalEnv) # Adjust jan4 jan4 <- head(jan4, -1) # Combine data into a single dataframe df2 <- data.frame( ano = rep(eixo.x.years, length.out = length(unlist(results))), valores = unlist(results), line = rep(names(results), lengths(results)), letra = rep(c("a", "b", "c"), times = c( sum(lengths(results[1:7])), sum(lengths(results[8:13])), sum(lengths(results[14:16])) )) ) # Recode `line` for better readability df2 <- df2 %>% mutate(line = recode(line, "jan1" = "Jan (1st quasi-week)", "jan2" = "Jan (2nd quasi-week)", "jan3" = "Jan (3rd quasi-week)", "jan4" = "Jan (4th quasi-week)", "jul1" = "Jul (1st quasi-week)", "jul2" = "Jul (2nd quasi-week)", "jul3" = "Jul (3rd quasi-week)", "jul4" = "Jul (4th quasi-week)", "jun4" = "Jun (4th quasi-week)", "sep1" = "Sep (1st quasi-week)", "sep3" = "Sep (3rd quasi-week)", "dec2" = "Dec (2nd quasi-week)", "dec4" = "Dec (4th quasi-week)", "mar2" = "Mar (2nd quasi-week)", "mar3" = "Mar (3rd quasi-week)" )) # Define colors cores <- c( "Jan (1st quasi-week)" = "blue", "Jan (2nd quasi-week)" = "gold", "Jan (3rd quasi-week)" = "firebrick", "Jan (4th quasi-week)" = "black", "Jul (3rd quasi-week)" = "cyan", "Jul (4th quasi-week)" = "red", "Dec (4th quasi-week)" = "#1f77b4", "Jun (4th quasi-week)" = "#7f7f7f", "Jul (1st quasi-week)" = "#008080", "Jul (2nd quasi-week)" = "#8B4513", "Sep (3rd quasi-week)" = "#4B0082", "Dec (2nd quasi-week)" = "#ffbb78", "Mar (2nd quasi-week)" = "#9467bd", "Mar (3rd quasi-week)" = "#bcbd22", "Sep (1st quasi-week)" = "#c49c94" ) # Function to create plots (updated with legend.position.inside) create_plot <- function(data, y_limits, y_breaks, y_lab, title) { ggplot(data) + aes(x = ano, y = valores, colour = line) + geom_line(linewidth = 1.0) + scale_color_manual(values = cores) + theme_bw() + theme( legend.position.inside = c(0.9, 0.3), # Updated here legend.title = element_blank(), legend.key.size = unit(0.1, "cm"), legend.box.spacing = unit(0.1, "cm"), legend.direction = "horizontal", legend.box = "horizontal", legend.background = element_blank(), legend.text = element_text(size = 9), legend.key.width = unit(1, "cm"), axis.title.y = element_text(size = 9, face = "bold"), strip.text = element_text(face = "bold", size = 12), plot.margin = margin(3, 3, 3, 3), axis.title.x = element_blank(), axis.text.x = element_blank() ) + scale_x_continuous(breaks = seq(1827, 2019, by = 25)) + scale_y_continuous(limits = y_limits, breaks = y_breaks) + ylab(y_lab) + ggtitle(title) + guides(colour = guide_legend(ncol = 2)) } # Create plots plot_a <- create_plot(df2 %>% filter(letra == "a"), c(0, 100), seq(0, 100, by = 20), "mean parameter \n (μ) GAMLSS", "A") plot_b <- create_plot(df2 %>% filter(letra == "b"), c(0, 1), seq(0, 1, by = 0.2), "dispersion parameter \n (δ) GAMLSS", "B") plot_c <- create_plot(df2 %>% filter(letra == "c"), c(0, 2), seq(0, 2, by = 0.4), "dispersion parameter \n (δ) GAMLSS", "C") + labs(caption = "Figure 2. Year-to-year changes in the dispersion parameter of gamma-based models estimated from monthly precipitation records of the Radcliffe Observatory site in Oxford-UK. (A) linear changes in the μ parameter (B) linear changes in the δ parameter (C) nonlinear changes in the δ parameter") + theme(axis.title.x = element_text(size = 10, face = "bold"), plot.caption = element_text(hjust = 0, size = 10, color = "black", lineheight = 1.0)) # Combine plots fig2 <- plot_a + plot_b + plot_c + plot_layout(ncol = 1, heights = c(1.3, 1.3, 1)) fig2 ## ----changes in drought frequency, fig.width=8, fig.height=7------------------ library(ggplot2) library(tidyr) library(dplyr) library(patchwork) Oxford.Changes$Changes.Freq.Drought <- as.data.frame(Oxford.Changes$Changes.Freq.Drought) # Define quasi-weeks and labels quasi_weeks <- list( "1st - Jan" = c(1,1), "2nd - Jan" = c(1,2), "3rd - Jan" = c(1,3), "4th - Jan" = c(1,4), "2nd - Mar" = c(3,2), "3rd - Mar" = c(3,3), "4th - Jun" = c(6,4), "1st - Jul" = c(7,1), "2nd - Jul" = c(7,2), "3rd - Jul" = c(7,3), "4th - Jul" = c(7,4), "1st - Sep" = c(9,1), "3rd - Sep" = c(9,3), "2nd - Dec" = c(12,2), "3rd - Dec" = c(12,4) ) # Function to extract data drought_data <- function(index) { sapply(quasi_weeks, function(qw) { Oxford.Changes$Changes.Freq.Drought %>% filter(Month == qw[1], quasiWeek == qw[2]) %>% pull(index) }) } # Create data frame df3 <- data.frame( x = names(quasi_weeks), moderate = drought_data(7), severe = drought_data(8), extreme = drought_data(9), stat = drought_data(5), nonstat = drought_data(6) ) %>% pivot_longer(-x, names_to = "category", values_to = "value") # Define colors color_map <- c( "moderate" = "#00243F", "severe" = "#517C9D", "extreme" = "#36A7C1", "stat" = "#1c340a", "nonstat" = "#3ba551" ) # Generate plot function generate_plot <- function(data, y_label, title, limits, colors, legend_labels) { ggplot(data, aes(x = x, y = value, fill = category)) + geom_bar(stat = "identity", position = "dodge") + scale_fill_manual(values = colors, labels = legend_labels) + scale_y_continuous(limits = limits) + labs(title = title, x = NULL, y = y_label) + theme_bw() + theme( axis.text.x = element_text(size = 10, angle = 45, hjust = 1, color = "black"), axis.text.y = element_text(size = 10, color = "black"), axis.title.x = element_text(size = 10, face = "bold", color = "black"), axis.title.y = element_text(size = 9, face = "bold", color = "black"), legend.position = "top", legend.title = element_blank(), legend.text = element_text(size = 10, color = "black"), legend.key.width = unit(1, "cm"), legend.key.size = unit(0.4, "cm"), plot.margin = margin(10, 10, 10, 10) ) } # Create plots plot_3a <- generate_plot( df3 %>% filter(category %in% c("moderate", "severe", "extreme")), "Percentual change (%)", "A", c(-20, 20), color_map[c("moderate", "severe", "extreme")], c("Change Moderate", "Change Severe", "Change Extreme") ) plot_3b <- generate_plot( df3 %>% filter(category %in% c("stat", "nonstat")), "mm", "B", c(0, 100), color_map[c("stat", "nonstat")], c("Stationary Normal Rain", "NonStationary Normal Rain") ) + labs(caption = "Figure 3. Year-to-year changes in drought frequency estimated from quasi-weekly precipitation records of the Radcliffe Observatory site in Oxford-UK. (A) Percentual changes in the frequency of moderate, severe, and extreme droughts. (B) Changes in the normal precipitation amount for stationary and nonstationary normality assumptions.") + theme(axis.title.x = element_text(size = 10, face = "bold"), plot.caption = element_text(hjust = 0, size = 10, color = "black", lineheight = 1.0)) # Combine plots Fig3 <- plot_3a + plot_3b + plot_layout(ncol = 1, heights = c(1.5, 1.2)) Fig3 ## ----Monitoring drought in Oxford-UK, 2019------------------------------------ Table1 <- Oxford.Changes$data.week[which(Oxford.Changes$data.week[, 1] == 2019), ] # Table 1. Output data field /$data.week obtained by applying the SPIChanges package at the 4-quasi week time scale to the precipitation series recorded at the Radcliffe Observatory site in Oxford-UK (January 1827 to January 2020): rain.at.TS is the precipitation amount accumulated at the time scale specified by the users, SPI is the Standardized Precipitation Index, Exp.Acum.Prob and Actual.Acum.Prob are the expected frequency of the SPI estimates calculated under the stationary and under non-stationary approaches, respectively and, ChangeFreq is the changes in the frequency of dry events (SPI < 0). It is the difference between Actual.Acum.Prob and Exp.Acum.Prob. Table1 ## ----selected linear models, fig.width=8, fig.height=5------------------------ eixo.x.Month <- Oxford.Changes.linear$model.selection[, 1] eixo.y.quasiWeek <- Oxford.Changes.linear$model.selection[, 2] valores <- Oxford.Changes.linear$model.selection[, 3] dados <- cbind(eixo.x.Month, eixo.y.quasiWeek, valores) df <- as.data.frame(dados) library(ggplot2) df$valores <- as.factor(df$valores) fig4 <- ggplot(df) + aes(x = eixo.x.Month, y = eixo.y.quasiWeek, fill = valores) + geom_tile() + scale_fill_manual(values = c( "1" = "#D3D3D3", "2" = "#ADD8E6", "3" = "#87CEEB", "4" = "#4682B4" )) + theme_bw() + labs(x = "Months", y = "quasi-week", fill = NULL, title = "only.linear = 'Yes'", caption = "Figure 4. Gamma-based models selected by the SPIChanges package. The plot corresponds to \n setting the `only.linear` argument of the SPIChanges() function to `Yes`. Monthly precipitation series \n recorded at the Radcliffe Observatory site in Oxford-UK (January 1827 to January 2020).") + scale_x_continuous(breaks = 1:12, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")) + geom_text(aes(label = valores), color = "black", size = 3, fontface = "bold") + theme( strip.text = element_text(color = "black", face = "bold"), axis.text.x = element_text(size = 10, color = "black"), axis.text.y = element_text(size = 10, color = "black"), axis.title.x = element_text(size = 10, face = "bold", color = "black"), axis.title.y = element_text(size = 10, face = "bold", color = "black"), plot.title = element_text(face = "bold", size = 14, color = "black"), plot.caption = element_text(hjust = 0, size = 10, color = "black", lineheight = 1.0), plot.margin = margin(10, 10, 20, 10) ) fig4 ## ----Monitoring drought in Oxford-UK, 2019, only linear models---------------- Table2 <- Oxford.Changes.linear$data.week[which(Oxford.Changes.linear$data.week[, 1] == 2019), ] # Table 2. Output data field /$data.week obtained by applying the SPIChanges package at the 4-quasi week time scale to the precipitation series recorded at the Radcliffe Observatory site in Oxford-UK (January 1827 to January 2020): rain.at.TS is the precipitation amount accumulated at the time scale specified by the users, SPI is the Standardized Precipitation Index, Exp.Acum.Prob and Actual.Acum.Prob are the expected frequency of the SPI estimates calculated under the stationary and under non-stationary approaches, respectively and, ChangeFreq is the changes in the frequency of dry events (SPI < 0). It is the difference between Actual.Acum.Prob and Exp.Acum.Prob. Table2 ## ----get-CPC-rain, eval=FALSE------------------------------------------------- # # this chunk may take a long time to get daily precipitation data for entire Brazil (1980-2024) # # Load necessary libraries # library(ncdf4) # library(curl) # library(SPIChanges) # # n <- nrow(lonlat) # # Configure curl with a longer timeout # h <- new_handle() # handle_setopt(h, timeout = 10800L) # # download_and_process <- function(year, lonlat) { # # Define file paths and URLs # rain_file <- file.path(tempdir(), paste0("cpcdata", year, ".nc")) # download_url <- paste0("https://psl.noaa.gov/thredds/fileServer/Datasets/cpc_global_precip/precip.", year, ".nc") # # Download the files # curl::curl_download( # url = download_url, # destfile = rain_file, # mode = "wb", # quiet = FALSE, # handle = h # ) # # # Open and process the netCDF file # cep.prec <- nc_open(rain_file) # lons <- ncvar_get(cep.prec, "lon") # lats <- ncvar_get(cep.prec, "lat") # prate <- ncvar_get(cep.prec, "precip") # mm/day # NN <- dim(prate)[3] # # Flip latitude and precipitation order # prate <- prate[, rev(seq_len(dim(prate)[2])), 1:NN] # lats <- rev(lats) # # Initialize matrix for the year # pdata <- matrix(NA, NN, n) # for (i in seq_len(n)) { # longitude <- lonlat[i, 1] + 360 # range 0, 360 # latitude <- lonlat[i, 2] # range -90, 90 # pdata[, i] <- prate[ # which.min((lons - longitude)^2), # which.min((lats - latitude)^2), # ] # } # # nc_close(cep.prec) # return(pdata) # } # # # Process data for multiple years # years <- 1980:2024 # pdata_list <- lapply(years, download_and_process, lonlat = lonlat) # # # Combine all years' data # pdata1 <- do.call(rbind, pdata_list) ## ----aggregate-CPC-rain, eval=FALSE------------------------------------------- # library(SPIChanges) # # This chunk takes a long time to aggregate precipitation data for entire Brazil (1980-2024) # pdata1[is.na(pdata1)] <- 0 # replacing NA by zero. # TS12 <- TSaggreg(daily.rain = pdata1[, 1], start.date = "1980-01-01", TS = 12) # N <- nrow(TS12) # rain.at.TS <- matrix(NA, n, 5) # rainTS12.bank <- matrix(NA, N, (n + 3)) # rainTS12.bank[, 1:4] <- as.matrix(TS12) # for (i in 2:n) { # TS12 <- TSaggreg(daily.rain = pdata1[, i], start.date = "1980-01-01", TS = 12) # rainTS12.bank[, (i + 3)] <- TS12[, 4] # } ## ----SPIChanges for entire Brazil, eval=FALSE--------------------------------- # # trying to speed up the things # library(foreach) # library(doParallel) # library(SPIChanges) # numCores <- detectCores() - 1 # Use one less core than the total available # cl <- makeCluster(numCores) # registerDoParallel(cl) # rain <- matrix(NA, N, 4) # rain[, 1:3] <- as.matrix(rainTS12.bank[, 1:3]) # final <- ncol(rainTS12.bank) # n <- nrow(lonlat) # Map <- matrix(NA, n, 18) # model <- matrix(NA, n, 3) # Map[, 1:2] <- as.matrix(lonlat) # model[, 1:2] <- as.matrix(lonlat) # changes <- matrix(NA, 1, 16) # # Perform parallel processing # results <- foreach(i = 4:final, .combine = rbind, .packages = "SPIChanges") %dopar% { # rain[, 4] <- rainTS12.bank[, i] # Local.Changes <- SPIChanges(rain.at.TS = rain, only.linear = "no") # changes[1, 1:3] <- Local.Changes$Changes.Freq.Drought[which(Local.Changes$Changes.Freq.Drought[, 1] == 2 & # Local.Changes$Changes.Freq.Drought[, 2] == 4), 7:9] # changes[1, 4] <- Local.Changes$model.selection[which(Local.Changes$model.selection[, 1] == 2 & # Local.Changes$model.selection[, 2] == 4), 3] # changes[1, 5:7] <- Local.Changes$Changes.Freq.Drought[which(Local.Changes$Changes.Freq.Drought[, 1] == 5 & # Local.Changes$Changes.Freq.Drought[, 2] == 4), 7:9] # changes[1, 8] <- Local.Changes$model.selection[which(Local.Changes$model.selection[, 1] == 5 & # Local.Changes$model.selection[, 2] == 4), 3] # changes[1, 9:11] <- Local.Changes$Changes.Freq.Drought[which(Local.Changes$Changes.Freq.Drought[, 1] == 8 & # Local.Changes$Changes.Freq.Drought[, 2] == 4), 7:9] # changes[1, 12] <- Local.Changes$model.selection[which(Local.Changes$model.selection[, 1] == 8 & # Local.Changes$model.selection[, 2] == 4), 3] # changes[1, 13:15] <- Local.Changes$Changes.Freq.Drought[which(Local.Changes$Changes.Freq.Drought[, 1] == 11 & # Local.Changes$Changes.Freq.Drought[, 2] == 4), 7:9] # changes[1, 16] <- Local.Changes$model.selection[which(Local.Changes$model.selection[, 1] == 11 & # Local.Changes$model.selection[, 2] == 4), 3] # return(changes) # } # # Combine results into the Map matrix # a <- 1 # for (i in 1:nrow(results)) { # Map[a, 3:18] <- results[i, ] # a <- a + 1 # } # # # Assign column names and save the results # colnames(Map) <- c( # "lon", "lat", "SummerModerate", "SummerSevere", "SummerExtreme", "SummerModel", # "AutomnModerate", "AutomnSevere", "AutomnExtreme", "AutomnModel", # "WinterModerate", "WinterSevere", "WinterExtreme", "WinterModel", # "SpringModerate", "SpringSevere", "SpringExtreme", "SpringModel" # ) # # Stop the cluster # stopCluster(cl) # head(Map) ## ----map, fig.width=8, fig.height=5------------------------------------------- library(ggplot2) library(sf) library(RColorBrewer) library(tidyr) library(dplyr) library(archive) library(SPIChanges) # Loading data rar_url <- "https://github.com/gabrielblain/Grid_Brazil/raw/main/regioes_2010.rar" temp_file <- tempfile(fileext = ".rar") temp_dir <- tempdir() download.file(rar_url, temp_file, mode = "wb") archive_extract(temp_file, dir = temp_dir) shp_path <- file.path(temp_dir, "regioes_2010.shp") # Adjust if files are extracted into a subdirectory limitere <- st_read(shp_path) # Tidying the data df <- Map %>% select(lat, lon, SpringSevere, AutomnSevere, WinterSevere, SummerSevere) %>% pivot_longer( cols = c(SpringSevere, AutomnSevere, WinterSevere, SummerSevere), names_to = "Station", values_to = "Change" ) df_ok <- df %>% mutate( Station = recode(Station, "SpringSevere" = "Spring - Severe", "SummerSevere" = "Summer - Severe", "AutomnSevere" = "Autumm - Severe", "WinterSevere" = "Winter - Severe" ) ) # Turning into sf limitere_sf <- st_as_sf(limitere) dados.sf <- st_as_sf(df_ok, coords = c("lon", "lat"), crs = 4326) # Fixing possible errors limitere_sf <- st_make_valid(limitere_sf) # Finding centroids limitere_centroids <- st_centroid(limitere_sf) # Adding siglas (annacron) limitere_centroids$sigla <- limitere$sigla # Finding coordinates coords <- st_coordinates(limitere_centroids) limitere_centroids$lon <- coords[, 1] limitere_centroids$lat <- coords[, 2] dados.sf$Station <- factor( dados.sf$Station, levels = c( "Spring - Severe", "Summer - Severe", "Autumm - Severe", "Winter - Severe" ) ) cores_roxo_branco_vermelho <- colorRampPalette(c("purple", "white", "red"))(n = 100) cores5 <- colorRampPalette(c("purple", "white", "red"))(7) # The Map map <- ggplot() + geom_sf(data = dados.sf, aes(color = Change), shape = 15, size = 3) + geom_sf(data = limitere_sf, fill = NA, color = "black", size = 1) + scale_colour_gradient2( low = "#6B238E", mid = "white", high = "red", midpoint = 0, limits = c(-30, 80) ) + theme_light() + labs(color = "Change (%)") + facet_wrap(vars(Station)) + theme( title = element_text(size = 12, face = "bold", color = "black"), text = element_text(size = 12, color = "black"), axis.text.x = element_text(size = 9, color = "black"), axis.text.y = element_text(size = 9, color = "black"), legend.position = "right", legend.direction = "vertical", legend.key.height = unit(.6, "cm"), legend.key.width = unit(0.5, "cm"), legend.title = element_text(size = 11, face = "bold"), legend.text = element_text(size = 11), axis.title.x = element_blank(), # Remove o rotulo do eixo x axis.title.y = element_blank(), # Remove o rotulo do eixo y strip.text = element_text(size = 11, face = "bold", color = "black") ) map <- map + geom_text( data = limitere_centroids, aes(x = lon, y = lat, label = sigla), color = "black", size = 3, fontface = "bold", position = position_nudge(y = 0.01), show.legend = FALSE ) print(map)