SPIChanges package: Monte Carlo Experiments and Case Studies

The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

Case Studies Applications

To further validate this package, we also applied the package’s functions to two distinct datasets. The first (Dataset 1) is from a long-running weather station in the UK, located at 51.7608°N and 1.2639°W (Radcliffe Observatory site in Oxford, covering the period from January 1827 to January 2020) [Burt2019]. These invaluable precipitation records are available at https://doi.org/10.6084/m9.figshare.11956239.v1 (Burt 2020) under a Creative Commons 4.0 License via FigShare. The percentage of missing data is less than 1%. Following the approach of Wu et al. (2006), we replaced all gaps in the series with zero, as this is the most probable value for a single day. We used this dataset to provide detailed information on the package’s outputs.

The second dataset (Dataset 2) comprises data from the CPC Global Unified Gauge-Based Analysis of Daily Precipitation (CPC data; NOAA/PSL) https://psl.noaa.gov/data/gridded/data.cpc.globalprecip.html for the entire country of Brazil. Brazil, the largest tropical country in the world, features a wide range of climates, from arid regions in the Northeast to equatorial climates, such as the Amazon Rainforest.

Case Study 1 – Oxford Rainfall

Analyse Data Looking for Changes

Apply the {SPIChanges} package’s SPIChanges() to daily precipitation data (OxfordRain) derived from Burt (2020).

library(SPIChanges)
data(OxfordRain)
rainTS4 <- TSaggreg(daily.rain = OxfordRain, start.date = "1827-01-01", TS = 4)
#> Done. Just ensure the last quasi-week is complete.
#>   The last day of your series is 31 and TS is 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")

Interpreting the Output Data Fields from the “Oxford.Changes”

We calculated the data frame object “Oxford.Changes” with the argument only.linear set to “No”. In this case, the SPIChanges() function considers 16 gamma-based models to describe potential changes in precipitation patterns in Oxford, UK. The function uses the second-order Akaike Information Criterion to select the best model among the 16 candidates and indicate whether, and how, the frequency of droughts has changed over time in Oxford, UK. Below we describe the four output data fields of the “Oxford.Changes” data frame object.

$model.selection

The $model.selection output field identifies the gamma-based models selected for each quasi-weekly period in Oxford, UK. This analysis reveals the following insights:

Stationary Models:
For 33 out of 48 quasi-weekly periods, the stationary gamma distribution (Model 1) was the best fitting model. This suggests that the probability of SPI values remained constant between 1827 and 2020 for these periods.
Non-Stationary Models:
- Model 2: Selected for six quasi-weekly periods (January and July). This model describes a linear change in the μ parameter, corresponding to the mean of the precipitation frequency distribution.
- Model 3: Selected for five quasi-weekly periods (June, July, September, and December). This model indicates a linear change in the δ parameter, representing the coefficient of variation, which suggests changes in dispersion while the mean remained constant.
- Model 4: Selected for the 4th quasi-week of December. This model shows linear changes in both μ and δ.
Nonlinear Models:
- Model 6: Selected for the 1st quasi-week of September, indicating a quadratic change in the δ parameter.
- Model 11: Selected for the 2nd and 3rd quasi-weeks of March, indicating a cubic change in the δ parameter.

In summary, the $model.selection output highlights the prevalence of stationary models while also capturing specific periods with notable changes in mean or dispersion of the rainfall frequency distributions. Figure 1 depicts the selected models for each quasi-week period.

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

$Statistics

While $model.selection identifies the quasi-weekly periods where the SPIChanges() function selected non-stationary models, the $Statistics output field provides detailed insights into how the parameters of the gamma-based models changed from 1827 to 2020. Specifically, $Statistics illustrates the year-to-year changes in the magnitude of these parameters (Figure 2).

Key Observations:

Quasi-Weekly Periods with Model 2:
- The four quasi-weekly periods in January (Winter season) exhibited changes toward rainier conditions, as reflected by an increase in the μ parameter.
- Conversely, the 3rd and 4th quasi-weeks of July (Summer season) demonstrated changes toward drier conditions, with a decrease in the μ parameter.
Quasi-Weekly Periods with Model 3:
- The 4th quasi-week of June, the 1st and 2nd quasi-weeks of July, and the 3rd quasi-week of September showed increasing δ parameter values. This indicates that the dispersion of these precipitation distributions expanded over time, while their mean values remained constant.
- In contrast, the 2nd quasi-week of December exhibited a decrease in the δ parameter, signifying a reduction in dispersion over time.
Quasi-Weekly Periods with Model 6:
- In the 1st quasi-week of September, the δ parameter increased from 1827 to the mid-1940s, followed by a slower decrease until 2019. Overall, the δ parameter showed an upward trend over the 193-year period.
Quasi-Weekly Periods with Model 11:
- For the 2nd and 3rd quasi-weeks of March, the δ parameter displayed a complex trajectory:
  - 1827 to ~1885: Decreased.
  - 1885 to Late 1960s: Increased.
  - Late 1960s to 2019: Decreased again.
- Over the entire series (1827–2019), the δ parameter experienced an overall decline, particularly from the middle of the series (approximately 1923) to 2019.

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
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)
#> <environment: R_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.Freq.Drought

The $Changes.Freq.Drought output field evaluates how temporal changes in the parameters of the GAMLSS models impacted the probability of drought occurrence in Oxford (Figure 3).

Key Observations:

Periods with Changes Toward Rainier Conditions:
- Quasi-weekly periods showing decreases in the expected frequency of drought events include:
  - The four quasi-weeks of January.
  - The 2nd and 3rd quasi-weeks of March.
  - The 1st quasi-week of September.
  - The 2nd and 4th quasi-weeks of December.
- These decreases in the expected frequency of moderate to extreme, severe to extreme, and extreme drought events are captured in the columns ChangeMod, ChangeSev, and ChangeExt of the output field.
- Consequently, the precipitation amounts obtained from the stationary approach for a cumulative probability of 0.5 (StatNormalRain) are lower than those derived from the nonstationary approach (NonStatNormalRain). This indicates that climatological expected precipitation for these nine quasi-week periods increased over time, reflecting a trend toward rainier conditions.
Periods with Changes Toward Drier Conditions:
- Quasi-weekly periods showing increases in the expected frequency of drought events include:
  - The 4th quasi-week of June.
  - The four quasi-weeks of July.
  - The 3rd quasi-week of September.
- These increases in the expected frequency of moderate to extreme, severe to extreme, and extreme drought events are captured in the columns ChangeMod, ChangeSev, and ChangeExt of the output field.
- As a result, the precipitation amounts obtained from the stationary approach for a cumulative probability of 0.5 (StatNormalRain) are higher than those derived from the nonstationary approach (NonStatNormalRain). This indicates that the climatological expected precipitation for these six quasi-week periods decreased over time, reflecting a trend toward drier conditions.

Summary:

The $Changes.Freq.Drought field provides valuable insights into the temporal dynamics of drought probabilities in Oxford. It highlights specific quasi-weekly periods with notable changes toward either rainier or drier conditions, emphasizing the importance of adopting nonstationary approaches for understanding change changes impacts over drought frequency.

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

$data.week

Finally, the output data field $data.week contains the precipitation amounts (rain.at.TS), SPI values, expected frequencies of the SPI estimates calculated under both the stationary approach (Exp.Acum.Prob) and the non-stationary approach (Actual.Acum.Prob), as well as the changes in the frequency of dry events (SPI < 0) caused by the changes in the precipitation patterns (ChangeDryFreq; in percentage). This output data field allowed us to assess the impact of the changes in the parameters of the gamma-based models (as demontrated by $Statistics) on the expected frequency of dry events over the entire data record (1827–2020).

As an example, let’s evaluate this output data field for the last year of the series (2019; Table 1).

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
#>      Year Month quasiWeek rain.at.TS    SPI Exp.Acum.Prob Actual.Acum.Prob
#> 9214 2019     1         1       41.5 -0.328         0.371            0.251
#> 9215 2019     1         2       41.2 -0.345         0.365            0.255
#> 9216 2019     1         3       18.1 -1.506         0.066            0.039
#> 9217 2019     1         4       28.6 -0.848         0.198            0.145
#> 9218 2019     2         1       51.0  0.141         0.556            0.556
#> 9219 2019     2         2       61.3  0.480         0.684            0.684
#> 9220 2019     2         3       52.8  0.341         0.633            0.633
#> 9221 2019     2         4       40.4  0.149         0.559            0.559
#> 9222 2019     3         1       30.3 -0.281         0.389            0.389
#> 9223 2019     3         2       37.9  0.156         0.562            0.476
#> 9224 2019     3         3       41.5  0.273         0.608            0.541
#> 9225 2019     3         4       35.5 -0.065         0.474            0.474
#> 9226 2019     4         1       36.2 -0.082         0.467            0.467
#> 9227 2019     4         2       29.2 -0.447         0.327            0.327
#> 9228 2019     4         3       24.7 -0.570         0.284            0.284
#> 9229 2019     4         4       29.9 -0.328         0.371            0.371
#> 9230 2019     5         1       28.2 -0.462         0.322            0.322
#> 9231 2019     5         2       21.2 -0.953         0.170            0.170
#> 9232 2019     5         3       21.2 -1.057         0.145            0.145
#> 9233 2019     5         4       21.3 -1.165         0.122            0.122
#> 9234 2019     6         1       37.9 -0.380         0.352            0.352
#> 9235 2019     6         2       80.3  0.818         0.793            0.793
#> 9236 2019     6         3       94.2  1.096         0.864            0.864
#> 9237 2019     6         4       99.4  1.294         0.902            0.861
#> 9238 2019     7         1       71.4  0.694         0.756            0.738
#> 9239 2019     7         2       24.4 -0.860         0.195            0.252
#> 9240 2019     7         3       26.9 -0.688         0.246            0.321
#> 9241 2019     7         4       42.9 -0.241         0.405            0.508
#> 9242 2019     8         1       47.2 -0.221         0.413            0.413
#> 9243 2019     8         2       74.5  0.554         0.710            0.710
#> 9244 2019     8         3       72.0  0.457         0.676            0.676
#> 9245 2019     8         4       46.9 -0.208         0.418            0.418
#> 9246 2019     9         1       45.0 -0.085         0.466            0.425
#> 9247 2019     9         2       18.1 -1.394         0.082            0.082
#> 9248 2019     9         3        4.9 -2.621         0.004            0.017
#> 9249 2019     9         4       77.9  0.732         0.768            0.768
#> 9250 2019    10         1       99.3  1.171         0.879            0.879
#> 9251 2019    10         2      157.1  2.156         0.984            0.984
#> 9252 2019    10         3      168.2  2.092         0.982            0.982
#> 9253 2019    10         4      120.0  1.303         0.904            0.904
#> 9254 2019    11         1      125.2  1.376         0.916            0.916
#> 9255 2019    11         2      116.2  1.264         0.897            0.897
#> 9256 2019    11         3      111.2  1.214         0.888            0.888
#> 9257 2019    11         4      105.6  1.279         0.900            0.900
#> 9258 2019    12         1       82.9  0.832         0.797            0.797
#> 9259 2019    12         2       74.4  0.654         0.743            0.755
#> 9260 2019    12         3      101.7  1.307         0.904            0.904
#> 9261 2019    12         4       87.3  0.932         0.824            0.743
#>      ChangeDryFreq
#> 9214           -12
#> 9215           -11
#> 9216          -2.7
#> 9217          -5.3
#> 9218     NoDrought
#> 9219     NoDrought
#> 9220     NoDrought
#> 9221     NoDrought
#> 9222             0
#> 9223     NoDrought
#> 9224     NoDrought
#> 9225             0
#> 9226             0
#> 9227             0
#> 9228             0
#> 9229             0
#> 9230             0
#> 9231             0
#> 9232             0
#> 9233             0
#> 9234             0
#> 9235     NoDrought
#> 9236     NoDrought
#> 9237     NoDrought
#> 9238     NoDrought
#> 9239           5.7
#> 9240           7.6
#> 9241          10.3
#> 9242             0
#> 9243     NoDrought
#> 9244     NoDrought
#> 9245             0
#> 9246          -4.1
#> 9247             0
#> 9248           1.2
#> 9249     NoDrought
#> 9250     NoDrought
#> 9251     NoDrought
#> 9252     NoDrought
#> 9253     NoDrought
#> 9254     NoDrought
#> 9255     NoDrought
#> 9256     NoDrought
#> 9257     NoDrought
#> 9258     NoDrought
#> 9259     NoDrought
#> 9260     NoDrought
#> 9261     NoDrought

For the 33 quasi-week periods where the stationary model was selected as the best fitting model (see $model.selection), the values for Exp.Acum.Prob and Actual.Acum.Prob are identical, and ChangeFreq is zero (Table 1). However, for the quasi-week periods where the SPIChanges package selected a nonstationary model (see $model.selection), Exp.Acum.Prob and Actual.Acum.Prob differed, indicating that the expected frequency of dry events in 2019 was influenced by changes in precipitation patterns. To better understand this, recall that the four quasi-weeks of January exhibited changes toward rainier conditions between 1827 and 2020 ($Changes.Freq.Drought). Therefore, the expected frequency of dry events in the later years of the series (e.g., January 2019; Table 1) was lower than in the earlier years. Specifically, the $data.week data field indicated that the expected frequency of the four dry events observed in January 2019, when estimated from the best-fitting model, was 13%, 11.5%, 2.8%, and 5.5% lower than those calculated using the stationary approach of the original SPI method. Similar results were found for the 1st quasi-week period in September (Table 1).

In contrast, July experienced changes toward drier conditions from 1827 to 2020 ($Changes.Freq.Drought). Accordingly, the $data.week data field indicated that the expected frequency of the three dry events observed in July 2019, when estimated from the best-fitting model, was 5.7%, 7.3%, and 10.6% higher than those calculated using the stationary approach of the original SPI. Similar results were seen in the 3rd quasi-week period of September (Table 1).

Interpreting the Output Data Fields from the “Oxford.Changes.linear”

The interpretation of the output data fields $model.selection, $Statistics, $Changes.Freq.Drought, and $data.week from the “Oxford.Changes.linear” data frame is similar to those from the “Oxford.Changes” data frame, described in Figure 1. However, since “Oxford.Changes.linear” was calculated with the argument only.linear of the SPIChanges() function set to "Yes", only 4 models (a stationary and three linear gamma-based models) were considered to describe changes in the rainfall patterns in Oxford (Figure 4), UK. As with the output data fields of the “Oxford.Changes” data frame, the outputs of “Oxford.Changes.linear” also indicate that most of the quasi-weekly periods did not exhibit significant changes in the frequency of meteorological droughts, see below.

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

$model.selection

This output data field indicates that 34 out of 48 quasi-weekly periods did not exhibit significant changes in the frequency of meteorological droughts, as the stationary gamma distribution (Model 1) was the best-fitting model for these precipitation series.

This data field also reveals the following:

Model 2, which describes a linear change only in the μ parameter of the gamma-based models, was selected for:
- The four quasi-weeks of January
- The 3rd and 4th quasi-weeks of July
Model 3, which shows a linear change only in the δ parameter of the gamma distribution, was selected for:
- The 3rd quasi-week of March
- The 4th quasi-week of June
- The 1st and 2nd quasi-weeks of July
- The 1st and 3rd quasi-weeks of September
- The 2nd quasi-week of December
Model 4, which describes a linear change in both the μ and δ parameters of the gamma-based models, was selected for:
- The 4th quasi-week of December

These findings indicate that while many periods were stationary, significant linear changes in rainfall patterns were observed in specific quasi-weekly periods.

$Statistics and $Changes.Freq.Drought

As one might have anticipated, the only quasi-weeks where these output data fields show different results from those of the “Oxford.Changes” data frame object are those where the SPIChanges() function with only.linear = No selected nonstationary models. In other words, different nonstationary models were selected for the 3rd quasi-week of March and the 1st quasi-week of September.

When only.linear = No, models 11 and 6 were chosen for these periods, respectively.
When only.linear = Yes, model 3 was selected for both periods.

It is important to note that these three models (11, 6, and 3) all assume that only the δ parameter changes over time. However, setting only.linear = No provided a more complete (non-linear) description of how the dispersion of the precipitation frequency distributions changed from 1827 to 2019.

The major difference between using only.linear = No and only.linear = Yes is observed in the 2nd quasi-week of March. With only.linear = No, the SPIChanges() function selected model 11, indicating changes toward drier conditions. In contrast, with only.linear = Yes, model 1 was selected, suggesting no changes in precipitation patterns.

This discrepancy highlights that when only stationary and linear models are considered, the SPIChanges() function may fail to capture and describe complex changes in the dispersion parameter of the gamma-based model.

$data.week

As described above, the output data field, $data.week contains the precipitation amounts (rain.at.TS), SPI values, expected frequencies of the SPI estimates calculated under both the stationary approach (Exp.Acum.Prob) and the non-stationary approach (Actual.Acum.Prob), as well as the changes in the frequency of dry events (SPI < 0) caused by the changes in the precipitation patterns (ChangeDryFreq; in percentage). This output data field allowed us to assess the impact of the linear changes in the parameters of the gamma-based models (as demonstrated by $Statistics) on the expected frequency of dry events over the entire data record (1827–2020).

Once again, let’s evaluate this output data field for the last year of the series (2019; Table 2).

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
#>      Year Month quasiWeek rain.at.TS    SPI Exp.Acum.Prob Actual.Acum.Prob
#> 9214 2019     1         1       41.5 -0.328         0.371            0.251
#> 9215 2019     1         2       41.2 -0.345         0.365            0.255
#> 9216 2019     1         3       18.1 -1.506         0.066            0.039
#> 9217 2019     1         4       28.6 -0.848         0.198            0.145
#> 9218 2019     2         1       51.0  0.141         0.556            0.556
#> 9219 2019     2         2       61.3  0.480         0.684            0.684
#> 9220 2019     2         3       52.8  0.341         0.633            0.633
#> 9221 2019     2         4       40.4  0.149         0.559            0.559
#> 9222 2019     3         1       30.3 -0.281         0.389            0.389
#> 9223 2019     3         2       37.9  0.156         0.562            0.562
#> 9224 2019     3         3       41.5  0.273         0.608            0.568
#> 9225 2019     3         4       35.5 -0.065         0.474            0.474
#> 9226 2019     4         1       36.2 -0.082         0.467            0.467
#> 9227 2019     4         2       29.2 -0.447         0.327            0.327
#> 9228 2019     4         3       24.7 -0.570         0.284            0.284
#> 9229 2019     4         4       29.9 -0.328         0.371            0.371
#> 9230 2019     5         1       28.2 -0.462         0.322            0.322
#> 9231 2019     5         2       21.2 -0.953         0.170            0.170
#> 9232 2019     5         3       21.2 -1.057         0.145            0.145
#> 9233 2019     5         4       21.3 -1.165         0.122            0.122
#> 9234 2019     6         1       37.9 -0.380         0.352            0.352
#> 9235 2019     6         2       80.3  0.818         0.793            0.793
#> 9236 2019     6         3       94.2  1.096         0.864            0.864
#> 9237 2019     6         4       99.4  1.294         0.902            0.861
#> 9238 2019     7         1       71.4  0.694         0.756            0.738
#> 9239 2019     7         2       24.4 -0.860         0.195            0.252
#> 9240 2019     7         3       26.9 -0.688         0.246            0.321
#> 9241 2019     7         4       42.9 -0.241         0.405            0.508
#> 9242 2019     8         1       47.2 -0.221         0.413            0.413
#> 9243 2019     8         2       74.5  0.554         0.710            0.710
#> 9244 2019     8         3       72.0  0.457         0.676            0.676
#> 9245 2019     8         4       46.9 -0.208         0.418            0.418
#> 9246 2019     9         1       45.0 -0.085         0.466            0.511
#> 9247 2019     9         2       18.1 -1.394         0.082            0.082
#> 9248 2019     9         3        4.9 -2.621         0.004            0.017
#> 9249 2019     9         4       77.9  0.732         0.768            0.768
#> 9250 2019    10         1       99.3  1.171         0.879            0.879
#> 9251 2019    10         2      157.1  2.156         0.984            0.984
#> 9252 2019    10         3      168.2  2.092         0.982            0.982
#> 9253 2019    10         4      120.0  1.303         0.904            0.904
#> 9254 2019    11         1      125.2  1.376         0.916            0.916
#> 9255 2019    11         2      116.2  1.264         0.897            0.897
#> 9256 2019    11         3      111.2  1.214         0.888            0.888
#> 9257 2019    11         4      105.6  1.279         0.900            0.900
#> 9258 2019    12         1       82.9  0.832         0.797            0.797
#> 9259 2019    12         2       74.4  0.654         0.743            0.755
#> 9260 2019    12         3      101.7  1.307         0.904            0.904
#> 9261 2019    12         4       87.3  0.932         0.824            0.743
#>      ChangeDryFreq
#> 9214           -12
#> 9215           -11
#> 9216          -2.7
#> 9217          -5.3
#> 9218     NoDrought
#> 9219     NoDrought
#> 9220     NoDrought
#> 9221     NoDrought
#> 9222             0
#> 9223     NoDrought
#> 9224     NoDrought
#> 9225             0
#> 9226             0
#> 9227             0
#> 9228             0
#> 9229             0
#> 9230             0
#> 9231             0
#> 9232             0
#> 9233             0
#> 9234             0
#> 9235     NoDrought
#> 9236     NoDrought
#> 9237     NoDrought
#> 9238     NoDrought
#> 9239           5.7
#> 9240           7.6
#> 9241          10.3
#> 9242             0
#> 9243     NoDrought
#> 9244     NoDrought
#> 9245             0
#> 9246           4.5
#> 9247             0
#> 9248           1.2
#> 9249     NoDrought
#> 9250     NoDrought
#> 9251     NoDrought
#> 9252     NoDrought
#> 9253     NoDrought
#> 9254     NoDrought
#> 9255     NoDrought
#> 9256     NoDrought
#> 9257     NoDrought
#> 9258     NoDrought
#> 9259     NoDrought
#> 9260     NoDrought
#> 9261     NoDrought

For the 34 quasi-week periods where the stationary model was selected as the best fitting model (see $model.selection), the values for Exp.Acum.Prob and Actual.Acum.Prob are identical, and ChangeFreq is zero (Table 2). However, for the quasi-week periods where the SPIChanges package selected linear nonstationary models (see $model.selection), Exp.Acum.Prob and Actual.Acum.Prob differed, indicating that the expected frequency of dry events in 2019 was influenced by changes in precipitation patterns.

To better understand this, recall that the four quasi-weeks of January exhibited changes toward rainier conditions between 1827 and 2020 ($Changes.Freq.Drought). Therefore, the expected frequency of dry events in the later years of the series (e.g., January 2019; Table 2) was lower than in the earlier years. Specifically, the $data.week data field indicated that the expected frequency of the four dry events observed in January 2019, when estimated from the best-fitting model, was 13%, 11.5%, 2.8%, and 5.5% lower than those calculated using the stationary approach of the original SPI method. Similar results were found for the 1st quasi-week period in September (Table 2).

In contrast, the 2nd, 3rd, and 4th quasi-weeks of July experienced changes toward drier conditions from 1827 to 2020 ($Changes.Freq.Drought). Accordingly, the $data.week data field indicated that the expected frequency of the three dry events observed in July 2019, when estimated from the best-fitting model, was 5.7%, 7.3%, and 10.7% higher than those calculated using the stationary approach of the original SPI. Similar results are observed in the 3rd quasi-week period of September (Table 2).

Case Study 2 - Brazil Drought Events

In this case study we apply the SPIChanges package for entire Brazil (1980-2024), downloading data from CPC Global Unified Gauge-Based Analysis of Daily Precipitation data provided by the NOAA PSL, Boulder, Colorado, USA, from their website at https://psl.noaa.gov.

# 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)

Aggregating daily data into a specified time scale:

For this analysis, we set the ‘TS’ argument in the TSaggreg() function to 12, which corresponds to a moving window of three months. Shorter time scales may result in precipitation frequency distributions where the percentage of zeros exceeds 80%, particularly in Brazil’s semi-arid regions.

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]
}

Applying the SPIChanges::SPIChanges() for entire Brazil

As previously descried, we applied the SPIChanges package to the entire country of Brazil. At the 12-quasi-week time scale, the 4th quasi-week of February, May, August, and November correspond to precipitation data accumulated over the Austral Summer, Autumn, Winter, and Spring seasons, respectively. Regarding changes in severe to extreme drought events (SPI ≤ -1.5), the results of the chunk below are consistent with previous studies indicating that Brazil, like other regions in South America, may become a drought hotspot due to its potential to respond drastically to excessive drying and warming.

Note that in this section, we use parallel processing to speed up the calculations using {foreach}, if you do not wish to use multiple cores, change numCores <- detectCores() - 1 to just read numCores <- 1.

# 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)

Assessing changes in severe to moderate drought events

This next chunk uses the file Map (from the previous chunk) to plot the changes in the expected frequency of severe to moderate drought events (SPI ≤ −1.5), calculated at 12- quasi-week (3-month) time scale in Brazil. The analysis of this map indicated that 50.8 (Summer), 45.5 (Autumn), 35.6 (Winter), and 49.2% (Spring) of Brazil experienced changes toward drier conditions between 1980 and 2024 as indicated by the red colours of the map. This change toward more severe to extreme drought events in almost half of Brazil is in line with the widespread occurrence of unprecedented drought events in the country after 2000.

library(ggplot2)
library(sf)
#> Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
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)
#> Reading layer `regioes_2010' from data source 
#>   `/private/var/folders/2y/gpykxwfs075fhp74nhn1skk00000gn/T/Rtmpc6qFoq/regioes_2010.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 5 features and 5 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -73.99024 ymin: -33.75136 xmax: -32.39088 ymax: 5.270972
#> Geodetic CRS:  WGS 84
# 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)
#> Warning: st_centroid assumes attributes are constant over geometries
# 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)

References

Burt, Stephen. 2020. Daily Meteorological Data from the Radcliffe Observatory (Now Radcliffe Meteorological Station), from January 1815: Updated to January 2020. Figshare. https://doi.org/10.6084/M9.FIGSHARE.11956239.V1.

Wu, Hong, Mark D. Svoboda, Michael J. Hayes, Donald A. Wilhite, and Fujiang Wen. 2006. “Appropriate Application of the Standardized Precipitation Index in Arid Locations and Dry Seasons.” International Journal of Climatology 27 (1): 65–79. https://doi.org/10.1002/joc.1371.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.