# eliminar messages and warnings
options(warn=-1, message=-1, echo=TRUE)Taller Estuario - Análisis de Series de Tiempo con datos Hidrológicos del USGS
¿Para qué sirve el análisis de series de tiempo?
El análisis de series de tiempo es una técnica estadística que se utiliza para analizar y modelar datos que varían con el tiempo. Se utiliza para comprender patrones, tendencias y ciclos en los datos a lo largo del tiempo. En el contexto de datos hidrológicos, el análisis de series de tiempo puede ayudar a predecir caudales, identificar patrones estacionales y evaluar la variabilidad de los flujos de agua en ríos y cuencas hidrográficas.
En un contexto más aplicado, el análisis de series de tiempo puede ayudar a los gestores de recursos hídricos a tomar decisiones informadas sobre la gestión de los recursos hídricos, la planificación de infraestructuras hidráulicas y la respuesta a eventos extremos como inundaciones y sequías, así como al cambio climático.
Objetivo del taller
El objetivo de este taller es realizar un análisis de series de tiempo de datos hidrológicos del USGS (Servicio Geológico de los Estados Unidos) para una estación de monitoreo en el Río Piedras en Hato Rey, Puerto Rico. Vamos a descargar los datos, visualizarlos, realizar un análisis de estacionalidad y tendencia, y predecir los flujos de agua futuros utilizando técnicas de series de tiempo.
Estación de monitoreo Hato Rey del Río Piedras
La estación Hato Rey de monitoreo del Río Piedras se encuentra en la intersección del río con la PR-17:
# insertar mapa
library(leaflet)
m <- leaflet() %>%
addTiles() %>%
setView(lng = -66.06894, lat = 18.40755, zoom = 100)
mVariables de interés
A partir de la estación de monitoreo Hato Rey del Río Piedras, se obtendrán los siguientes datos:
Discharge (ft³/s): código USGS: 00060
Gage height (ft): código USGS: 00065
Los códigos de parámetros y estadísticos se obtienen de la páginas de USGS:
Uso del paquete dataRetrieval para descargar datos de la estación de monitoreo
El paquete dataRetrieval en R permite descargar datos hidrológicos del USGS de forma sencilla. Vamos a utilizar este paquete para descargar los datos de la estación de monitoreo Hato Rey del Río Piedras.
# activación del paquete
library(dataRetrieval)
# parámetros de la descarga
siteNo <- "50049100"
pCode <- c("00060","00065")
start.date <- "1994-01-28"
end.date <- "2024-08-31"
# comando para descargar
hato_rey <- readNWISdata(siteNumber = siteNo, parameterCd = pCode, startDate = start.date, endDate = end.date)
# inicio de los datos
head(hato_rey) agency_cd site_no dateTime X_00060_00003 X_00060_00003_cd X_00065_00003
1 USGS 50049100 1994-01-28 20 A e NA
2 USGS 50049100 1994-01-29 75 A 5.74
3 USGS 50049100 1994-01-30 14 A 5.44
4 USGS 50049100 1994-01-31 13 A 5.42
5 USGS 50049100 1994-02-01 17 A 5.45
6 USGS 50049100 1994-02-02 15 A 5.44
X_00065_00003_cd tz_cd
1 <NA> UTC
2 A UTC
3 A UTC
4 A UTC
5 A UTC
6 A UTC# fin de los datos
tail(hato_rey) agency_cd site_no dateTime X_00060_00003 X_00060_00003_cd
11169 USGS 50049100 2024-08-26 27.5 P
11170 USGS 50049100 2024-08-27 22.0 P
11171 USGS 50049100 2024-08-28 23.6 P
11172 USGS 50049100 2024-08-29 28.6 P
11173 USGS 50049100 2024-08-30 47.5 P e
11174 USGS 50049100 2024-08-31 23.8 P e
X_00065_00003 X_00065_00003_cd tz_cd
11169 4.02 P UTC
11170 3.95 P UTC
11171 3.97 P UTC
11172 4.04 P UTC
11173 NA <NA> UTC
11174 NA <NA> UTCVariables obtenidas
Vamos a ver las variables obtenidas en el archivo descargado.
# variables obtenidas
names(hato_rey)[1] "agency_cd" "site_no" "dateTime" "X_00060_00003"
[5] "X_00060_00003_cd" "X_00065_00003" "X_00065_00003_cd" "tz_cd" Vamos a renombrar las variables para facilitar su uso y obtener la metadata del archivo.
# renombrar las variables
hato_rey <- renameNWISColumns(hato_rey)Metadata del archivo
Vamos a ver la metadata del archivo descargado para obtener información sobre la estación de monitoreo, las variables y los estadísticos utilizados.
# atributos del archivo
names(attributes(hato_rey))[1] "names" "row.names" "url" "siteInfo"
[5] "variableInfo" "disclaimer" "statisticInfo" "queryTime"
[9] "class" # ver atributos de la metadata
attr(hato_rey, "siteInfo") station_nm site_no agency_cd timeZoneOffset
1 RIO PIEDRAS AT HATO REY, PR 50049100 USGS -04:00
timeZoneAbbreviation dec_lat_va dec_lon_va srs siteTypeCd hucCd
1 AST 18.40755 -66.06894 EPSG:4326 ST 21010005
stateCd countyCd network
1 72 72127 NWISattr(hato_rey, "names")[1] "agency_cd" "site_no" "dateTime" "Flow" "Flow_cd" "GH"
[7] "GH_cd" "tz_cd" attr(hato_rey, "variableInfo") variableCode variableName variableDescription
1 00060 Streamflow, ft³/s Discharge, cubic feet per second
2 00065 Gage height, ft Gage height, feet
valueType param_unit options noDataValue
1 Derived Value ft3/s Mean NA
2 Derived Value ft Mean NAattr(hato_rey, "statisticInfo") statisticCd statisticName
1 00003 MeanGráficos de serie de tiempo de los datos
Usamos el paquete ggplot2 para realizar un gráfico de la serie de tiempo de la descarga.
# Load required libraries
library(ggplot2)
library(lubridate) # For handling date-time data
# Assuming 'hato_rey' is your dataframe
# First, ensure the 'dateTime' column is in the correct format
hato_rey$dateTime <- as.POSIXct(hato_rey$dateTime)
# Create the base plot
ts <- ggplot(data = hato_rey,
aes(dateTime, Flow)) +
geom_line() +
xlab("Año") +
ylab("Descarga (ft³/s)")
# Add more ticks on x-axis
ts + scale_x_datetime(
date_breaks = "4 years",
date_labels = "%Y",
date_minor_breaks = "1 year"
)
# save plot
ggsave("hato_rey_ts_disch.png")Análisis de Serie de Tiempo
Para un primer análisis de serie de tiempo, vamos a utilizar el paquete tseries para realizar una prueba para saber si la serie de tiempo es estacionaria (sin tendencia) o no. La prueba Dickey-Fuller aumentada (ADF) es una prueba estadística que determina si una serie temporal es estacionaria.
# análisis de serie de tiempo
library(tseries)
# Perform ADF test
adf_result <- adf.test(hato_rey$Flow)
# Print ADF test results
print(adf_result)
Augmented Dickey-Fuller Test
data: hato_rey$Flow
Dickey-Fuller = -17.309, Lag order = 22, p-value = 0.01
alternative hypothesis: stationary# Interpret results
cat("\nInterpretation:\n")
Interpretation:if (adf_result$p.value < 0.05) {
cat("The time series is stationary (reject the null hypothesis).\n")
} else {
cat("The time series is non-stationary (fail to reject the null hypothesis).\n")
}The time series is stationary (reject the null hypothesis).Análisis de estacionalidad, tendencia y predicción
En esta sección utilizaremos el paquete forecast para realizar un análisis de estacionalidad, tendencia y predicción de la serie de tiempo.
# análisis de estacionalidad y tendencia
library(forecast)
# periodo de estacionalidad anual (365 días)
hato_rey_ts <- ts(hato_rey$Flow, frequency = 365)
# descomposición de la serie de tiempo
hato_rey_decomp <- decompose(hato_rey_ts)
# Create a vector of years
years <- 1994:2023
# gráfico de la descomposición con Time = año
plot(hato_rey_decomp, type = "l", xaxt = "n", xlab = "Año")
axis(side = 1, at = 1:30, labels = years, cex.axis = 0.7)
Estacionalidad promedio anual
Vamos a calcular la estacionalidad promedio anual de la serie de tiempo y visualizarla en un gráfico.
monthplot(hato_rey_ts, xlab = "Time (days)", ylab = "Caudal (ft³/s)")
seasonplot(hato_rey_ts, year.labels = TRUE, main = "", type = "l")
Promedio mensual de la serie de tiempo
Vamos a calcular el promedio mensual de la serie de tiempo y visualizarlo en un gráfico, con boxplots y puntos de media por mes.
# Cargar las bibliotecas necesarias
library(ggplot2)
library(lubridate)
library(dplyr)
# Convertir hato_rey_ts a un dataframe
start_date <- as.Date("1994-01-01") # Asumimos que los datos comienzan en 1994
dates <- seq(start_date, by="day", length.out=length(hato_rey_ts))
df <- data.frame(
fecha = dates,
caudal = as.vector(hato_rey_ts)
)
# Verificar el rango de fechas
print(paste("Rango de fechas:", min(df$fecha), "a", max(df$fecha)))[1] "Rango de fechas: 1994-01-01 a 2024-08-04"# Extraer el mes y calcular la media mensual
df <- df %>%
mutate(
mes = month(fecha, label = TRUE, abbr = TRUE),
anio = year(fecha)
) %>%
group_by(mes, anio) %>%
summarise(caudal_medio = mean(caudal, na.rm = TRUE), .groups = 'drop')
# Verificar los años únicos en el dataset
print("Años únicos en el dataset:")[1] "Años únicos en el dataset:"print(sort(unique(df$anio))) [1] 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008
[16] 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023
[31] 2024# Encontrar el año con el caudal máximo para cada mes
max_caudal_por_mes <- df %>%
group_by(mes) %>%
slice(which.max(caudal_medio)) %>%
ungroup()
# Imprimir los años con caudal máximo por mes
print("Años con caudal máximo por mes:")[1] "Años con caudal máximo por mes:"print(max_caudal_por_mes)# A tibble: 12 × 3
mes anio caudal_medio
<ord> <dbl> <dbl>
1 Jan 2022 103.
2 Feb 1995 77.6
3 Mar 2007 83.6
4 Apr 2004 216.
5 May 2011 127.
6 Jun 2017 125.
7 Jul 2011 257.
8 Aug 1996 246.
9 Sep 2017 144.
10 Oct 2001 165.
11 Nov 2007 110.
12 Dec 2009 96.1# Crear el gráfico
p <- ggplot(df, aes(x = mes, y = caudal_medio, group = mes)) +
geom_boxplot(fill = "lightblue", color = "blue", alpha = 0.7) +
geom_jitter(width = 0.2, alpha = 0.5, color = "darkblue") +
stat_summary(fun = mean, geom = "point", shape = 23, size = 3, fill = "red") +
geom_text(data = max_caudal_por_mes,
aes(label = anio, y = caudal_medio),
vjust = -0.5, color = "red", fontface = "bold") +
labs(x = "Mes", y = "Caudal medio (ft³/s)") +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 0, hjust = 0.5),
legend.position = "none",
plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5, size = 10, color = "darkgrey")
)
p
# save plot
ggsave("hato_rey_caudal_mes.png")Predicción con la serie de tiempo
Creación de archivo csv
Necesitamos guardar los datos en un archivo csv para poder realizar otros análisis de serie de tiempo en R.
write.csv(hato_rey, file = "hato_rey.csv")Versión de análisis de serie de tiempo con STL y ETS
ETS (Error, Tendencia, Estacional) y STL (Descomposición Estacional y de Tendencia usando Loess)
ETS: se utiliza principalmente para realizar pronósticos STL: se utiliza principalmente para descomponer y comprender los componentes de la serie
Incluye un análisis de autocorrelación y autocorrelación parcial.
# Load required libraries
library(tidyverse)
library(lubridate)
library(forecast)
# Read the CSV file
data <- read.csv("hato_rey.csv")
# Convert dateTime to Date type and Flow to numeric
data$dateTime <- as.Date(data$dateTime)
data$Flow <- as.numeric(data$Flow)
# Create a time series object
ts_data <- ts(data$Flow, frequency = 365, start = c(year(min(data$dateTime)), yday(min(data$dateTime))))
# Decompose the time series
decomposed <- stl(ts_data, s.window = "periodic")
# Plot the decomposition
# plot(decomposed)
# Forecast using STL + ETS(A,N,N)
forecast_result <- forecast(decomposed, method = "ets", etsmodel = "ANN", h = 780)
# Plot the forecast
plot(forecast_result)
# Plot ACF and PACF
par(mfrow = c(1, 1))
acf(ts_data, main = "Autocorrelation Function (ACF)")
pacf(ts_data, main = "Partial Autocorrelation Function (PACF)")
# Reset plot layout
par(mfrow = c(1, 1))Interpretación de los gráficos de autocorrelación
ACF mide la correlación entre una serie temporal y sus valores rezagados (“lags”). Muestra cómo se relaciona un punto de datos con sus valores pasados en distintos rezagos temporales. La disminución lenta de la ACF sugiere una tendencia en los datos. En el caso de las series estacionarias, la ACF debería caer rápidamente a cero.
Los picos significativos de PACF indican qué rezagos son más informativos para predecir el valor actual.
Cálculo de la descarga de inundación de 100 años
El cálculo de la descarga de inundación de 100 años es un análisis común en hidrología para evaluar la probabilidad de que ocurra una inundación de cierta magnitud en un período de 100 años. Vamos a utilizar el método de USGS para calcular la descarga de inundación de 100 años para la estación de monitoreo Hato Rey del Río Piedras.
# Install and load required packages
library(dataRetrieval)
library(fitdistrplus)
library(evd)
library(dplyr)
# Set the USGS site number (example: Potomac River at Point of Rocks, MD)
site_number <- "50049100"
# Fetch annual peak streamflow data
peak_data <- readNWISpeak(siteNumbers = site_number)
# Prepare data for analysis
peak_flows <- peak_data$peak_va
# Fit a log-normal distribution
fit_ln <- fitdist(peak_flows, "lnorm")
# Calculate the 100-year flood (1% annual exceedance probability)
flood_100yr <- qlnorm(0.99, meanlog = fit_ln$estimate[1], sdlog = fit_ln$estimate[2])
# Print the result
cat("The estimated 100-year flood discharge is", flood_100yr, "cubic feet per second.\n")The estimated 100-year flood discharge is 16072.88 cubic feet per second.# Calculate confidence intervals using the delta method
ci <- confint(fit_ln, level = 0.95)
ci_100yr <- qlnorm(0.99, meanlog = ci[,1], sdlog = ci[,2])
cat("95% Confidence Interval:", min(ci_100yr), "to", max(ci_100yr), "cubic feet per second.\n")95% Confidence Interval: 4.693922 to 4.131216e+12 cubic feet per second.# Plot the frequency curve
return_periods <- c(2, 5, 10, 25, 50, 100, 200, 500)
exceedance_probs <- 1 - 1 / return_periods
plot_data <- data.frame(
Return_Period = return_periods,
Discharge = qlnorm(exceedance_probs, meanlog = fit_ln$estimate[1], sdlog = fit_ln$estimate[2])
)
plot(plot_data$Return_Period, plot_data$Discharge,
log = "x", type = "l", col = "black",
xlab = "Return Period (years)",
ylab = "Discharge (cfs)",
main = "Flood Frequency Curve"
)