The Agronomic Institute and the Weather Station of Campinas-SP, Brazil

Founded in 1886 by Emperor Dom Pedro II, the Agronomic Institute of Campinas (IAC), located in the state of Sao Paulo, is one of Brazil’s pioneering research centers.
Since 1890, the institution has maintained daily records of air temperature and precipitation, creating one of the longest meteorological series in South America (Blain et al. 2009 and Blain et al. 2018).
We created this R-Markdown file to facilitate the reproducibility of the results found in the study A CENTURY AND BEYOND: TRENDS AND IMPLICATIONS OF CLIMATE CHANGE FROM AGRONOMIC INSTITUTE’S HISTORICAL WEATHER DATASET by Prela-Pantano et al. (2025).
The goal of this latter study was to evaluate this long-term meteorological series to detect signs of climate change and place recent weather events within a historical context beginning in the late 19th century.
The hypothesis of the study was that these 135 years of continuous meteorological records provide unique and valuable insights into the climate evolution of this tropical region.
The goals of the study were (i) to detect and describe signs of climate change in this historical data records and, (ii) to make the historical records readily and public available.
The study also applied distinct data quality assessment methods to ensure the reliability of the meteorological series, before making it available in this document.

Please click on the topics of the menu (above) the see the data and analysis. Enjoy!

Authors

Angelica Prela-Pantano⁽¹⁾ ORCID
Orivaldo Brunin⁽²⁾ ORCID
Graciela da Rocha Sobierajski⁽¹⁾ ORCID
Ludmila Bardin-Camparotto⁽¹⁾ ORCID
Letícia Lopes Martins⁽¹⁾ ORCID
Gabriel Constantino Blain⁽¹⁾* ORCID
*Autor correspondente: gabriel.blain@sp.gov.br

1 Instituto Agronômico de Campinas (IAC/APTA/SAA-SP)

2 Fundação de Apoio A Pesquisa Agrícola

Required R-packages

Installing and loading the packages required to run this R-Markdown file.

  if (!require(gamlss, quietly = TRUE)) install.packages("gamlss")
  if (!require(gamlss.dist, quietly = TRUE)) install.packages("gamlss.dist")
  if (!require(MuMIn, quietly = TRUE)) install.packages("MuMIn")
  if (!require(spsUtil, quietly = TRUE)) install.packages("spsUtil")
  if (!require(trend, quietly = TRUE)) install.packages("trend")
  if (!require(pracma, quietly = TRUE)) install.packages("pracma")
  if (!require(lubridate, quietly = TRUE)) install.packages("lubridate")
  if (!require(dplyr, quietly = TRUE)) install.packages("dplyr")
  if (!require(KScorrect, quietly = TRUE)) install.packages("KScorrect")
  if (!require(ggplot2, quietly = TRUE)) install.packages("ggplot2")
  if (!require(gt, quietly = TRUE)) install.packages("gt")
  if (!require(reactable, quietly = TRUE)) install.packages("reactable")
  if (!require(htmltools, quietly = TRUE)) install.packages("htmltools")
  if (!require(SPEI, quietly = TRUE)) install.packages("SPEI")
  library(gamlss)
  library(gamlss.dist)
  library(MuMIn)
  library(spsUtil)
  library(trend)
  library(pracma)
  library(lubridate)
  library(dplyr)
  library(KScorrect)
  library(ggplot2)
  library(gt)
  library(reactable)
  library(htmltools)
  library(SPEI)

Custom functions

We created the R-functions Changepoint(), test.Good() and select.model() to help implementing the methods proposed in this study.

The Changepoint() function was designed to find abrupt changes in the meteorological series. It starts by evaluating the presence of trends in the data using the Mann-Kendall test. If the trend is found significant, the function performs a locally weighted regression (loess) to detrend the series. Then, the Pettitt change point test is applied to detect abrupt changes that may be caused by changes in instruments or in the weather station sites.

Changepoint <- function (dataseries){
  result <- as.data.frame(matrix(NA,1,3))
  result[1,3] <- "NoChange"
  n <- length(dataseries)
  dataseries.med <- mean(dataseries)
  dataseries.sd <- sd(dataseries)
  dataseries <- (dataseries-dataseries.med)/dataseries.sd
  mk <- mk.test(dataseries)
  mk$p.value <- round(mk$p.value,3)
  if (mk$p.value < 0.05) {
    time <- 1L:n
    loess_fit <- loess(dataseries ~ time, span = 0.5)
    dataseries <- dataseries - predict(loess_fit)}
    Pettitt <- pettitt.test(dataseries)
    Pettitt$p.value <- round(Pettitt$p.value,3)
    if (mk$p.value < 0.001) {mk$p.value <- c("< 0.001")}
    if (Pettitt$p.value == 1) {Pettitt$p.value <- c("> 0.999")}
    result[1,1:2] <- c(round(mk$statistic,2),mk$p.value)
    result[1,3] <- Pettitt$p.value
    colnames(result) <- c("MK_Stat","MK_Pvalue","Pet_Pvalue")
    return(result)
}

The test.Good() function applies the Lilliefors test and the second-order Akaike Information Criterion to assess how well the normal, log-normal, gamma, and Weibull distributions fit the data. It is important to emphasize that this function was designed specifically for tropical and subtropical conditions, where annual and monthly average air temperature values are consistently greater than zero. Further information on the Lilliefors test can be found in Blain, 2014.

test.Good <- function(series){
  series.nozero <- series[series > 0]
  Lilli <- matrix(NA,1,8)
  Normal <-  quiet(LcKS(series.nozero, cdf="pnorm", nreps=5000))
  LogNormal <-  quiet(LcKS(series.nozero, cdf="plnorm", nreps=5000))
  Gamma <-  quiet(LcKS(series.nozero, cdf="pgamma", nreps=5000, parallel = TRUE))
  Weibull <-  quiet(LcKS(series.nozero, cdf="pweibull", nreps=5000, parallel = TRUE))
  Lilli[1,1:4] <- c(Normal$p.value, LogNormal$p.value, Gamma$p.value, Weibull$p.value)
  t.nor <- gamlss::gamlss(series.nozero ~ 1, family = NO, mu.link = "log", sigma.link = "log")
  t.lnor <- gamlss::gamlss(series.nozero ~ 1, family = LNO, mu.link = "log", sigma.link = "log")
  t.gam <- gamlss::gamlss(series.nozero ~ 1, family = GA, mu.link = "log", sigma.link = "log")
  t.wei <- gamlss::gamlss(series.nozero ~ 1, family = WEI, mu.link = "log", sigma.link = "log")
  Akaik <- MuMIn::AICc(t.nor,t.lnor,t.gam,t.wei, k = 4.4)
  Lilli[1,5:8] <- Akaik$AICc
  colnames(Lilli) <- c("LilliNor","LilliLogNor","LilliGam","LilliWei",
                       "AICcNor","AICcLogNor","AICcGam","AICcWei")
  Lilli <- round(Lilli,3)
  return(Lilli)
}

We developed the select.model() function under the framework of the Generalized Linear Additive models combined with Generalized Additive Models for distribution’s parameters (GAMLSS). The GAMLSS method allows great flexibility in modeling the frequency distribution of response variables and are able to incorporate both linear and nonlinear functions Rigby and Stasinopoulos, 2005. The GAMLSS models considered in this study are described as follow.

Model 1 (stationary): The mean (µ) and dispersion (δ) of the distribution do not change over time.
Model 2: Only µ varies linearly over time.
Model 3: Only δ varies linearly over time.
Model 4: Both µ and δ vary linearly over time.
Model 5: Only µ varies non-linearly over time as a quadratic polynomial function.
Model 6: µ varies non-linearly over time as a quadratic polynomial function and δ vary linearly over time.

The AICc was again used to select the best-fitting GAMLSS model for each series. For the cases where the AICc selected model 1 (the stationary model) we accepted the hypothesis that data’s frequency distribution did not change over time. In other words, we regarded the series as free from climate change signs. On the other hand, we regarded the selection of a nonstationary model (models 2 to 6) as evidence that data’s frequency distribution significantly changes over time. In other words, we regarded the series as having significant climate change signs. As described in Blain et al. (2022), this parametric trend analysis is not only able to detect changes in data’s frequency distribution; it is also capable of isolating its effect on the central tendency and/or dispersion of the meteorological series. This potentially leads to a broad understanding of effect of the current climate change on the frequency and intensity of weather events.

select.model <- function(series, dist) {
  if (dist != "GA" & dist != "NO" & dist != "LNO" & dist != "WEI") {
    stop("argument dist must be a chacter variable equal to NO, LNO, GA or WEI")
  }
  series <- series[series>0]
  n <- length(series)
  time <- 1L:n
  model <- quiet(gamlss::gamlss(series ~ 1, family = dist,
                          mu.link = "log", sigma.link = "log"))
  model.ns10 <- quiet(gamlss::gamlss(series ~ time, family = dist,
                               mu.link = "log", sigma.link = "log"))
  model.ns01 <- quiet(gamlss::gamlss(series ~ 1, sigma.formula = ~ time,
                               family = dist, mu.link = "log", sigma.link = "log"))
  model.ns11 <- quiet(gamlss::gamlss(series ~ time, sigma.formula = ~ time,
                               family = dist, mu.link = "log", sigma.link = "log"))
  model.ns20 <- quiet(gamlss::gamlss(series ~ time + I(time^2), family = dist, mu.link = "log",
                               sigma.link = "log"))
  model.ns21 <- quiet(gamlss::gamlss(series ~ time + I(time^2), sigma.formula = ~ time, family = dist,
                               mu.link = "log", sigma.link = "log"))
  
  #Selection among the candidate models  
  best <- which.min(list(
    MuMIn::AICc(model, k = 4.4),
    MuMIn::AICc(model.ns10, k =  4.4),
    MuMIn::AICc(model.ns01, k =  4.4),
    MuMIn::AICc(model.ns11, k =  4.4),
    MuMIn::AICc(model.ns20, k =  4.4),
    MuMIn::AICc(model.ns21, k =  4.4)
  ))
  return(best)
}

Loading the historical data

This chunk loads the historical daily records and calculates the potential evapotranspiration (PE) and the difference between precipitation and EP. The daily EP is estimated by the Hargreaves & Samani (1985) equation from the CropWaterBalane package (Blain et al. 2024). We used this EP estimation method because of the lack of necessary variables for calculating this agro-meteorological parameter by methods such as the FAO Penman-Monteith (Allen et al., 1998).

# Load data
url <- "https://raw.githubusercontent.com/gabrielblain/IAC_Hist_Weather/refs/heads/main/IAC_Hist_Weather.csv"
rawdata <- read.csv(url,sep=",", header = TRUE)
tmax <- rawdata$tmax
rainfall <- rawdata$rain
tmin <- rawdata$tmin
tavg <- (tmin+tmax)/2
# Create dates and years
dates <- seq.Date(from = as.Date("1890-01-01"), to = as.Date("2024-12-31"), by = "days")
n <- length(dates)
years <- as.numeric(format(dates, "%Y"))
month <- as.numeric(format(dates, "%m"))
days <- as.numeric(format(dates, "%d"))
# Create data frame
data <- data.frame(
  dates = as.Date(dates),
  years = years,
  month = month,
  days = days,
  rain = rainfall,
  tmin = tmin,
  tmax = tmax,
  tavg = round(tavg, 2)
)
colnames(data) <- c("dates", "years", "month", "days", "rain", "tmin", "tmax", "tavg")
initial.year <- year(dates[1])
final.year <- year(dates[n])
All.years <- seq(from = initial.year, to=final.year)
base64_csv <- base64enc::base64encode(url)
# Render data table with Reactable
reactable::reactable(
  data,
  defaultColDef = reactable::colDef(
    style = list(textAlign = "center"),        # Align cell contents
    headerStyle = list(textAlign = "center")  # Align header text
  ),
  highlight = TRUE,  # Add highlight on row hover
  bordered = TRUE,   # Add borders around cells
  striped = TRUE     # Add alternating row colors for readability
)

This data series can be downloaded under the tag Final Remarks and Data Download.

Plots: Daily values

Daily precipitation records

# Convert 'dates' to Date format
data$dates <- as.Date(data$dates)
start_date <- as.Date("1890-01-01")
end_date <- as.Date("2024-12-31")
# Create the daily precipitation plot
DailyRain <- ggplot(data, aes(x = dates, y = rain)) +
  geom_point(color = "blue", size = 1) +      # Optional: points for daily values
  labs(
    title = "Daily Precipitation - Campinas, SP",
    x = "Date",
    y = "Precipitation (mm)"
  ) +
  scale_x_date(
    limits = c(start_date, end_date),
    breaks = seq(start_date, end_date, by = "4 years"),
    date_labels = "%Y"                       # Format as year (e.g., 1890, 1900)
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1), # Rotate x-axis labels
    plot.title = element_text(hjust = 0.5)            # Center the title
  )

# Display the plot
DailyRain

Daily minimum air temperature records

# Convert 'dates' to Date format
data$dates <- as.Date(data$dates)
start_date <- as.Date("1890-01-01")
end_date <- as.Date("2024-12-31")
# Create the daily Tmax plot
DailyTmin <- ggplot(data, aes(x = dates, y = tmin)) +
  geom_point(color = "blue", size = 1) +      # Optional: points for daily values
  labs(
    title = "Daily Minimum Air temperature - Campinas, SP",
    x = "Date",
    y = "Tmin (\u00BAC)"
  ) +
  scale_x_date(
    limits = c(start_date, end_date),
    breaks = seq(start_date, end_date, by = "4 years"),                # Breaks every 4 years
    date_labels = "%Y"                       # Format as year (e.g., 1890, 1900)
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1), # Rotate x-axis labels
    plot.title = element_text(hjust = 0.5)            # Center the title
  )

# Display the plot
DailyTmin

Daily maximm air temperature records

# Convert 'dates' to Date format
data$dates <- as.Date(data$dates)
start_date <- as.Date("1890-01-01")
end_date <- as.Date("2024-12-31")
# Create the daily Tmax plot
DailyTmax <- ggplot(data, aes(x = dates, y = tmax)) +
  geom_point(color = "red", size = 1) +      # Optional: points for daily values
  labs(
    title = "Daily Maximum Air temperature - Campinas, SP",
    x = "Date",
    y = "Tmax (\u00BAC)"
  ) +
  scale_x_date(
    limits = c(start_date, end_date),
    breaks = seq(start_date, end_date, by = "4 years"), 
    date_labels = "%Y"                       # Format as year (e.g., 1890, 1900)
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1), # Rotate x-axis labels
    plot.title = element_text(hjust = 0.5)            # Center the title
  )

# Display the plot
DailyTmax

Daily average air temperature records

# Convert 'dates' to Date format
data$dates <- as.Date(data$dates)
start_date <- as.Date("1890-01-01")
end_date <- as.Date("2024-12-31")
# Create the daily Tmax plot
DailyTavg <- ggplot(data, aes(x = dates, y = tmax)) +
  geom_point(color = "black", size = 1) +      # Optional: points for daily values
  labs(
    title = "Daily Average Air temperature - Campinas, SP",
    x = "Date",
    y = "tavg (\u00BAC)"
  ) +
  scale_x_date(
    limits = c(start_date, end_date),
    breaks = seq(start_date, end_date, by = "4 years"), 
    date_labels = "%Y"                       # Format as year (e.g., 1890, 1900)
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1), # Rotate x-axis labels
    plot.title = element_text(hjust = 0.5)            # Center the title
  )

# Display the plot
DailyTavg

Adjusting time scales:annual

Transforming daily to Yearly

This code aggregates daily records into annual data

Annual <- data %>%
  group_by(years) %>%
  summarise(
    Rainfall = round(sum(rain, na.rm = TRUE),1),
    Tmin = round(mean(tmin, na.rm = TRUE),1),
    Tmax = round(mean(tmax, na.rm = TRUE),1),
    Tavg = round(mean(tavg, na.rm = TRUE),1)
    )
reactable::reactable(Annual,
                     style = list(textAlign = "center"))

Plots: Annual values

Annual precipitation records

# Create the Annual precipitation plot
AnnualRain <- ggplot(Annual, aes(x = years, y = Rainfall)) +
  geom_point(color = "blue", size = 1.5) +      # Points for annual values
  geom_line(color = "blue", linetype = "solid", linewidth = 0.8) +  # Line to connect points
  labs(
    title = "Annual Precipitation - Campinas, SP",
    x = "Year",
    y = "Precipitation (mm)"
  ) +
  scale_x_continuous(
    breaks = seq(1890, 2024, by = 4),  # Breaks every 4 years
    expand = c(0, 0)                   # Remove extra space on axis
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1), # Optional: Rotate x-axis labels
    plot.title = element_text(hjust = 0.5),           # Center the title
    text = element_text(size = 12)                    # Adjust font size
  )

# Display the plot
AnnualRain

Annual minimum air temperature records

# Create the Annual precipitation plot
AnnualTmin <- ggplot(Annual, aes(x = years, y = Tmin)) +
  geom_point(color = "blue", size = 1.5) +      # Points for annual values
  geom_line(color = "blue", linetype = "solid", linewidth = 0.8) +  # Line to connect points
  labs(
    title = "Annual Means: Minimum air temperature - Campinas, SP",
    x = "Year",
    y = "Tmin (\u00BAC)"
  ) +
  scale_x_continuous(
    breaks = seq(1890, 2024, by = 4),  # Breaks every 4 years
    expand = c(0, 0)                   # Remove extra space on axis
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1), # Optional: Rotate x-axis labels
    plot.title = element_text(hjust = 0.5),           # Center the title
    text = element_text(size = 12)                    # Adjust font size
  )

# Display the plot
AnnualTmin

Annual maximum air temperature records

# Create the Annual precipitation plot
AnnualTmax <- ggplot(Annual, aes(x = years, y = Tmax)) +
  geom_point(color = "red", size = 1.5) +      # Points for annual values
  geom_line(color = "red", linetype = "solid", linewidth = 0.8) +  # Line to connect points
  labs(
    title = "Annual Means: Maximum air temperature - Campinas, SP",
    x = "Year",
    y = "Tmax (\u00BAC)"
  ) +
  scale_x_continuous(
    breaks = seq(1890, 2024, by = 4),  # Breaks every 4 years
    expand = c(0, 0)                   # Remove extra space on axis
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1), # Optional: Rotate x-axis labels
    plot.title = element_text(hjust = 0.5),           # Center the title
    text = element_text(size = 12)                    # Adjust font size
  )

# Display the plot
AnnualTmax

Annual average air temperature records

# Create the Annual precipitation plot
AnnualTavg <- ggplot(Annual, aes(x = years, y = Tavg)) +
  geom_point(color = "black", size = 1.5) +      # Points for annual values
  geom_line(color = "black", linetype = "solid", linewidth = 0.8) +  # Line to connect points
  labs(
    title = "Annual Means: Average air temperature - Campinas, SP",
    x = "Year",
    y = "Tavg (\u00BAC)"
  ) +
  scale_x_continuous(
    breaks = seq(1890, 2024, by = 4),  # Breaks every 4 years
    expand = c(0, 0)                   # Remove extra space on axis
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1), # Optional: Rotate x-axis labels
    plot.title = element_text(hjust = 0.5),           # Center the title
    text = element_text(size = 12)                    # Adjust font size
  )

# Display the plot
AnnualTavg

Adjusting time scales:Monthly

Transforming daily to Monthly

This code aggregates daily records into monthly data

data$month <- floor_date(data$dates, "month")
Monthly <- data %>%
  group_by(month) %>%
  summarise(
    Rainfall = round(sum(rain, na.rm = TRUE),1),
    Tmin = round(mean(tmin, na.rm = TRUE),1),
    Tmax = round(mean(tmax, na.rm = TRUE),1),
    Tavg = round(mean(tavg, na.rm = TRUE),1)
  )
Month <- as.numeric(format(Monthly$month, "%m"))
PE <- round(thornthwaite(Monthly$Tavg, -22.9),1)

[1] "Checking for missing values (`NA`): all the data must be complete. Input type is vector. Assuming the data are monthly time series starting in January, all regular (non-leap) years."

PPE <- round((Monthly$Rainfall-PE),1)
Monthly <- cbind(Monthly,PE,PPE,Month)
reactable::reactable(Monthly,
                     style = list(textAlign = "center"))

Evaluating abrupt changes

The Pettitt test (Pettitt, 1979) is a non-parametric method for detecting a single change-point in a time series. It tests the null hypothesis of no change-point versus the alternative hypothesis that a change-point exists. The test statistic is based on the rank sums of the observations before and after the potential change-point. Among the existing tests, the Pettitt test is one of the most widely used methods in agro-hydro-meteorological studies (e.g. Wijngaard et al. 2003, Liu et al. 2012, Mallakpour & Gabriele Villarini 2016). As any other change-point test, the performance of the Pettitt test for detecting non-climatic/artificial changes in data records may be affected by true climate change, such as a point between two phases of long quasi-periodic variation (Wang 2008) or trends imposed by the current climate change. Therefore, before applying the Pettitt test to the precipitation and air temperature series (aggregated at the annual and monthly scales) we applied the widely used Mann-Kendall test (MK; Kendall and Stuart, 1967). Whenever a significant trend was detected, we used the locally estimated scatterplot smoothing (LOESS; Cleveland, 1979 and Cleveland and Devlin 1988) to remove (decompose) long-term trends of the series and then applied the Pettitt test to the detrended series.

Change-point detection (annual scale)

Annual.Quality <- as.data.frame(matrix(NA,3,3))
Annual.Quality[1,] <- Changepoint(dataseries = Annual$Rainfall)
Annual.Quality[2,] <- Changepoint(dataseries = Annual$Tmin)
Annual.Quality[3,] <- Changepoint(dataseries = Annual$Tmax)
colnames(Annual.Quality) <- c("MK_Stat","MK_Pvalue","Pet_Pvalue")
rownames(Annual.Quality) <- c("Rain","Tmin","Tmax")
reactable::reactable(Annual.Quality,
                     style = list(textAlign = "center"))

Change-point detection (Monthly scale)

Rain.Monthly.Quality <- as.data.frame(matrix(NA,12,3))
Tmin.Monthly.Quality <- as.data.frame(matrix(NA,12,3))
Tmax.Monthly.Quality <- as.data.frame(matrix(NA,12,3))
for (i in 1:12){
  x <- Monthly[which(Monthly$Month == i),2:4]
  Rain.Monthly.Quality[i,] <- Changepoint(dataseries = x[,1])
  Tmin.Monthly.Quality[i,] <- Changepoint(dataseries = x[,2])
  Tmax.Monthly.Quality[i,] <- Changepoint(dataseries = x[,3])
}

rownames(Rain.Monthly.Quality) <- c("Rain_Jan","Rain_Feb","Rain_Mar","Rain_Apr","Rain_May","Rain_Jun",
                                   "Rain_Jul","Rain_Ago","Rain_Sep","Rain_Oct","Rain_Nov","Rain_Dec")
rownames(Tmin.Monthly.Quality) <- c("Tmin_Jan","Tmin_Feb","Tmin_Mar","Tmin_Apr","Tmin_May","Tmin_Jun",
                                   "Tmin_Jul","Tmin_Ago","Tmin_Sep","Tmin_Oct","Tmin_Nov","Tmin_Dec")
rownames(Tmax.Monthly.Quality) <- c("Tmax_Jan","Tmax_Feb","Tmax_Mar","Tmax_Apr","Tmax_May","Tmax_Jun",
                                   "Tmax_Jul","Tmax_Ago","Tmax_Sep","Tmax_Oct","Tmax_Nov","Tmax_Dec")
Monthly.Quality <- rbind(Rain.Monthly.Quality,Tmin.Monthly.Quality,Tmax.Monthly.Quality)
colnames(Monthly.Quality) <- c("MK_Stat","MK_Pvalue","Pet_Pvalue")
reactable::reactable(Monthly.Quality,style = list(textAlign = "center"))

Reference Values (1961-1900) and climate change indices

Efforts such as those of the World Meteorological Organization Commission for Climatology/Climate Variability (WMO- CCL/CLIVAR) recommends the use of feasible indices to evaluate regional signs of climate change (e.g. Karl et al., 1999 and Peterson et al., 2001). Table 1 describes the indices adopted in this study estimated at annual and monthly scales. The 10th, 90th, 95th and, 99th percentiles described in Table 1 were calculated considering the 1961-1990 standard period. Further information on these indices can be found at https://etccdi.pacificclimate.org/list_27_indices.shtml. Considering the lack of knowledge regarding the distinct frequency distributions of each of these indices (Table 1), we applied the MK test – instead of a parametric test – to detected trends in their time series.

Table 1. Climate change Indices

Index	Description	Units
	Annual Scale
warm.nights	Number of nights where Tmin > 90th percentile	Days
cool.nights	Number of nights where Tmin < 10th percentile	Days
warm.days	Number of days where Tmax > 90th percentile	Days
cool.days	Number of days where Tmax < 10th percentile	Days
SU	Number of days where Tmax > 25°C	Days
TR	Number of days where Tmin > 20°C	Days
PRCPTOT	Total number of wet days (Rainfall > 1 mm)	Days
SDII	Simple Daily Intensity Index (mean rainfall per wet day)	mm/day
R10	Average rainfall on days with Rainfall > 10 mm	mm
R20	Average rainfall on days with Rainfall > 20 mm	mm
R95	Number of days with Rainfall > 95th percentile	Days
R99	Number of days with Rainfall > 99th percentile	Days
CDD	Maximum number of consecutive dry days (Rainfall < 1 mm)	Days
CWD	Maximum number of consecutive wet days (Rainfall ≥ 1 mm)	Days
	Monthly Scale
TXx	Largest daily maximum temperature for each month	°C
TNx	Largest daily minimum temperature for each month	°C
TXn	Lowest daily maximum temperature for each month	°C
TNn	Lowest daily minimum temperature for each month	°C

Annual

ref.period <- data[which(data$years>=1961 & data$years <= 1990),]
ref.prec.mean <- round(mean(ref.period$rain),2)
ref.Tmin.mean <- round(mean(ref.period$tmin),2)
ref.Tmax.mean <- round(mean(ref.period$tmax),2)
ref.PE.mean <- round(mean(ref.period$PE),2)
ref.PPE.mean <- round(mean(ref.period$PPE),2)
ref.prec.95 <- round(quantile(ref.period$rain,0.95),2)
ref.prec.99 <- round(quantile(ref.period$rain,0.99),2)
ref.Tmin.10 <- round(quantile(ref.period$tmin,0.10),2)
ref.Tmin.90 <- round(quantile(ref.period$tmin,0.90),2)
ref.Tmin.95 <- round(quantile(ref.period$tmin,0.95),2)
ref.Tmax.10 <- round(quantile(ref.period$tmax,0.10),2)
ref.Tmax.90 <- round(quantile(ref.period$tmax,0.90),2)
ref.Tmax.95 <- round(quantile(ref.period$tmax,0.95),2)
Annual.Indices <- data %>%
  group_by(years) %>%
  summarise(
    warm.nights =sum(tmin > ref.Tmin.90),
    cool.nights = sum(tmin < ref.Tmin.10),
    warm.days = sum(tmax > ref.Tmax.90),
    cool.days = sum(tmax < ref.Tmax.10),
    SU = sum(tmax > 25),
    TR = sum(tmin > 20),
    PRCPTOT = sum(rain > 1),
    SDII = round(sum(rain)/PRCPTOT,2),
    R10 = round(sum(rain)/sum(rain > 10),2),
    R20 = round(sum(rain)/sum(rain > 20),2),
    R95 = round(sum(rain > ref.prec.95),2),
    R99 = round(sum(rain > ref.prec.99),2),
    CDD = max(rle(rain < 1)$lengths[rle(rain < 1)$values], na.rm = TRUE),
    CWD= max(rle(rain >= 1)$lengths[rle(rain >= 1)$values], na.rm = TRUE),
  )
reactable::reactable(Annual.Indices,
                     style = list(textAlign = "center"))

Monthly

data$month <- floor_date(data$dates, "month")
Monthly.Indices <- data %>%
  group_by(month) %>%
  summarise(
    TXx = round(max(tmax, na.rm = TRUE),1),
    TNx = round(max(tmin,na.rm = TRUE),1),
    TXn = round(min(tmax, na.rm = TRUE),1),
    TNn = round(min(tmin, na.rm = TRUE),1)
  )
Month <- as.numeric(format(Monthly$month, "%m"))
Monthly.Indices <- cbind(Monthly.Indices,
                         Monthly$Rainfall,
                         Monthly$Tmin,
                         Monthly$Tmax,
                         Monthly$Tavg,
                         Monthly$PE,
                         Monthly$PPE,
                         Month)
colnames(Monthly.Indices) <- c("Date","TXx","TNx","TXn","TNn",
                               "Rain","Tmin","Tmax","Tavg","PE","PPE","Month")
reactable::reactable(Monthly.Indices,
                     style = list(textAlign = "center"))

Non parametric trend tests: Annual scale

Considering the lack of knowledge regarding the distinct frequency distributions of each of the annual indices of Table 1, we applied the Mann-Kendall (MK) test to detected trends in their time series. Once a significant trend was detected we applied the Sen-Theil slope (Sen 1968) and (Theil 1950).

# Create the analysis matrix
Annual.Analysis <- as.data.frame(matrix(NA, 2, 18))
Annual.Analysis[2, ] <- "NoTrend"
rownames(Annual.Analysis) <- c("MK_pvalue", "Slope")
colnames(Annual.Analysis) <- c("warm.nights", "cool.nights", "warm.days", "cool.days",
                                "SU", "TR", "PRCPTOT", "SDII", "R10", "R20", 
                                "R95", "R99", "CDD", "CWD","Rain", "Tmin","Tmax","Tavg")

Annual.Indices <- cbind(Annual.Indices,Annual$Rainfall,Annual$Tmin,Annual$Tmax,Annual$Tavg)
Annual.trend <- apply(Annual.Indices[2:19],2,mk.test)
Annual.Analysis[1,] <- round(c(Annual.trend$warm.nights$p.value,
                        Annual.trend$cool.nights$p.value,
                        Annual.trend$warm.days$p.value,
                        Annual.trend$cool.days$p.value,
                        Annual.trend$SU$p.value,
                        Annual.trend$TR$p.value,
                        Annual.trend$PRCPTOT$p.value,
                        Annual.trend$SDII$p.value,
                        Annual.trend$R10$p.value,
                        Annual.trend$R20$p.value,
                        Annual.trend$R95$p.value,
                        Annual.trend$R99$p.value,
                        Annual.trend$CDD$p.value,
                        Annual.trend$CWD$p.value,
                        Annual.trend$`Annual$Rainfall`$p.value,
                        Annual.trend$`Annual$Tmin`$p.value,
                        Annual.trend$`Annual$Tmax`$p.value,
                        Annual.trend$`Annual$Tavg`$p.value),2)
for (i in 1:18){
  if (as.numeric(Annual.Analysis[1,i]) < 0.10) {
    sen <- sens.slope(as.matrix(Annual.Indices[,(i+1)]))
    Annual.Analysis[2,i] <- round(sen$estimates,2)
    if(as.numeric(Annual.Analysis[1,i]) < 0.001) {Annual.Analysis[1,i] <- "< 0.001"}
  }
}
# Convert to data frame for better table handling
Annual.Analysis <- as.data.frame(Annual.Analysis)

# Create and style table using gt
Annual.Analysis %>%
  t() %>%  # Transpose for better readability
  as.data.frame() %>%
  gt(rownames_to_stub = TRUE) %>%
  tab_header(
    title = "Annual Trend Analysis",
    subtitle = "Mann-Kendall Test and Sen's Slope"
  ) %>%
  tab_style(
    style = cell_text(align = "center"),  # Center-align text
    locations = cells_body(columns = everything())  # Apply to all body columns
  )

	MK_pvalue	Slope
Annual Trend Analysis
Mann-Kendall Test and Sen's Slope
warm.nights	< 0.001	0.48
cool.nights	< 0.001	-0.46
warm.days	< 0.001	0.34
cool.days	< 0.001	-0.2
SU	< 0.001	0.46
TR	< 0.001	0.4
PRCPTOT	0.77	NoTrend
SDII	0.15	NoTrend
R10	0.55	NoTrend
R20	0.15	NoTrend
R95	0.64	NoTrend
R99	0.22	NoTrend
CDD	0.33	NoTrend
CWD	0.79	NoTrend
Rain	0.36	NoTrend
Tmin	< 0.001	0.02
Tmax	< 0.001	0.02
Tavg	< 0.001	0.02

Non parametric trend tests: Monthly scale

Considering the lack of knowledge regarding the distinct frequency distributions of each of the monthly indices of Table 1, we applied the Mann-Kendall (MK) test to detected trends in their time series. Once a significant trend was detected we applied the Sen-Theil slope.

Monthly.Trend <- as.data.frame(matrix(NA,12,10))
Monthly.Slope <- as.data.frame(matrix(NA,12,10))
colnames(Monthly.Trend) <- c("Rain","Tmin","Tmax","Tavg","TXx","TNx","TXn","TNn","PE","PPE")
colnames(Monthly.Slope) <- c("Rain","Tmin","Tmax","Tavg","TXx","TNx","TXn","TNn","PE","PPE")

for (i in 1:12){
Month.trend <- apply(Monthly.Indices[which(Monthly.Indices$Month == i),2:11],2,mk.test)

Monthly.Trend[i,1:10] <- round(as.numeric(c(Month.trend$Rain$p.value,
                          Month.trend$Tmin$p.value,
                          Month.trend$Tmax$p.value,
                          Month.trend$Tavg$p.value,
                          Month.trend$TXx$p.value,
                          Month.trend$TNx$p.value,
                          Month.trend$TXn$p.value,
                          Month.trend$TNn$p.value,
                          Month.trend$PE$p.value,
                          Month.trend$PPE$p.value)),3)

Month.slope <- apply(Monthly.Indices[which(Monthly.Indices$Month == i),2:11],2,sens.slope)
Monthly.Slope[i,] <- round(as.numeric(c(Month.slope$Rain$estimates,
                                  Month.slope$Tmin$estimates,
                                  Month.slope$Tmax$estimates,
                                  Month.slope$Tavg$estimates,
                                  Month.slope$TXx$estimates,
                                  Month.slope$TNx$estimates,
                                  Month.slope$TXn$estimates,
                                  Month.slope$TNn$estimates,
                                  Month.slope$PE$estimates,
                                  Month.slope$PPE$estimates)),3)
Monthly.Slope[i,which(Monthly.Trend[i,] > 0.10)] <- c("NoChange")}
Monthly.Trend[Monthly.Trend < 0.001] <- "< 0.001"
# Convert to data frame for better table handling
Monthly.Trend <- as.data.frame(Monthly.Trend)

# Create and style table using gt
Monthly.Trend %>%
  as.data.frame() %>%
  gt(rownames_to_stub = TRUE) %>%
  tab_header(
    title = "Monthly Trend Analysis",
    subtitle = "Mann-Kendall test"
  ) %>%
  tab_style(
    style = cell_text(align = "center"),  # Center-align text
    locations = cells_body(columns = everything())  # Apply to all body columns
  )

	Rain	Tmin	Tmax	Tavg	TXx	TNx	TXn	TNn	PE	PPE
Monthly Trend Analysis
Mann-Kendall test
1	0.487	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	0.035	< 0.001	< 0.001	0.777
2	0.064	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	0.001	< 0.001	< 0.001	0.005
3	0.748	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	0.005	< 0.001	< 0.001	0.138
4	0.602	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	0.254
5	0.584	< 0.001	0.003	< 0.001	< 0.001	< 0.001	0.093	0.012	< 0.001	0.604
6	0.579	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	0.002	0.004	< 0.001	0.058
7	0.567	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	0.849	< 0.001	< 0.001	0.117
8	0.227	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	0.151	< 0.001	< 0.001	0.003
9	0.145	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	0.001
10	0.667	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	0.037
11	0.690	< 0.001	< 0.001	< 0.001	< 0.001	< 0.001	0.011	< 0.001	< 0.001	0.068
12	0.512	< 0.001	< 0.001	< 0.001	0.004	< 0.001	0.006	< 0.001	< 0.001	0.097

# Create and style table using gt
Monthly.Slope %>%
  as.data.frame() %>%
  gt(rownames_to_stub = TRUE) %>%
  tab_header(
    title = "Monthly Trend Analysis",
    subtitle = "Sen's Slope"
  ) %>%
  tab_style(
    style = cell_text(align = "center"),  # Center-align text
    locations = cells_body(columns = everything())  # Apply to all body columns
  )

	Rain	Tmin	Tmax	Tavg	TXx	TNx	TXn	TNn	PE	PPE
Monthly Trend Analysis
Sen's Slope
1	NoChange	0.018	0.014	0.017	0.011	0.016	0.008	0.021	0.198	NoChange
2	-0.36	0.017	0.019	0.018	0.015	0.015	0.015	0.022	0.191	-0.538
3	NoChange	0.018	0.017	0.017	0.017	0.018	0.013	0.019	0.190	NoChange
4	NoChange	0.021	0.016	0.018	0.014	0.017	0.016	0.018	0.163	NoChange
5	NoChange	0.018	0.008	0.013	0.013	0.014	0.008	0.014	0.092	NoChange
6	NoChange	0.019	0.015	0.017	0.012	0.017	0.017	0.019	0.107	-0.157
7	NoChange	0.023	0.012	0.017	0.015	0.020	NoChange	0.023	0.114	NoChange
8	NoChange	0.021	0.016	0.018	0.019	0.017	NoChange	0.029	0.138	-0.2
9	NoChange	0.025	0.026	0.025	0.024	0.026	0.018	0.027	0.210	-0.364
10	NoChange	0.027	0.024	0.026	0.023	0.027	0.022	0.033	0.263	-0.32
11	NoChange	0.021	0.015	0.017	0.012	0.018	0.011	0.023	0.187	-0.263
12	NoChange	0.019	0.016	0.018	0.009	0.016	0.014	0.026	0.215	-0.338

Goodness-of-Fit tests

The first step to apply any parametric trend this is to find suitable theoretical distributions for performing such probabilistic assessments. The Kolmogorov-Smirnov/Lilliefors Goodness-of-Fit test (Lilliefors, H. W. 1967) has widespread use in agro-hydro-meteorological studies where the parameters of the probability functions and the Goodness-of-Fit test are calculated from the same series. Different from the original version of the Kolmogorov-Smirnov test, the critical values for the Kolmogorov-Smirnov/Lilliefors test must be derived from statistical simulation techniques. Further information regarding this latter test can be found in Blain, 2014. For those cases where the Kolmogorov-Smirnov/Lilliefors test indicated two or more distributions as suitable for a particular series, we applied the second-order Akaike Information Criterium (AICc) to rank their performances.

Annual Scale

GoodFit.Annual <- matrix(NA,3,3)
GoodFit.Annual <- apply(Annual[,2:4], 2, test.Good)
rownames(GoodFit.Annual) <- c("Pvalue Normal","Pvalue Log-normal","Pvalue Gamma","Pvalue Weibull",
                             "AICc Normal","AICc Log-normal","AICc Gamma","AICc Weibull")
GoodFit.Annual <- as.data.frame(GoodFit.Annual)

# Create and style table using gt
GoodFit.Annual %>%
  as.data.frame() %>%
  gt(rownames_to_stub = TRUE) %>%
  tab_header(
    title = "Goodness-of-Fit"
  ) %>%
  tab_style(
    style = cell_text(align = "center"),  # Center-align text
    locations = cells_body(columns = everything())  # Apply to all body columns
  )

Variable	Distribution
Best Distribution
Pre	Gamma
Tmin	Normal
Tmax	Normal

Monthly Scale

GoodFit.Monthly <- list()
for (i in 1:12) {
  test <- apply(Monthly.Indices[which(Monthly.Indices$Month == i),6:8], 2, test.Good)
  rownames(test) <- c("Pvalue Normal","Pvalue Log-normal","Pvalue Gamma","Pvalue Weibull",
                             "AICc Normal","AICc Log-normal","AICc Gamma","AICc Weibull")
  GoodFit.Monthly[[i]] <- test}
names(GoodFit.Monthly) <- c("jan","feb","mar","apr","may","jun",
                            "jul","ago","sep","oct","nov","dec")

Month	Distribution
Precipitation
Jan	Gamma
Feb	Weibull
Mar	Weibull
Apr	Weibull
May	Weibull
Jun	Gamma
Jul	Weibull
Ago	Weibull
Sep	Weibull
Oct	Weibull
Nov	Gamma
Dec	Gamma

Month	Distribution
Tmin
Jan	Normal
Feb	Weibull
Mar	Weibull
Apr	Normal
May	Normal
Jun	Normal
Jul	Normal
Ago	Normal
Sep	Normal
Oct	LogNormal
Nov	LogNormal
Dec	Normal

Month	Distribution
Tmax
Jan	LogNormal
Feb	Gamma
Mar	Normal
Apr	LogNormal
May	Normal
Jun	LogNormal
Jul	Normal
Ago	Normal
Sep	LogNormal
Oct	LogNormal
Nov	LogNormal
Dec	Normal

Parametric trend Assessment: GAMLSS models and AICc

Nonstationary parametric distribution can be used as trend tests because they detect changes in the probabilistic structure of a time series (Blain et al., 2022). The idea behind this approach is to verify if nonstationary versions of a parametric distribution led to a better description of data variability in respect to those of the stationary version. If no significant improvement is found, the series is assumed to present no change in its frequency-distribution and no sign of climate change is detected. If a nonstationary version is found to be the best fitting model, the series is assumed to present significant changes in its frequency-distribution and the hypothesis of no climate change is not accepted.

The GAMLSS models considered in this study are described as follow.

Model 1 (stationary): The mean (µ) and dispersion (δ) of the distribution do not change over time.
Model 2: Only µ varies linearly over time.
Model 3: Only δ varies linearly over time.
Model 4: Both µ and δ vary linearly over time.
Model 5: Only µ varies non-linearly over time as a quadratic polynomial function.
Model 6: µ varies non-linearly over time as a quadratic polynomial function and δ vary linearly over time.

Annual

# Create the matrix
Best.Model.Annual <- matrix(NA, 3, 1)
rownames(Best.Model.Annual) <- c("Rain", "Tmin", "Tmax")
colnames(Best.Model.Annual) <- c("BestModelAnnual")
# Fill the matrix with selected models
Best.Model.Annual[1, 1] <- select.model(Annual$Rainfall, dist = "GA")
Best.Model.Annual[2, 1] <- select.model(Annual$Tmin, dist = "NO")
Best.Model.Annual[3, 1] <- select.model(Annual$Tmax, dist = "NO")
# Convert the matrix to a data frame
Best.Model.Annual_df <- as.data.frame(Best.Model.Annual)
Best.Model.Annual_df$Index <- rownames(Best.Model.Annual_df)
# Use gt to render the table
gt::gt(Best.Model.Annual_df) %>%
  gt::cols_label(Index = "Index", BestModelAnnual = "Best Model (Annual)") %>%
  gt::tab_header(title = "Parametric Trend Analysis (Annual)")  %>%
  tab_style(
    style = cell_text(align = "center"),  # Center-align text
    locations = cells_body(columns = everything())  # Apply to all body columns
  )

Best Model (Annual)	Index
Parametric Trend Analysis (Annual)
1	Rain
2	Tmin
2	Tmax

Monthly scale

Best.Model.Monthly <- matrix(NA,12,3)
rownames(Best.Model.Monthly) <- c("jan","fev","mar","apr","may","jun",
                                  "jul","ago","sep","oct","nov","dec")
colnames(Best.Model.Monthly) <- c("Rain","Tmin","Tmax")
#January
Best.Model.Monthly[1,1] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 1),6], dist = "GA")
Best.Model.Monthly[1,2] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 1),7], dist = "NO")
Best.Model.Monthly[1,3] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 1),8], dist = "LNO")
#Februay
Best.Model.Monthly[2,1] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 2),6], dist = "GA")
Best.Model.Monthly[2,2] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 2),7], dist = "WEI")
Best.Model.Monthly[2,3] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 2),8], dist = "NO")
#March
Best.Model.Monthly[3,1] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 3),6], dist = "WEI")
Best.Model.Monthly[3,2] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 3),7], dist = "WEI")
Best.Model.Monthly[3,3] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 3),8], dist = "NO")
#April
Best.Model.Monthly[4,1] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 4),6], dist = "WEI")
Best.Model.Monthly[4,2] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 4),7], dist = "NO")
Best.Model.Monthly[4,3] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 4),8], dist = "LNO")
#May
Best.Model.Monthly[5,1] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 5),6], dist = "WEI")
Best.Model.Monthly[5,2] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 5),7], dist = "NO")
Best.Model.Monthly[5,3] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 5),8], dist = "LNO")
#June
Best.Model.Monthly[6,1] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 6),6], dist = "GA")
Best.Model.Monthly[6,2] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 6),7], dist = "NO")
Best.Model.Monthly[6,3] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 6),8], dist = "LNO")
#July
Best.Model.Monthly[7,1] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 7),6], dist = "WEI")
Best.Model.Monthly[7,2] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 7),7], dist = "NO")
Best.Model.Monthly[7,3] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 7),8], dist = "NO")
#Ago
Best.Model.Monthly[8,1] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 8),6], dist = "WEI")
Best.Model.Monthly[8,2] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 8),7], dist = "NO")
Best.Model.Monthly[8,3] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 8),8], dist = "NO")
#Sep
Best.Model.Monthly[9,1] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 9),6], dist = "WEI")
Best.Model.Monthly[9,2] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 9),7], dist = "NO")
Best.Model.Monthly[9,3] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 9),8], dist = "LNO")
#Oct
Best.Model.Monthly[10,1] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 10),6], dist = "WEI")
Best.Model.Monthly[10,2] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 10),7], dist = "LNO")
Best.Model.Monthly[10,3] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 10),8], dist = "LNO")
#Nov
Best.Model.Monthly[11,1] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 11),6], dist = "GA")
Best.Model.Monthly[11,2] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 11),7], dist = "LNO")
Best.Model.Monthly[11,3] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 11),8], dist = "NO")
#Dec
Best.Model.Monthly[12,1] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 12),6], dist = "GA")
Best.Model.Monthly[12,2] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 12),7], dist = "NO")
Best.Model.Monthly[12,3] <- select.model(Monthly.Indices[which(Monthly.Indices$Month == 12),8], dist = "NO")
# Convert the matrix to a data frame
Best.Model.Monthly_df <- as.data.frame(Best.Model.Monthly)
Best.Model.Monthly_df$Month <- rownames(Best.Model.Monthly_df)
# Use gt to render the table
gt::gt(Best.Model.Monthly_df) %>%
  gt::cols_label(Month = "Month", Rain = "Rain Model", Tmin = "Tmin Model", Tmax = "Tmax Model") %>%
  gt::tab_header(title = "Parametric Trend Analysis (Monthly)") %>%
  tab_style(
    style = cell_text(align = "center"),  # Center-align text
    locations = cells_body(columns = everything())  # Apply to all body columns
  )

Rain Model	Tmin Model	Tmax Model	Month
Parametric Trend Analysis (Monthly)
1	2	2	jan
1	4	2	fev
1	4	2	mar
1	2	2	apr
1	2	2	may
1	2	2	jun
1	2	5	jul
1	2	2	ago
1	2	2	sep
1	2	2	oct
1	2	2	nov
1	4	5	dec

Final Remarks and Data Download

This document, generated in the R-software environment, is part of a broader study, which has been submitted to publication in the journal Bragantia.
The study was based on the hypothesis that long-term continuous meteorological records provide unique and valuable insights into climate evolution. This is the reason why this document made publicly available one of the longest meteorological series in South America: The Agronomic Institute’s historical meteorological dataset.
As described in this document, rigorous data-quality tests were applied, confirming the dataset’s high reliability.
The trend analyses revealed significant indicators of climate change, with a pronounced shift toward warmer conditions in both minimum and maximum air temperature records.
The increased air temperatures resulted in significant upward trends in potential evapotranspiration rates (a measure of atmospheric water demand). Consequently, a notable decrease in the difference between precipitation and potential evapotranspiration was observed, strongly suggesting intensified regional drought conditions.
This drying trend was particularly evident in February, June, August, September, October, November, and December, indicating a lengthening of the local dry season.