Curso-Taller: Análisis Hidrológico de Cuencas con R
Clase 2: Procesamiento de Datos Climáticos Grillados (parte 1)
Harold Llauca
📦 1. Preparación del entorno de trabajo
Se configura el entorno en R, definiendo el directorio de trabajo y cargando los paquetes necesarios para asegurar un flujo de trabajo ordenado y reproducible.
Directorio de trabajo
Instalar y cargar paquetes
# Establecer mirror de CRAN
options(repos = c(CRAN = "https://cloud.r-project.org"))
# Lista de paquetes necesarios
packages <- c(
"terra", # manejo de datos raster y vectoriales
"sf", # datos vectoriales espaciales
"exactextractr", # extracción de valores por polígonos
"zoo", # manejo de series de tiempo
"hydroTSM", # agregación temporal (diario, mensual, anual)
"openxlsx", # lectura y escritura de archivos Excel
"ggplot2", # visualización de datos
"tidyr" # transformación a formato largo
)
# Instalar paquetes si no están instalados
for (pkg in packages) {
if (!requireNamespace(pkg, quietly = TRUE)) {
install.packages(pkg)
}
}
# Cargar paquetes
invisible(lapply(packages, function(pkg) {
suppressPackageStartupMessages(library(pkg, character.only = TRUE))
}))## Warning: package 'ggplot2' was built under R version 4.5.2🗺️ 2. Datos espaciales y visualización
Se visualiza la zona de estudio para identificar la ubicación y extensión de las subcuencas.
# Lectura de datos espaciales
cat <- vect('./Data/SIG/subbasins.shp') # subcuencas (polígonos)
dem <- rast('./Data/SIG/dem.tif') # elevación (raster)## class : SpatVector
## geometry : polygons
## dimensions : 9, 3 (geometries, attributes)
## extent : -75.72391, -75.04077, -13.74518, -12.87817 (xmin, xmax, ymin, ymax)
## source : subbasins.shp
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## names : Area COMID Region
## type : <num> <num> <chr>
## values : 226.8 501 A
## 298.6 510 A
## 597.8 521 A## class : SpatRaster
## size : 3120, 2459, 1 (nrow, ncol, nlyr)
## resolution : 0.0002777778, 0.0002777778 (x, y)
## extent : -75.72403, -75.04097, -13.74486, -12.87819 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : dem.tif
## name : file32784a1815ed
## min value : 730
## max value : 5208## [1] "GEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]]"## [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n MEMBER[\"World Geodetic System 1984 (G2296)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"# Visualización
plot(dem,
main = "Subcuencas sobre el DEM",
col = terrain.colors(50))
plot(cat, add = TRUE,
border = 'black',
lwd = 1)
text(cat, cat$COMID)
🌧️ 3. Datos climáticos grillados
3.1 Exploración de datos PISCOp v2.1
PISCO (Peruvian Interpolated data of SENAMHI’s Climatological and Hydrological Observations) es un producto climático grillado que combina datos de estaciones con interpolación espacial, proporcionando información continua de precipitación a escala diaria y alta resolución (~0.1°).
Descarga archivo NetCDF en: https://iridl.ldeo.columbia.edu/SOURCES/.SENAMHI/.HSR/.PISCO/.Prec/.v2p1/.stable/.daily/.Prec/datafiles.html.
Nota: El archivo se encuentra en formato NetCDF, que permite almacenar datos multidimensionales organizados en dimensiones como tiempo, latitud y longitud.
# Secuencia de fechas del producto (diario)
dates <- seq(as.Date('1981-01-01'),
as.Date('2016-12-31'),
by = 'day')# Leer archivo NetCDF como raster multicapa (cada capa = 1 día)
pisco <- rast('./Data/PISCO/PISCOp_v2.1.nc')## class : SpatRaster
## size : 198, 133, 13149 (nrow, ncol, nlyr)
## resolution : 0.09999998, 0.1 (x, y)
## extent : -81.3, -68, -18.8, 1 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
## source : PISCOp_v2.1.nc
## varname : Prec (Precipitation)
## names : Prec_~671.5, Prec_~672.5, Prec_~673.5, Prec_~674.5, Prec_~675.5, Prec_~676.5, ...
## unit : mm/day
## depth : 7671.5 to 20819.5 (T [days since 1960-01-01]: 13149 steps)## [1] 13149## [1] 198 133 13149## [1] TRUEExplorar raster multicapa
# Visualizar una capa específica (ejemplo: día 100)
plot(pisco[[100]], alpha=0.9,
main = "Precipitación diaria (capa 100)")
Recortar el raster a la zona de estudio (subcuencas)
# Se usa un factor 1.1 para ampliar ligeramente la extensión para evitar que el recorte quede muy ajustado al borde
pisco_cat <- crop(pisco, 1.1 * ext(cat))
# Comparar tamaño del raster antes y después del recorte
dim(pisco)## [1] 198 133 13149## [1] 10 8 13149## SpatExtent : -81.2999999884403, -68.0000030633175, -18.7999999999697, 0.999999988048815 (xmin, xmax, ymin, ymax)## SpatExtent : -75.8000012600061, -75.0000014449611, -13.8000000029954, -12.8000000036005 (xmin, xmax, ymin, ymax)# Visualizar una capa del raster recortado
plot(pisco_cat[[100]], alpha=0.9,
main = "Precipitación recortada (capa 100)")
# Superponer las subcuencas
plot(cat, add = TRUE,
border = 'magenta',
lwd = 2)
🧮 4. Procesamiento espacio-temporal
4.1 Extracción por subcuenca
A partir del raster multicapa se calcula la precipitación media en cada subcuenca, obteniendo datos espacio-temporales.
# Convertir a formato compatible y extraer promedio por polígono
pr_cat <- t(exact_extract(pisco_cat,
as(cat, 'Spatial'),
fun = 'mean',
weights ='area'))## | | | 0% | |======== | 11% | |================ | 22% | |======================= | 33% | |=============================== | 44% | |======================================= | 56% | |=============================================== | 67% | |====================================================== | 78% | |============================================================== | 89% | |======================================================================| 100%# Extraer ID de cada subcuenca (ej. COMID)
cat_id <- cat$COMID
# Asignar nombres de columnas a la matriz de precipitación
colnames(pr_cat) <- cat_id
# Nota:
# Cada columna ahora representa una subcuenca identificada por su COMID## 501 510 521 525 527 533
## [1,] 0.0000000 0.000000 0.000000000 0.00000000 0.0000000 0.0000000
## [2,] 3.6393113 4.073518 0.552590847 0.43440908 0.6454489 0.5282456
## [3,] 4.2326713 4.979711 0.873725295 2.31732988 3.5731442 1.7681453
## [4,] 5.4923182 9.815919 0.947537184 1.89248145 3.7520499 5.3061824
## [5,] 0.2459272 2.369083 0.003837289 0.03376509 0.8552083 2.8505914
## [6,] 18.1395130 21.440689 5.309626102 3.96379280 4.8613534 7.5137649
## 534 543 550
## [1,] 0.000000000 0.0000000 0.0000000000
## [2,] 0.001838161 0.1237550 0.0000000000
## [3,] 0.017125918 0.2629671 0.0000000000
## [4,] 0.021997111 0.4067754 0.0000000000
## [5,] 0.010838236 0.1085778 0.0007849804
## [6,] 0.031090088 0.2624809 0.0270864312## [1] 13149 9# Nota:
# - exact_extract() calcula estadísticas del raster dentro de cada polígono
# - fun = "mean" → promedio areal de precipitación
# - El resultado se transpone (t) para tener:
# filas = tiempo (días)
# columnas = subcuencasOrganización en formato de serie temporales.
# Asociar los datos de precipitación con las fechas
dpr_cat <- zoo(pr_cat, dates)
# Ver primeras filas
head(dpr_cat)## 501 510 521 525 527 533
## 1981-01-01 0.0000000 0.000000 0.000000000 0.00000000 0.0000000 0.0000000
## 1981-01-02 3.6393113 4.073518 0.552590847 0.43440908 0.6454489 0.5282456
## 1981-01-03 4.2326713 4.979711 0.873725295 2.31732988 3.5731442 1.7681453
## 1981-01-04 5.4923182 9.815919 0.947537184 1.89248145 3.7520499 5.3061824
## 1981-01-05 0.2459272 2.369083 0.003837289 0.03376509 0.8552083 2.8505914
## 1981-01-06 18.1395130 21.440689 5.309626102 3.96379280 4.8613534 7.5137649
## 534 543 550
## 1981-01-01 0.000000000 0.0000000 0.0000000000
## 1981-01-02 0.001838161 0.1237550 0.0000000000
## 1981-01-03 0.017125918 0.2629671 0.0000000000
## 1981-01-04 0.021997111 0.4067754 0.0000000000
## 1981-01-05 0.010838236 0.1085778 0.0007849804
## 1981-01-06 0.031090088 0.2624809 0.0270864312# Visualizar series (todas las subcuencas)
plot(dpr_cat,
type = 'h',
col = 'dodgerblue',
main = "Precipitación diaria por subcuenca (mm/d)")
# Nota:
# zoo permite manejar series temporales con fechas reales
# filas = tiempo, columnas = subcuencas# Convertir objeto zoo a data.frame
df <- data.frame(
Date = index(dpr_cat), # Fechas de la serie temporal
coredata(dpr_cat), # Valores de precipitación por subcuenca
check.names = FALSE # Mantener nombres originales de subcuencas
)
# Convertir a formato largo (una fila por fecha y subcuenca)
df_long <- pivot_longer(
df,
cols = -Date, # Todas las columnas excepto la fecha
names_to = "subbasin", # Identificador de subcuenca
values_to = "precip" # Valores de precipitación
)
# Obtener valor máximo global para usar una misma escala en todos los paneles
ymax_global <- max(df_long$precip, na.rm = TRUE)
# Construir gráfico con ggplot2
plt <- ggplot(df_long, aes(x = Date, y = precip)) +
# Serie temporal como línea
geom_line(
color = "dodgerblue", # Color de la serie
linewidth = 0.35 # Línea delgada para series largas
) +
# Un panel por subcuenca
facet_wrap(
~ subbasin, # Variable de agrupación
scales = "fixed", # Misma escala en todos los paneles
ncol = 3 # Número de columnas en la grilla
) +
# Limitar eje Y usando el máximo global
coord_cartesian(
ylim = c(0, ymax_global)
) +
# Etiquetas del gráfico
labs(
title = "Precipitación diaria por subcuenca",
subtitle = "Series extraídas desde PISCO v2.1 (1981-2016)",
x = "Fecha",
y = "Precipitación (mm/día)"
) +
# Tema base
theme_bw() +
# Ajustes estéticos
theme(
strip.text = element_text(face = "bold", size = 9), # Títulos de cada subcuenca
axis.text.x = element_text(size = 7), # Texto eje X
axis.text.y = element_text(size = 7), # Texto eje Y
axis.title = element_text(size = 10), # Títulos de ejes
plot.title = element_text(size = 13, face = "bold"), # Título principal
plot.subtitle = element_text(size = 10), # Subtítulo
panel.grid.minor = element_blank(), # Quitar grilla menor
legend.position = "none" # Sin leyenda
)
# Mostrar gráfico
plt
Agregación de datos diarios a mensuales
# Convertir precipitación diaria a mensual (acumulado)
mpr_cat <- daily2monthly(dpr_cat, FUN = 'sum')
# Visualizar series mensuales
plot(mpr_cat,
type = 'h',
col = 'dodgerblue',
main = "Precipitación mensual acumulada (mm/mes)")
# Nota:
# FUN = "sum" se usa porque la precipitación es acumulativa
# (a diferencia de temperatura, donde se usaría "mean")Agregación de datos mensuales a anuales
# Convertir precipitación mensual a anual (acumulado)
ypr_cat <- monthly2annual(mpr_cat, FUN = 'sum')
# Visualizar series anuales
plot(ypr_cat,
type = 'h',
col = 'dodgerblue',
main = "Precipitación anual acumulada (mm/año)")
# Nota:
# Se mantiene FUN = "sum" porque la precipitación sigue siendo acumulativa
# El resultado representa el total de lluvia por año en cada subcuenca# Exportar datos en formato xlsx
write.xlsx(x = list(
Diario = data.frame(Fecha = index(dpr_cat), coredata(dpr_cat)),
Mensual = data.frame(Fecha = index(mpr_cat), coredata(mpr_cat)),
Anual = data.frame(Fecha = index(ypr_cat), coredata(ypr_cat))
),
rowNames = FALSE,
file = './Modules/Clase_2/outputs/precip_pisco_subcuencas.xlsx')4.1 Extracción por estación
# Definir coordenadas (longitud, latitud)
xy <- data.frame(
x = c(-75.5, -75.2),
y = c(-13.6, -13.1)
)
# Convertir a objeto espacial (puntos)
pts <- vect(xy,
geom = c("x", "y"),
crs = "EPSG:4326")
plot(cat)
plot(pts, col='red', add=TRUE)
# Extraer valores del raster en esos puntos
pr_pts <- t(terra::extract(pisco_cat, pts))
rownames(pr_pts) <- NULL
# Nota:
# - Cada punto obtiene una serie temporal completa del raster
# - El resultado se transpone para tener:
# filas = tiempo
# columnas = puntos# Asociar los datos de precipitación con las fechas
dpr_pts <- zoo(pr_pts, dates)
# Ver primeras filas
head(dpr_pts)##
## 1981-01-01 1.00000000 2.000000
## 1981-01-02 0.00000000 0.000000
## 1981-01-03 0.00000000 6.130610
## 1981-01-04 0.00000000 5.622951
## 1981-01-05 0.00000000 13.433528
## 1981-01-06 0.02360236 3.692428# Visualizar series (todas las subcuencas)
plot(dpr_pts,
type = 'h',
col = 'dodgerblue',
main = "Precipitación diaria por estación (mm/d)")
# Nota:
# zoo permite manejar series temporales con fechas reales
# filas = tiempo, columnas = subcuencas# Convertir precipitación diaria a mensual (acumulado)
mpr_pts <- daily2monthly(dpr_pts, FUN = 'sum')
# Visualizar series mensuales
plot(mpr_pts,
type = 'h',
col = 'dodgerblue',
main = "Precipitación mensual acumulada (mm/mes)")
# Nota:
# FUN = "sum" se usa porque la precipitación es acumulativa
# (a diferencia de temperatura, donde se usaría "mean")# Convertir precipitación mensual a anual (acumulado)
ypr_pts <- monthly2annual(mpr_pts, FUN = 'sum')
# Visualizar series anuales
plot(ypr_pts,
type = 'h',
col = 'dodgerblue',
main = "Precipitación anual acumulada (mm/año)")
# Nota:
# Se mantiene FUN = "sum" porque la precipitación sigue siendo acumulativa
# El resultado representa el total de lluvia por año en cada subcuenca# Exportar datos en formato xlsx
write.xlsx(x = list(
Diario = data.frame(Fecha = index(dpr_pts), coredata(dpr_pts)),
Mensual = data.frame(Fecha = index(mpr_pts), coredata(mpr_pts)),
Anual = data.frame(Fecha = index(ypr_pts), coredata(ypr_pts))
),
rowNames = FALSE,
file = './Modules/Clase_2/outputs/precip_pisco_estaciones.xlsx')