Taller Estuario datos USGS

Author

Denny S. Fernández del Viso

Eliminar messages and warnings

# eliminar messages and warnings
options(warn=-1, message=-1, echo=TRUE)

Descripción del lugar

Quebrada Josefina

Insertar mapa donde se encuentra la estación de monitoreo

# insertar mapa
library(leaflet)
m <- leaflet() %>%
  addTiles() %>%
  setView(lng = -66.07517, lat = 18.40717, zoom = 100)
m

Quebrada Josefina

Variables de interés

Discharge (ft³/s): código USGS: 00060
Gage height (ft): código USGS: 00065

Usando el paquete `dataRetrieval` para descargar datos de la estación de monitoreo

# activación del paquete
library(dataRetrieval)

# parámetros de la descarga
siteNo <- "50049310"
pCode <- c("00060","00065")
start.date <- "1989-01-01"
end.date <- "2024-08-31"

# comando para descargar
josefina <- readNWISdata(siteNumber = siteNo, parameterCd = pCode, startDate = start.date, endDate = end.date)

# inicio de los datos
head(josefina)

  agency_cd  site_no   dateTime X_..2.._00060_00003 X_..2.._00060_00003_cd
1      USGS 50049310 1989-01-01                 3.0                      A
2      USGS 50049310 1989-01-02                 3.4                      A
3      USGS 50049310 1989-01-03                 4.7                      A
4      USGS 50049310 1989-01-04                 4.3                      A
5      USGS 50049310 1989-01-05                 2.8                      A
6      USGS 50049310 1989-01-06                 4.1                      A
  X_00060_00003 X_00060_00003_cd X_00065_00003 X_00065_00003_cd tz_cd
1          3.10                A            NA             <NA>   UTC
2          3.44                A            NA             <NA>   UTC
3          4.71                A            NA             <NA>   UTC
4          4.31                A            NA             <NA>   UTC
5          2.88                A            NA             <NA>   UTC
6          4.16                A            NA             <NA>   UTC

# fin de los datos
tail(josefina)

     agency_cd  site_no   dateTime X_..2.._00060_00003 X_..2.._00060_00003_cd
2801      USGS 50049310 2024-08-26                  NA                   <NA>
2802      USGS 50049310 2024-08-27                  NA                   <NA>
2803      USGS 50049310 2024-08-28                  NA                   <NA>
2804      USGS 50049310 2024-08-29                  NA                   <NA>
2805      USGS 50049310 2024-08-30                  NA                   <NA>
2806      USGS 50049310 2024-08-31                  NA                   <NA>
     X_00060_00003 X_00060_00003_cd X_00065_00003 X_00065_00003_cd tz_cd
2801          9.09                P          1.53                P   UTC
2802          7.74                P          1.51                P   UTC
2803          8.95                P          1.52                P   UTC
2804          9.67                P          1.53                P   UTC
2805         21.80                P          1.62                P   UTC
2806         13.70                P          1.56                P   UTC

Variables obtenidas

# variables obtenidas
names(josefina)

 [1] "agency_cd"              "site_no"                "dateTime"              
 [4] "X_..2.._00060_00003"    "X_..2.._00060_00003_cd" "X_00060_00003"         
 [7] "X_00060_00003_cd"       "X_00065_00003"          "X_00065_00003_cd"      
[10] "tz_cd"

Vamos a renombrar las variables para facilitar su uso y obtener la metadata del archivo.

# renombrar las variables
josefina <- renameNWISColumns(josefina)

# atributos del archivo
names(attributes(josefina))

[1] "names"         "row.names"     "url"           "siteInfo"     
[5] "variableInfo"  "disclaimer"    "statisticInfo" "queryTime"    
[9] "class"

# ver atributos de la metadata
attr(josefina, "siteInfo")

                              station_nm  site_no agency_cd timeZoneOffset
1 QUEBRADA JOSEFINA AT PINERO AVENUE, PR 50049310      USGS         -04:00
  timeZoneAbbreviation dec_lat_va dec_lon_va       srs siteTypeCd    hucCd
1                  AST   18.40717  -66.07517 EPSG:4326         ST 21010005
  stateCd countyCd network
1      72    72127    NWIS

attr(josefina, "names")

 [1] "agency_cd"     "site_no"       "dateTime"      "..2.._Flow"   
 [5] "..2.._Flow_cd" "Flow"          "Flow_cd"       "GH"           
 [9] "GH_cd"         "tz_cd"

attr(josefina, "variableInfo")

  variableCode           variableName              variableDescription
1        00060 Streamflow, ft&#179;/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          NA

attr(josefina, "statisticInfo")

  statisticCd statisticName
1       00003          Mean

Gráficos de los datos

# gráfico de la descarga
library(ggplot2)
ts <- ggplot(data = josefina,
             aes(dateTime, Flow)) +
      geom_line()
ts

Hay un gran “gap” en los datos, por lo que vamos a seleccionar un rango de fechas más recientes y con mediciones continuas para realizar un análisis de serie de tiempo.

Seleccionar nuevos datos y graficar

# seleccionar datos
josefina2 <- josefina[josefina$dateTime > "2018-11-21",]

# gráfico de la descarga
ts2 <- ggplot(data = josefina2,
             aes(dateTime, Flow)) +
      geom_line()
ts2

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 de estacionariedad de la serie de tiempo.

# análisis de serie de tiempo
library(tseries)

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo

# Perform ADF test
adf_result <- adf.test(josefina2$Flow)

# Print ADF test results
print(adf_result)


    Augmented Dickey-Fuller Test

data:  josefina2$Flow
Dickey-Fuller = -10.746, Lag order = 12, 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)
josefina2_ts <- ts(josefina2$Flow, frequency = 365)

# descomposición de la serie de tiempo
josefina2_decomp <- decompose(josefina2_ts)

# gráfico de la descomposición con Time = año
plot(josefina2_decomp, type = "l", xaxt = "n", xlab = "Año")
axis(side = 1, at = 1:7, labels = c(2019, 2020, 2021, 2022, 2023, 2024,2025), cex.axis = 0.7)

# análisis de la tendencia
josefina2_trend <- josefina2_decomp$trend
plot(josefina2_trend, type = "l", xaxt = "n", xlab = "Año")
axis(side = 1, at = 1:7, labels = c(2019, 2020, 2021, 2022, 2023, 2024,2025), cex.axis = 0.7)

# análisis de la estacionalidad
josefina2_seasonal <- josefina2_decomp$seasonal
plot(josefina2_seasonal, type = "l", xaxt = "n", xlab = "Año")
axis(side = 1, at = 1:7, labels = c(2019, 2020, 2021, 2022, 2023, 2024,2025), cex.axis = 0.7)

# análisis de la estacionalidad
josefina2_random <- josefina2_decomp$random
plot(josefina2_random, type = "l", xaxt = "n", xlab = "Año")
axis(side = 1, at = 1:7, labels = c(2019, 2020, 2021, 2022, 2023, 2024,2025), cex.axis = 0.7)

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(josefina2, file = "josefina2.csv")

Otra versión de análisis de serie de tiempo

Incluye un análisis de autocorrelación y autocorrelación parcial.

# Load required libraries
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ lubridate 1.9.3     ✔ tibble    3.2.1
✔ purrr     1.0.2     ✔ tidyr     1.3.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(lubridate)
library(forecast)

# Read the CSV file
data <- read.csv("josefina2.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 = 180)

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

Otro código

Similar a las gráficas anteriores, incluye la línea de tendencia y la serie de tiempo original. También se realiza un ajuste estacional.

# Load necessary libraries
library(tidyverse)
library(lubridate)
library(tseries)
library(forecast)

# Read the CSV file
data <- read.csv("josefina2.csv")

# Convert dateTime to Date type
data$dateTime <- as.Date(data$dateTime)

# Create a time series object for Flow
flow_ts <- ts(data$Flow, frequency = 365, start = c(year(min(data$dateTime)), yday(min(data$dateTime))))

# Perform time series decomposition
flow_decomp <- stl(flow_ts, s.window = "periodic")

# Plot the decomposition
plot(flow_decomp)

# Extract components
trend <- flow_decomp$time.series[, "trend"]
seasonal <- flow_decomp$time.series[, "seasonal"]
remainder <- flow_decomp$time.series[, "remainder"]

# Calculate some basic statistics
cat("Trend range:", range(trend), "\n")

Trend range: 9.13621 21.56197

cat("Seasonal range:", range(seasonal), "\n")

Seasonal range: -10.33703 118.2158

cat("Remainder range:", range(remainder), "\n")

Remainder range: -127.2727 420.6568

# Plot the original series with the trend
plot(flow_ts, main = "Flow Time Series with Trend", ylab = "Flow")
lines(trend, col = "red", lwd = 2)
legend("topright", legend = c("Original", "Trend"), col = c("black", "red"), lty = 1, lwd = c(1, 2))

# Plot the seasonal component
plot(seasonal, main = "Seasonal Component of Flow", ylab = "Seasonal Effect")

# Perform seasonal adjustment
seasonally_adjusted <- flow_ts - seasonal
plot(seasonally_adjusted, main = "Seasonally Adjusted Flow", ylab = "Adjusted Flow")
lines(trend, col = "red", lwd = 2)
legend("topright", legend = c("Seasonally Adjusted", "Trend"), col = c("black", "red"), lty = 1, lwd = c(1, 2))

Descripción del lugar

Quebrada Josefina

Insertar mapa donde se encuentra la estación de monitoreo

Quebrada Josefina

Variables de interés

Usando el paquete dataRetrieval para descargar datos de la estación de monitoreo

Variables obtenidas

Gráficos de los datos

Seleccionar nuevos datos y graficar

Análisis de Serie de Tiempo

Análisis de estacionalidad, tendencia y predicción

Creación de archivo csv

Otra versión de análisis de serie de tiempo

Otro código

Usando el paquete `dataRetrieval` para descargar datos de la estación de monitoreo