Inicio ← Clase 2 - Parte 1 Descargar .Rmd Clase 3 - Parte 1 →

Curso-Taller: Análisis Hidrológico de Cuencas con R

Clase 2: Procesamiento de Datos Climáticos Grillados (parte 2)

📦 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.

# 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)
  "hydroGOF",       # métricas de desempeño hidrológico
  "hyfo",           # corrección de sesgo (bias correction)
  "dplyr",          # manipulación de datos
  "tidyr",          # transformación de datos (formato largo)
  "ggplot2",        # visualización de datos
  "ggspatial",      # elementos cartográficos en mapas
  "tidyterra",      # integración de terra con ggplot2
  "ggnewscale",     # múltiples escalas en gráficos
  "openxlsx",       # lectura y escritura de archivos Excel
  "getPass"         # manejo de credenciales (opcional)
)

# 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. Lectura de base de datos de estaciones pluviométricas

# leer archivo de datos de precipitacion
raw  <- read.xlsx('./Data/OBS/data_precip.xlsx',
                  sheet=1, colNames=FALSE, detectDates = TRUE)

# previsualizar datos
head(raw, 10) 

##            X1                  X2                  X3                  X4
## 1         nom          Cusicancha             Hornada            Huancano
## 2         abr                 CUS                 HOR                 HUA
## 3          id              113086              151509              113009
## 4         lon -75.299989999999994 -75.099999999999994 -75.616690000000006
## 5         lat           -13.47791           -13.42793           -13.59456
## 6         alt                3263                4530                1005
## 7  1981-01-01                   0                   0                   0
## 8  1981-01-02                   0                   0                   0
## 9  1981-01-03                   0                   0                   0
## 10 1981-01-04                   2                   0                   0
##                     X5        X6                  X7
## 1             Letrayoc   Ticrapo              Totora
## 2                  LET       TIC                 TOT
## 3               113153    156117              156119
## 4  -75.744979999999998 -75.43329 -75.316689999999994
## 5  -13.652240000000001 -13.37795           -13.11145
## 6                  756      2174                3966
## 7               -999.9         0                   0
## 8               -999.9         0                   0
## 9               -999.9         0                   0
## 10              -999.9         0                   0

# Extraer metadatos de estaciones
abr  <- as.character(raw[2,-1])
print(abr)

## [1] "CUS" "HOR" "HUA" "LET" "TIC" "TOT"

cod  <- as.character(raw[3,-1])
print(cod)

## [1] "113086" "151509" "113009" "113153" "156117" "156119"

lon  <- as.numeric(raw[4,-1])
print(lon)

## [1] -75.29999 -75.10000 -75.61669 -75.74498 -75.43329 -75.31669

lat  <- as.numeric(raw[5,-1])
print(lat)

## [1] -13.47791 -13.42793 -13.59456 -13.65224 -13.37795 -13.11145

npts <- length(abr)
print(npts)

## [1] 6

Mostrar mapa de estaciones pluviométricas en la cuenca.

# crear objeto espacial de estaciones
xy     <- data.frame(x=lon, y=lat)
pts    <- vect(xy, geom=c("x","y"), crs="EPSG:4326")
pts_sf <- st_as_sf(pts)
pts_sf$nombre <- abr

# cargar objetos espacialesde la cuenca
cat <- vect('./Data/SIG/subbasins.shp')
riv <- vect('./Data/SIG/rivers.shp')
dem <- rast('./Data/SIG/dem.tif')

# generar raster de relieve
slp <- terra::terrain(dem, "slope", unit="radians")
asp <- terra::terrain(dem, "aspect", unit="radians")
hsd <- terra::shade(slp, asp, angle=45, direction=300, normalize=TRUE)
dem.df <- as.data.frame(dem, xy=TRUE)
names(dem.df)[3] <- "alt"
hsd.df <- as.data.frame(hsd, xy=TRUE)

# crear mapa
plt <- ggplot() +
  geom_raster(data = hsd.df,
              aes(x, y, fill = hillshade),
              show.legend = FALSE) +
  scale_fill_distiller(palette = "Greys") +
  new_scale_fill() +
  scale_x_continuous(n.breaks = 4) +
  geom_raster(data = dem.df,
              aes(x, y, fill = alt),
              alpha = 0.7) +
  scale_fill_gradientn(
    colours = terrain.colors(100),
    name = "msnm"
  ) +
  geom_sf(data = st_as_sf(cat), color = 'black', fill = NA) +
  geom_sf(data = st_as_sf(riv), color = 'dodgerblue', fill = NA) +
  geom_sf(data = pts_sf, color = 'red', size = 2) +
  geom_text(data = pts_sf,
            aes(x = st_coordinates(pts_sf)[,1],
                y = st_coordinates(pts_sf)[,2],
                label = abr),
            color = "black", size = 3, hjust = -0.2) +   # ← aquí faltaba +

  theme_bw() +
  
  guides(
    fill = guide_colorbar(
      barheight = unit(6, "cm"),
      barwidth  = unit(0.4, "cm")
    )) +
  
  theme(legend.position = "right") + 
  
  annotation_north_arrow(
    location = "tl",
    which_north = "true",
    pad_x = unit(0.1, "in"),
    pad_y = unit(0.1, "in"),
    style = north_arrow_fancy_orienteering
  ) +
  
  annotation_scale(location = "br", width_hint = 0.5) +
  
  labs(
    title = 'Cuenca de estudio',
    subtitle = 'Mapa de estaciones pluviométricas',
    x = 'Longitud',
    y = 'Latitud',
    fill = 'msnm'
  )

# mostrar mapa
print(plt)

Ahora extraemos los datos de precipitación diaria del archivo Excel de la base de datos.

# Extraer fechas
time <- as.Date(raw[-1:-6,1])
head(time)

## [1] "1981-01-01" "1981-01-02" "1981-01-03" "1981-01-04" "1981-01-05"
## [6] "1981-01-06"

# Extraer datos (notar que datos faltantes se muestran como -999.9)
prec <- data.frame(sapply(raw[-1:-6,-1], 'as.numeric'))
head(prec)

##   X2  X3 X4     X5  X6 X7
## 1  0 0.0  0 -999.9 0.0  0
## 2  0 0.0  0 -999.9 0.0  0
## 3  0 0.0  0 -999.9 0.0  0
## 4  2 0.0  0 -999.9 0.0  0
## 5  4 1.8  0 -999.9 0.0  0
## 6  0 6.6  0 -999.9 0.5  5

# Reconocer datos faltantes y crear series de tiempo
prec[prec==-999.9] <- NA

# Transformar datos a formato de serie de tiempo
colnames(prec) <- abr
prec_obs       <- zoo(prec, time)
head(prec_obs)

##            CUS HOR HUA LET TIC TOT
## 1981-01-01   0 0.0   0  NA 0.0   0
## 1981-01-02   0 0.0   0  NA 0.0   0
## 1981-01-03   0 0.0   0  NA 0.0   0
## 1981-01-04   2 0.0   0  NA 0.0   0
## 1981-01-05   4 1.8   0  NA 0.0   0
## 1981-01-06   0 6.6   0  NA 0.5   5

Mostramos las series de tiempo de precipitación diaria observada en estaciones.

# crear dataframe con datos
df <- data.frame(
  Date = index(prec_obs),
  coredata(prec_obs),
  check.names = FALSE
)

# transformar a formato largo
df_long <- pivot_longer(
  df,
  cols = -Date,
  names_to = "station",
  values_to = "precip"
)

# eje Y global
ymax_global <- max(df_long$precip, na.rm = TRUE)

# crear figura de serie de tiempo
plt <- ggplot(df_long, aes(Date, precip)) +
        geom_line(color = "dodgerblue", linewidth = 0.4) +
        facet_wrap(~ station, scales = "fixed", ncol = 3) +
        coord_cartesian(ylim = c(0, ymax_global)) +
        labs(
          title = "Series de precipitación diaria (1981-2016)",
          subtitle = "Estaciones pluviométricas en la cuenca",
          x = "Tiempo",
          y = "(mm/día)"
        ) +
        theme_bw() +
        theme(
          strip.text = element_text(face = "bold", size = 9),
          axis.text.x = element_text(size = 7),
          axis.text.y = element_text(size = 7),
          axis.title = element_text(size = 10),
          plot.title = element_text(size = 13, face = "bold"),
          plot.subtitle = element_text(size = 10),
          panel.grid.minor = element_blank(),
          legend.position = "none"
        )

# mostrar figura
plt

## Warning: Removed 35829 rows containing missing values or values outside the scale range
## (`geom_line()`).

Se exploran los estadísticos descriptivos de las series observadas, incluyendo medidas de tendencia central, dispersión y extremos, con el fin de caracterizar su comportamiento y variabilidad.

# estadisticos descriptivos
smry(prec_obs)

##               Index        CUS        HOR        HUA        LET        TIC
## Min.     1981-01-01     0.0000     0.0000     0.0000     0.0000     0.0000
## 1st Qu.  1990-01-01     0.0000     0.0000     0.0000     0.0000     0.0000
## Median   1999-01-01     0.0000     0.0000     0.0000     0.0000     0.0000
## Mean     1999-01-01     1.0480     1.9470     0.0346     0.0258     0.7294
## 3rd Qu.  2008-01-01     0.0000     2.4000     0.0000     0.0000     0.0000
## Max.     2016-12-31    54.5000    37.9000    12.9000    13.4000    33.0000
## IQR            <NA>     0.0000     2.4000     0.0000     0.0000     0.0000
## sd             <NA>     2.9133     3.7814     0.3938     0.3172     2.4976
## cv             <NA>     2.7793     1.9425    11.3701    12.2751     3.4241
## Skewness       <NA>     4.2379     2.8947    17.7792    22.1931     5.5866
## Kurtosis       <NA>    28.2207    11.5668   412.6224   678.2928    43.1577
## NA's           <NA>  1334.0000  7914.0000  5079.0000  7024.0000 10656.0000
## n              <NA> 13149.0000 13149.0000 13149.0000 13149.0000 13149.0000
##                 TOT
## Min.         0.0000
## 1st Qu.      0.0000
## Median       0.0000
## Mean         1.6240
## 3rd Qu.      0.0000
## Max.        37.0000
## IQR          0.0000
## sd           3.8107
## cv           2.3463
## Skewness     3.2932
## Kurtosis    13.7098
## NA's     10350.0000
## n        13149.0000

Se revisa la disponibilidad de los datos observados, verificando la presencia de vacíos, la continuidad temporal y la calidad general de los registros.

# revisar porcentaje de disponibilidad
disp <- colSums(!is.na(prec_obs)) / nrow(prec_obs) * 100
print(disp)

##      CUS      HOR      HUA      LET      TIC      TOT 
## 89.85474 39.81291 61.37349 46.58149 18.95962 21.28679

# mostrar histograma de disp. de datos
hist(disp,
     col ='dodgerblue',
     main = "Disponibilidad de datos de precip. (1981-2016)",
     xlab = "Porcentaje de disponibilidad (%)",
     ylab = 'Número de estaciones')

# revisar datos completos por año
dwi(prec_obs)

##      CUS HOR HUA LET TIC TOT
## 1981 365 365 365   0 365 365
## 1982 365 365  59   0 365 365
## 1983 365 304   0   0 365 365
## 1984 366 274   0   0 366 366
## 1985 365 365   0   0 365 335
## 1986 273 365   0   0 334 365
## 1987 274 335 276   0 151 364
## 1988 183 366 274   0 182 274
## 1989 245 365  61   0   0   0
## 1990 124 365 153   0   0   0
## 1991 365 365   0   0   0   0
## 1992 216 274   0   0   0   0
## 1993 365 306  92   0   0   0
## 1994 365 365 276   0   0   0
## 1995 365 365 365   0   0   0
## 1996 366  91   0   0   0   0
## 1997 365   0   0   0   0   0
## 1998 365   0   0   0   0   0
## 1999 365   0 245  35   0   0
## 2000 366   0 366 335   0   0
## 2001 365   0 365 365   0   0
## 2002 335   0 365 365   0   0
## 2003 365   0 365 365   0   0
## 2004 366   0 366 366   0   0
## 2005 365   0 365 365   0   0
## 2006 365   0 365 365   0   0
## 2007 365   0 365 365   0   0
## 2008 366   0 366 366   0   0
## 2009 365   0 365 365   0   0
## 2010 365   0 365 365   0   0
## 2011 365   0 365 365   0   0
## 2012 366   0 366 308   0   0
## 2013 365   0 365 365   0   0
## 2014 365   0 365 365   0   0
## 2015 273   0 365 365   0   0
## 2016  31   0  60 335   0   0

Se acotan los datos al periodo de análisis 2000–2001, en función de la disponibilidad de información de CHIRPS (archivos descargados), con el objetivo de realizar una comparación consistente del sesgo entre las series.

# extraer datos observados para el periodo de analisis
prec_obs <- window(prec_obs, start='2000-01-01', end='2001-12-31')

head(prec_obs)

##             CUS HOR HUA LET TIC TOT
## 2000-01-01  5.4  NA 0.0 0.0  NA  NA
## 2000-01-02  4.4  NA 0.0 0.0  NA  NA
## 2000-01-03  9.0  NA 0.5 0.0  NA  NA
## 2000-01-04  9.7  NA 0.0 0.8  NA  NA
## 2000-01-05 14.9  NA 0.6 0.2  NA  NA
## 2000-01-06 13.8  NA 0.0 0.0  NA  NA

☁️ 3. Lectura de datos grillados CHIRPS (v2 y v3)

3.1 Datos CHIRPS v2

Se leen los archivos en formato NetCDF descargados desde https://www.chc.ucsb.edu/data/chirps. Cabe señalar que la versión CHIRPS v2 será descontinuada en diciembre de 2026, por lo que su uso debe considerarse con fines comparativos o de transición hacia versiones más recientes. Para la descarga y manejo de los datos no se utiliza el paquete de R chirps, debido a limitaciones de acceso a los servidores, optándose en su lugar por métodos directos de descarga.

# lectura archivos netCDF chirps v2
files <- list.files('./Data/CHIRPS/examples/v2',
                    pattern = '*.nc',
                    full.names = TRUE)
head(files)

## [1] "./Data/CHIRPS/examples/v2/chirps-v2.0.2000.days_p05.nc"
## [2] "./Data/CHIRPS/examples/v2/chirps-v2.0.2001.days_p05.nc"

# leer archivos netcdf como raster multicapas
dates_chirps <- seq(as.Date('2000-01-01'),
                    as.Date('2001-12-31'),
                    by='days')
r <- rast(files)

# mostrar metadatos del raster
print(r)

## class       : SpatRaster 
## size        : 2000, 7200, 731  (nrow, ncol, nlyr)
## resolution  : 0.05, 0.05  (x, y)
## extent      : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
## sources     : chirps-v2.0.2000.days_p05.nc  (366 layers) 
##               chirps-v2.0.2001.days_p05.nc  (365 layers) 
## varnames    : precip (Climate Hazards group InfraRed Precipitation with Stations) 
##               precip (Climate Hazards group InfraRed Precipitation with Stations) 
## names       : precip_1, precip_2, precip_3, precip_4, precip_5, precip_6,    ... 
## unit        : mm/day 
## time (days) : 2000-01-01 to 2001-12-31 (731 steps)

# mostrar como ejemplo el mapa raster para la capa 100
plot(r[[100]], main = "CHIRPS v2 - capa 100")

# extraer series de tiempo para puntos de estaciones
prec_chirps_v2 <- zoo(t(terra::extract(r, pts, ID = FALSE)), dates_chirps)
colnames(prec_chirps_v2) <- abr
head(prec_chirps_v2)

##                  CUS      HOR      HUA      LET       TIC       TOT
## 2000-01-01  7.903306 9.299727 3.639511 2.631507 10.367119 13.213758
## 2000-01-02 15.806612 6.199819 7.279023 5.263014 20.734238 26.427515
## 2000-01-03  3.951653 9.299727 0.000000 0.000000  5.183559  0.000000
## 2000-01-04  0.000000 3.099909 0.000000 0.000000  0.000000  0.000000
## 2000-01-05  0.000000 0.000000 0.000000 0.000000  0.000000  6.606879
## 2000-01-06 11.109358 9.686875 0.000000 3.888637 11.012827  0.000000

3.2 Datos CHIRPS v3 (mezclados con IMERG)

Existen dos variantes de CHIRPS v3 para la generación del producto diario: una versión “sat”, que utiliza información satelital (como IMERG) para desagregar la precipitación pentadal, y una versión “rnl”, que emplea reanálisis (ERA5) con el mismo propósito. Ambas parten de una base pentadal común, pero difieren en la fuente utilizada para representar la variabilidad temporal diaria. Los datos se encuentran disponibles en: https://www.chc.ucsb.edu/data/chirps3.

# lectura de datos chirps v3
files <- list.files('./Data/CHIRPS/examples/v3',
                    pattern = '*.tif$',
                    full.names = TRUE)
head(files)

## [1] "./Data/CHIRPS/examples/v3/chirps-v3.0.sat.2000.01.01.tif"
## [2] "./Data/CHIRPS/examples/v3/chirps-v3.0.sat.2000.01.02.tif"
## [3] "./Data/CHIRPS/examples/v3/chirps-v3.0.sat.2000.01.03.tif"
## [4] "./Data/CHIRPS/examples/v3/chirps-v3.0.sat.2000.01.04.tif"
## [5] "./Data/CHIRPS/examples/v3/chirps-v3.0.sat.2000.01.05.tif"
## [6] "./Data/CHIRPS/examples/v3/chirps-v3.0.sat.2000.01.06.tif"

# leer archivos netcdf como raster multicapas
dates_chirps <- seq(as.Date('2000-01-01'),
                    as.Date('2001-12-31'),
                    by='days')
r <- rast(files)

# mostrar metadatos del archivo
print(r)

## class       : SpatRaster 
## size        : 2400, 7200, 731  (nrow, ncol, nlyr)
## resolution  : 0.05, 0.05  (x, y)
## extent      : -180, 180, -60, 60  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : chirps-v3.0.sat.2000.01.01.tif  
##               chirps-v3.0.sat.2000.01.02.tif  
##               chirps-v3.0.sat.2000.01.03.tif  
##               ... and 728 more sources
## names       : chirp~01.01, chirp~01.02, chirp~01.03, chirp~01.04, chirp~01.05, chirp~01.06, ...

# mostrar mapa raster para la capa 100 como ejemplo
x <- r[[100]]
x[x == -9999] <- NA  # el tif tiene pixeles NA como -9999
plot(x, main = "CHIRPS v3 - capa 100")

# extraer series de tiempo para puntos de estaciones
prec_chirps_v3 <- zoo(t(terra::extract(r, pts, ID = FALSE)), dates_chirps)
colnames(prec_chirps_v3) <- abr
head(prec_chirps_v3)

##                  CUS       HOR       HUA       LET       TIC       TOT
## 2000-01-01  3.991274 16.360737 0.5604521 0.2443996  5.239185 21.232407
## 2000-01-02 13.962648 10.810702 3.1621339 3.7474608 16.277941 25.442059
## 2000-01-03  7.097113  5.755259 3.7528841 0.0000000 14.763000  4.879518
## 2000-01-04  3.623477 20.531244 1.4632002 0.0000000  8.152246  2.976121
## 2000-01-05  1.607407  4.772933 0.3137329 0.1629331  4.795808 17.604809
## 2000-01-06  3.636910  7.423131 0.0000000 3.5026150  2.704137  1.378378

🏞️ 4. Comparativo de datos observados vs CHIRPS (sin corrección)

Evaluar comparativamente la precipitación diaria observada en estaciones pluviométricas frente a los datos de CHIRPS en sus distintas versiones, con el fin de analizar su capacidad para reproducir la variabilidad temporal, los extremos y el sesgo respecto a las observaciones.

# crear dataframe con datos observados y chirps
df <- data.frame(
  Date = index(prec_obs),
  Obs = coredata(prec_obs),
  CHIRPS_v2 = coredata(prec_chirps_v2),
  CHIRPS_v3 = coredata(prec_chirps_v3),
  check.names = FALSE
)

# transformar a formato largo
df_long <- df %>%
  pivot_longer(
    cols = -Date,
    names_to = c("type", "station"),
    names_sep = "\\.",
    values_to = "precip"
  ) %>%
  mutate(
    type = recode(
      type,
      "Obs" = "Observado",
      "CHIRPS_v2" = "CHIRPS v2",
      "CHIRPS_v3" = "CHIRPS v3"
    )
  ) %>%
  filter(!is.na(precip))

# ordenar estaciones para mostrar
station_order <- colnames(prec_obs)
df_long$station <- factor(df_long$station, levels = station_order)

# crear figura comparativa de series de tiempo
plt <- ggplot(df_long, aes(x = Date, y = precip,
                           color = type, linetype = type)) +
  geom_line(linewidth = 0.35) +
  scale_color_manual(values = c(
    "Observado" = "black",
    "CHIRPS v2" = "dodgerblue",
    "CHIRPS v3" = "red"
  )) +
  scale_linetype_manual(values = c(
    "Observado" = "solid",
    "CHIRPS v2" = "solid",
    "CHIRPS v3" = "solid"
  )) +
facet_wrap(~ station, ncol = 3, scales = "free_y") + 
  labs(
    title = "Comparación de precipitación diaria por estación",
    subtitle = "Observado vs CHIRPS v2 y CHIRPS v3",
    x = "Fecha",
    y = "Precipitación (mm/día)"
  ) +
  theme_bw() +
  theme(
    strip.text = element_text(face = "bold", size = 9),
    axis.text.x = element_text(size = 6),
    axis.text.y = element_text(size = 7),
    axis.title = element_text(size = 11),
    plot.title = element_text(size = 14, face = "bold"),
    plot.subtitle = element_text(size = 10),
    panel.grid.minor = element_blank(),
    legend.position = "top",
    legend.title = element_blank()
  )

# mostrar figura
plt

Se calculan las métricas de desempeño KGE, PBIAS y nRMSE para cada conjunto de datos, con el fin de evaluar la capacidad de las series simuladas para reproducir las observaciones en términos de ajuste global, sesgo y error relativo.

# Convertir a matrices (mismas dimensiones)
obs  <- coredata(prec_obs)
v2   <- coredata(prec_chirps_v2)
v3   <- coredata(prec_chirps_v3)

stations <- colnames(obs)

# Inicializar resultados
gof_raw <- data.frame(
  station = stations,
  
  KGE_v2 = NA, PBIAS_v2 = NA, nRMSE_v2 = NA,
  KGE_v3 = NA, PBIAS_v3 = NA, nRMSE_v3 = NA
)

# Calcular metricas de desempeño (KGE, PBIAS y nRMSE)
for (i in seq_along(stations)) {
  o  <- obs[, i]
  s2 <- v2[, i]
  s3 <- v3[, i]
  
  # Remover NA conjuntos
  idx2 <- complete.cases(o, s2)
  idx3 <- complete.cases(o, s3)
  
  # CHIRPS v2
  gof_raw$KGE_v2[i]   <- suppressWarnings(KGE(sim = s2[idx2], obs = o[idx2]))
  gof_raw$PBIAS_v2[i] <- suppressWarnings(pbias(sim = s2[idx2], obs = o[idx2]))
  gof_raw$nRMSE_v2[i] <- suppressWarnings(nrmse(sim = s2[idx2], obs = o[idx2], norm = "sd"))
  
  # CHIRPS v3
  gof_raw$KGE_v3[i]   <- suppressWarnings(KGE(sim = s3[idx3], obs = o[idx3]))
  gof_raw$PBIAS_v3[i] <- suppressWarnings(pbias(sim = s3[idx3], obs = o[idx3]))
  gof_raw$nRMSE_v3[i] <- suppressWarnings(nrmse(sim = s3[idx3], obs = o[idx3], norm = "sd"))
}

# mostrar tabla de resultados
print(gof_raw)

##   station      KGE_v2 PBIAS_v2 nRMSE_v2     KGE_v3 PBIAS_v3 nRMSE_v3
## 1     CUS   0.5370129      5.6     94.5   0.548565     -8.8     85.6
## 2     HOR          NA       NA       NA         NA       NA       NA
## 3     HUA -15.6481042   1623.7    489.2  -3.486001    433.7    194.2
## 4     LET -25.3000026   2548.2    785.7 -11.810734   1233.0    442.1
## 5     TIC          NA       NA       NA         NA       NA       NA
## 6     TOT          NA       NA       NA         NA       NA       NA

# transformar resultados de desempeño a formato largo
gof_long <- gof_raw %>%
  pivot_longer(
    cols = -station,
    names_to = c("metric", "product"),
    names_pattern = "(.*)_(v2|v3)",
    values_to = "value"
  ) %>%
  mutate(
    product = recode(product,
                     "v2" = "CHIRPS v2",
                     "v3" = "CHIRPS v3"),
    metric = recode(metric,
                    "PBIAS" = "PBIAS (%)",
                    "nRMSE" = "nRMSE (%)")
  )

# crear figura de desempeño
plt <- ggplot(gof_long, aes(x = station, y = value, fill = product)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  facet_wrap(~ metric, scales = "free_y", ncol = 3) +
  labs(
    title = "Métricas de desempeño por estación",
    subtitle = "Comparación entre CHIRPS v2 y CHIRPS v3",
    x = "Estación",
    y = "Valor de la métrica",
    fill = NULL
  ) +
  theme_bw() +
  theme(
    strip.text = element_text(face = "bold", size = 10),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 7),
    axis.text.y = element_text(size = 8),
    axis.title = element_text(size = 11),
    plot.title = element_text(size = 14, face = "bold"),
    plot.subtitle = element_text(size = 10),
    panel.grid.minor = element_blank(),
    legend.position = "top"
  )

# mostrar figura
plt

## Warning: Removed 18 rows containing missing values or values outside the scale range
## (`geom_col()`).

🏞️ 5. Corrección del sesgo de datos CHIRPS respecto de observaciones

Se aplica el método de mapeo de cuantiles empíricos (EQM) para la corrección del sesgo en los datos CHIRPS. Este método ajusta la distribución acumulada empírica de los datos simulados a la distribución de referencia observada, corrigiendo no solo el valor medio, sino también la variabilidad y la frecuencia de ocurrencia de los eventos. Al tratarse de una corrección distribucional, permite mejorar la representación de toda la gama de precipitaciones, incluyendo eventos extremos, garantizando una mejor correspondencia estadística entre las series modeladas y observadas.

# Función para corregir CHIRPS usando EQM
correct_eqm <- function(prec_sim, prec_obs, prThreshold = 0.1) {
  
  tim <- index(prec_obs)
  obs <- coredata(prec_obs)
  sim <- coredata(prec_sim)
  
  stations <- colnames(prec_obs)
  
  sim_bc <- sim
  stations_no_bc <- c()
  
  for (i in seq_along(stations)) {
    
    obs_i <- obs[, i]
    sim_i <- sim[, i]
    
    idx <- complete.cases(obs_i, sim_i)
    
    # Si no hay suficientes datos observados, no se corrige
    if (sum(idx) < 30) {
      stations_no_bc <- c(stations_no_bc, stations[i])
      next
    }
    
    obs_df <- data.frame(date = tim[idx], val = obs_i[idx])
    sim_df <- data.frame(date = tim[idx], val = sim_i[idx])
    
    # Se corrige solo el periodo con datos coincidentes
    sim_bc[idx, i] <- biasCorrect(
      frc = sim_df,
      hindcast = sim_df,
      obs = obs_df,
      method = "eqm",
      preci = TRUE,
      prThreshold = prThreshold,
      extrapolate = "no"
    )[, 2]
  }
  
  sim_bc <- zoo(sim_bc, index(prec_sim))
  colnames(sim_bc) <- colnames(prec_sim)
  
  return(list(
    bc = sim_bc,
    stations_no_bc = stations_no_bc
  ))
}

# Aplicar corrección a CHIRPS v2
res_v2 <- correct_eqm(
  prec_sim = prec_chirps_v2,
  prec_obs = prec_obs
)

prec_chirps_v2_bc <- res_v2$bc

head(prec_chirps_v2_bc)

##                  CUS      HOR      HUA LET       TIC       TOT
## 2000-01-01  8.200000 9.299727 0.000000 0.0 10.367119 13.213758
## 2000-01-02 18.447005 6.199819 2.142857 0.8 20.734238 26.427515
## 2000-01-03  3.876959 9.299727 0.000000 0.0  5.183559  0.000000
## 2000-01-04  0.000000 3.099909 0.000000 0.0  0.000000  0.000000
## 2000-01-05  0.000000 0.000000 0.000000 0.0  0.000000  6.606879
## 2000-01-06 11.000000 9.686875 0.000000 0.0 11.012827  0.000000

# Aplicar corrección a CHIRPS v3
res_v3 <- correct_eqm(
  prec_sim = prec_chirps_v3,
  prec_obs = prec_obs
)

prec_chirps_v3_bc <- res_v3$bc

head(prec_chirps_v2_bc)

##                  CUS      HOR      HUA LET       TIC       TOT
## 2000-01-01  8.200000 9.299727 0.000000 0.0 10.367119 13.213758
## 2000-01-02 18.447005 6.199819 2.142857 0.8 20.734238 26.427515
## 2000-01-03  3.876959 9.299727 0.000000 0.0  5.183559  0.000000
## 2000-01-04  0.000000 3.099909 0.000000 0.0  0.000000  0.000000
## 2000-01-05  0.000000 0.000000 0.000000 0.0  0.000000  6.606879
## 2000-01-06 11.000000 9.686875 0.000000 0.0 11.012827  0.000000

# creamos un dataframe con los datos
df <- data.frame(
  Date = index(prec_obs),
  Obs = coredata(prec_obs),
  CHIRPS_v2 = coredata(prec_chirps_v2_bc),
  CHIRPS_v3 = coredata(prec_chirps_v3_bc),
  check.names = FALSE
)

# transformamos el dataframe en formato largo
df_long <- df %>%
  pivot_longer(
    cols = -Date,
    names_to = c("type", "station"),
    names_sep = "\\.",
    values_to = "precip"
  ) %>%
  mutate(
    type = recode(
      type,
      "Obs" = "Observado",
      "CHIRPS_v2" = "CHIRPS v2",
      "CHIRPS_v3" = "CHIRPS v3"
    )
  ) %>%
  filter(!is.na(precip))

# ordenamos en funcion de estaciones
station_order <- colnames(prec_obs)
df_long$station <- factor(df_long$station, levels = station_order)

# crear figura comparativa
plt <- ggplot(df_long, aes(x = Date, y = precip,
                           color = type, linetype = type)) +
  geom_line(linewidth = 0.35) +
  scale_color_manual(values = c(
    "Observado" = "black",
    "CHIRPS v2" = "dodgerblue",
    "CHIRPS v3" = "red"
  )) +
  scale_linetype_manual(values = c(
    "Observado" = "solid",
    "CHIRPS v2" = "solid",
    "CHIRPS v3" = "solid"
  )) +
facet_wrap(~ station, ncol = 3, scales = "free_y") + 
  labs(
    title = "Comparación de precipitación diaria por estación",
    subtitle = "Observado vs CHIRPS v2 (bc) y CHIRPS v3 (BC)",
    x = "Fecha",
    y = "Precipitación (mm/día)"
  ) +
  theme_bw() +
  theme(
    strip.text = element_text(face = "bold", size = 9),
    axis.text.x = element_text(size = 6),
    axis.text.y = element_text(size = 7),
    axis.title = element_text(size = 11),
    plot.title = element_text(size = 14, face = "bold"),
    plot.subtitle = element_text(size = 10),
    panel.grid.minor = element_blank(),
    legend.position = "top",
    legend.title = element_blank()
  )

#mostrar figura comparativa
plt

Se recalculan las métricas de desempeño para las series corregidas, con el fin de evaluar la mejora en la representación de la variabilidad, el sesgo y el ajuste respecto a las observaciones.

# Convertir a matrices (mismas dimensiones)
obs  <- coredata(prec_obs)
v2   <- coredata(prec_chirps_v2_bc)
v3   <- coredata(prec_chirps_v3_bc)

stations <- colnames(obs)

# Inicializar resultados
gof_bc <- data.frame(
  station = stations,
  
  KGE_v2 = NA, PBIAS_v2 = NA, nRMSE_v2 = NA,
  KGE_v3 = NA, PBIAS_v3 = NA, nRMSE_v3 = NA
)

# Calcular metricas de desempeño
for (i in seq_along(stations)) {
  o  <- obs[, i]
  s2 <- v2[, i]
  s3 <- v3[, i]
  
  # Remover NA conjuntos
  idx2 <- complete.cases(o, s2)
  idx3 <- complete.cases(o, s3)
  
  # CHIRPS v2
  gof_bc$KGE_v2[i]   <- suppressWarnings(KGE(sim = s2[idx2], obs = o[idx2]))
  gof_bc$PBIAS_v2[i] <- suppressWarnings(pbias(sim = s2[idx2], obs = o[idx2]))
  gof_bc$nRMSE_v2[i] <- suppressWarnings(nrmse(sim = s2[idx2], obs = o[idx2], norm = "sd"))
  
  # CHIRPS v3
  gof_bc$KGE_v3[i]   <- suppressWarnings(KGE(sim = s3[idx3], obs = o[idx3]))
  gof_bc$PBIAS_v3[i] <- suppressWarnings(pbias(sim = s3[idx3], obs = o[idx3]))
  gof_bc$nRMSE_v3[i] <- suppressWarnings(nrmse(sim = s3[idx3], obs = o[idx3], norm = "sd"))
}

# mostrar tabla de resultados
print(gof_bc)

##   station      KGE_v2 PBIAS_v2 nRMSE_v2     KGE_v3 PBIAS_v3 nRMSE_v3
## 1     CUS  0.53662553      1.9     96.7 0.58617959      1.0     91.1
## 2     HOR          NA       NA       NA         NA       NA       NA
## 3     HUA  0.03822173     -6.0    137.4 0.01824547      3.0    141.2
## 4     LET -0.01082602     -4.6    141.9 0.13763188     -2.8    131.1
## 5     TIC          NA       NA       NA         NA       NA       NA
## 6     TOT          NA       NA       NA         NA       NA       NA

# transformar a formato largo
gof_long <- gof_bc %>%
  pivot_longer(
    cols = -station,
    names_to = c("metric", "product"),
    names_pattern = "(.*)_(v2|v3)",
    values_to = "value"
  ) %>%
  mutate(
    product = recode(product,
                     "v2" = "CHIRPS v2",
                     "v3" = "CHIRPS v3"),
    metric = recode(metric,
                    "PBIAS" = "PBIAS (%)",
                    "nRMSE" = "nRMSE (%)")
  )

# crear figura comparativa del desempeño
plt <- ggplot(gof_long, aes(x = station, y = value, fill = product)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.7) +
  facet_wrap(~ metric, scales = "free_y", ncol = 3) +
  labs(
    title = "Métricas de desempeño por estación",
    subtitle = "Comparación entre CHIRPS v2 (bc) y CHIRPS v3 (bc)",
    x = "Estación",
    y = "Valor de la métrica",
    fill = NULL
  ) +
  theme_bw() +
  theme(
    strip.text = element_text(face = "bold", size = 10),
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 7),
    axis.text.y = element_text(size = 8),
    axis.title = element_text(size = 11),
    plot.title = element_text(size = 14, face = "bold"),
    plot.subtitle = element_text(size = 10),
    panel.grid.minor = element_blank(),
    legend.position = "top"
  )

# mostrar figura
plt

## Warning: Removed 18 rows containing missing values or values outside the scale range
## (`geom_col()`).

Se exportan las series de datos CHIRPS v2 y v3 corregidas por sesgo en formato Excel para cada estación, facilitando su análisis, visualización y uso en aplicaciones posteriores.

# Exportar datos
write.xlsx(list(
    Observado = data.frame(Fecha=index(prec_obs), coredata(prec_obs)),
    CHIRPS_v2_bc = data.frame(Fecha=index(prec_chirps_v2_bc), coredata(prec_chirps_v2_bc)),
    CHIRPS_v3_bc = data.frame(Fecha=index(prec_chirps_v3_bc), coredata(prec_chirps_v3_bc))
    ),
    './Modules/Clase_2/outputs/precip_estaciones_chirps_bc.xlsx',
    rowNames= FALSE)