1 Detalle de la práctica

Este es un R Notebook donde se aplicará la ecuación de balance hídrico a largo plazo partiendo de información de sensores remotos. La entrada de precipitación corresponde al producto MSWEP y la de evapotranspiración potencial al producto GLEAM. El periodo seleccionado corresponde al comprendido entre 1980 y 2016.

La precipitación mensual se descargó siguiendo los pasos publicados en la página oficial http://www.gloh2o.org/mswep/ Este conjunto de datos tiene una resolución espacial de 0.1°X 0.1°

La evapotranspiración se descargó siguiendo los pasos publicados en la página oficial https://www.gleam.eu/ Este conjunto de datos tiene una resolución espacial de 0.25°X 0.25°

Desarrolle este ejercicio en un R Script

Los datos de entrada se encuentran en el siguiente link https://drive.google.com/file/d/1j7C2usHWJ1dXxAUdQ_0bmwvEpFQssBsG/view?usp=sharing

2 Configuraciones iniciales e información de entrada

Se cargan los paquetes que se van a utilizar para desarrollar el ejercicio y la ruta en la cual se encuentran los archivos.

rm(list = ls())
library(rgdal) #para manejar datos espaciales
library(raster) #para manejar datos espaciales
library(knitr) #para hacer tablas
library(ggplot2) #para hacer mapas

# actualizar la ruta según su caso personal
ruta_info <- "Z:/My Drive/PERSONALES/TRABAJOS/2022 CURSO HIDROLOGIA UNAL/3. CLASES/BALANCE HÍDRICO/EJERCICIO PRACTICO BH LARGO PLAZO"

Se carga el polígono de la cuenca que se va a evaluar, en este caso la cuenca del río Magdalena hasta la estación Pto Berrio Automática. Tenga en cuenta que la mayoría de productos satelitales y de reanálisis provienen en coordenadas geográficas, utilice todo en un mismo sistema coordenado para evitar problemas.

cuenca <- readOGR(dsn = paste0(ruta_info,"/SHP"), layer = "CUENCA[23097030]_EPSG4326")

## OGR data source with driver: ESRI Shapefile 
## Source: "Z:\My Drive\PERSONALES\TRABAJOS\2022 CURSO HIDROLOGIA UNAL\3. CLASES\BALANCE HÍDRICO\EJERCICIO PRACTICO BH LARGO PLAZO\SHP", layer: "CUENCA[23097030]_EPSG4326"
## with 1 features
## It has 18 fields
## Integer64 fields read as strings:  CODIGO_CAT

Al trabajar con datos espaciales es preferible utilizar un mismo sistema de coordenadas. Se procede a verificar que el sistema de la cuenca sea WGS84.

print(crs(cuenca)) #coordenadas geográficas

## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

Se calcula el área de la cuenca ya que será usada para calculos posteriores.

area_cuenca <- area(cuenca)/1000000 #unidades km2
print(paste0("El area de la cuenca es: ", round(area_cuenca,2), " km2"))

## [1] "El area de la cuenca es: 74402.1 km2"

plot(cuenca)

3 Precipitación

La información se encuentra disponible en formato .ncdf a nivel mensual. Por tanto, se debe agregar a escala anual para después calcular el valor total anual multianual.

ruta_p <- paste0(ruta_info,"/RASTER/P/MENSUAL/")
files_p <- list.files(ruta_p)

#discretizo cuales archivos corresponden a cada año
anos <- substring(files_p, 0,4)
anos_unicos <- unique(anos)

#En este ciclo se calcula la precipitación de cada año como la suma de cada uno de los meses
anual_p <- list()
for(a in 1:length(anos_unicos)){
  mensual_tmp <- list()
  #en este ciclo leo el archivo correspondiente a cada mes y lo recorto con la forma de la  cuenca
  for(m in 1:12){
    
    #cargar el raster del año a y mes m
    mensual_tmp[[m]] <- raster(paste0(ruta_info,"/RASTER/P/MENSUAL/", anos_unicos[a],
                                      sprintf("%02d", m),".nc"))
    
    #Cortar y enmascarar el raster con la forma de la cuenca
    mensual_tmp[[m]] <- crop(mensual_tmp[[m]], cuenca)
    mensual_tmp[[m]] <- mask(mensual_tmp[[m]], cuenca)
  }
  #sumo todos los recortes mensuales de cada año
  anual_p[[a]] <- sum(stack(mensual_tmp))
}

#promedio los campos anuales
anual_p <- mean(stack(anual_p))



plot(anual_p)

4 Evapotranspiración

La información se encuentra disponible en formato .ncdf a nivel anual. Por tanto, se debe calcular el valor total anual multianual.

ruta_et <-paste0(ruta_info,"/RASTER/ET/ANUAL/")

#cargar achivo que contiene la evapotranspiración anual 1980-2018
anual_et <- brick(paste0(ruta_info,"/RASTER/ET/ANUAL/E_1980_2018_GLEAM_v3.3a_YR.nc"))

#verificar la cantidad de capas contenidas en el ncdf (en este caso cada capa corresponde a un año)
anos_anual_et <- dim(anual_et)[3]
plot(anual_et[[1]])

Esta fuente de información se encuentra rotada, por tanto debo corregir esto.

#Necesito rotar los raster
anual_et <- t(flip(anual_et, direction='y' ))
plot(anual_et[[1]])

Necesito solo tomar los años entre 1980 y 2016.

anos_et_gleam <- 1980:2018
anos_et_requeridos <- 1980:2016

#obtengo las posiciones de los años que se requieren
anos_et_match <- unlist(lapply(anos_et_requeridos, function(x){which(x==anos_et_gleam)}))


anual_et <- anual_et[[anos_et_match]]

#calculo el valor medio de los rasters anuales y los recorto con el polígono de la cuenca
anual_et <- mean(anual_et)

anual_et <- crop(anual_et, cuenca)
anual_et <- mask(anual_et, cuenca)

plot(anual_et)

5 Caudal a partir del balance a largo plazo

Se aplica la ecuación de balance hídrico

#Con los ejercicios anteriores obtuvimos campos de precipitación y evapotranspiración así que calculo el valor
#medio de la cuenca
p_cuenca <- cellStats(anual_p, "mean")
et_cuenca <- cellStats(anual_et, "mean")

q_bh<- (p_cuenca - et_cuenca)*area_cuenca*0.001*(1/(365*24*60*60))*1000000

print(paste0("El caudal calculado a partir del balance a largo plazo es: " , round(q_bh,0), " m3/s"))

## [1] "El caudal calculado a partir del balance a largo plazo es: 1693 m3/s"

6 Caudal medido en la estación Pto Berrío Automática [23097030]

Utilizo la serie de caudales medios registrados en la estación.

q_estacion <- read.csv(paste0(ruta_info,"/Q_MED_M_23097030-PTO_BERRIO_AUTOMATICA.csv"),skip = 1)
q_estacion <- mean(q_estacion$Valor..Cubic.Metres.Per.Second.)

print(paste0("El caudal calculado con los registros de la estación es: " , round(q_estacion,0), " m3/s"))

## [1] "El caudal calculado con los registros de la estación es: 2299 m3/s"

tabla <- c("Fuente", "Caudal m3/s",
           "Estación", round(q_estacion,0),
           "BH Largo plazo",round(q_bh,0))
tabla <- matrix(tabla, nrow = 3, byrow = T)
colnames(tabla) <- tabla[1,]
tabla <- tabla[-1,]

kable(tabla, caption="Resultados")

Table 6.1: Resultados
Fuente	Caudal m3/s
Estación	2299
BH Largo plazo	1693

BALANCE HÍDRICO A LARGO PLAZO CON INFORMACIÓN DE SENSORES REMOTOS

Camila García-Echeverri