1. IntroducciÓn

La evapotranspiración es un componente clave para el estudio del ciclo hidrológico (Díaz et al., 2011), con su estimación se logra explicar el intercambio de agua y energía entre el suelo, la superficie terrestre y la atmósfera (Huang et al., 2021), es decir, que depende de múltiples parámetros para su estimación (Martínez, 2002). Medir la EvapotranspiraciÓn es difícil, la informaciÓn obtenida en campo genera algunas estimaciones puntuales para interpolar/extrapolar en un Área determinada, incrementando errores en la aplicaciÓn del mÉtodo (Ochoa-Sánchez et al., 2019).

En ecosistemas semiÁridos la evapotranspiraciÓn representa un gran porcentaje del presupuesto de agua, teniendo en cuenta que el balance hídrico en estos ecosistemas estÁ regulado por la radiaciÓn solar y eventos pluviomÉtricos, provocando una humedad en el suelo que restringe las interacciones suelo-planta-atmÓsfera (Olivera-Guerra et al., 2014). El departamento de La Guajira presenta condiciones excepcionales, tiene influencia de la costa caribe, así como la presencia de sistemas montanosos y ecosistemas semidesÉrticos, lo cual influye en el comportamiento del ciclo hidrolÓgico (CORPOGUAJIRA & ASOCARIBE, 2018). Tradicionalmente estÁ dividido en Alta, Media y Baja Guajira (CORPOGUAJIRA, 2009) presentando condiciones ambientales y zonas de vida diferentes. Un estudio desarrollado por Fedesarrollo, (2019) indica que: En el municipio de Uribia (Alta Guajira), las lluvias oscilan alrededor de los 300 mm/ano, originando sequía estacional durante gran parte del ano. Hacia el sur, en Manaure y Maicao (Media Guajira), las lluvias se incrementan ligeramente alcanzando en algunos sectores cantidades cercanas a los 1000 mm/año. En la mayor parte de Riohacha y hacia los límites con los departamento del Cesar, en el extremo sur, las lluvias pueden alcanzar los 1500 mm/año, a pesar de esto, la tasa de evaporaciÓn es muy alta, posicionÁndolo como uno de los departamentos con mÁs alto riesgo de sequía (IDEAM, 2012).

El presente trabajo tiene como objetivo estimar la evapotranspiración potencial a traves del metodo Hargreaves. (Hargreaves, 1983) mediante el uso de programaciÓn en RStudio y su aplicaciÓn en el departamento de La Guajira.

2. Área de estudio y datos

2.1. Área de estudio

El Área de estudio corresponde al departamento de La Guajira, ubicado al extremo nororiental de Colombia (Fig.1). El territorio Guajiro tiene una superficie de 20.848 km2 y cuenta con la mayor extensiÓn de costa marina sobre el Mar Caribe (GobernaciÓn de La Guajira, 2020), muestra una alta variabilidad espacial y temporal de sus recursos hídricos, debido a sus particulares condiciones climÁticas, geolÓgicas y morfolÓgicas, registrando los menores volúmenes de lluvias en el país (CORPOGUAJIRA, 2016; IDEAM, 2014).

El departamento de La Guajira tiene muchas particularidades en su territorio, alberga zonas muy cÁlidas con desiertos semiÁridos hasta sistemas montañosos como la Sierra Nevada de Santa Marta. Cuenta con una operaciÓn minera a cielo abierto de exportaciÓn de carbÓn ubicada en la zona baja que tiene influencia directa en el flujo de agua para el río de mayor importancia del departamento El Ranchería (Díaz et al., 2020). Sus ecosistemas predominantes son: desiertos, humedales en zonas montañosas y bosque tropical seco propenso a la aridez, baja vegetaciÓn y una fuerte estacionalidad de lluvias y adaptada a condiciones de estrÉs hídrico (GobernaciÓn de La Guajira, 2020)

Figura 1. UbicaciÓn del departamento de La Guajira

2.2. Datos

Se utilizaron datos de temperatura media, mínima y máxima; así como de precipitación total a resolución diaria, los datos fueron solicitados al IDEAM (Instituto de Hidrología, Meteorología y Estudios Ambientales) para un periodo de analisis desde el año 1990-2010. Se realizó un buffer de 10km sobre el departamento de La Guajira, obteniendo información cercana de los departamentos Magdalena y Cesar.

3. Enfoque metodolÓgico

Colombia ha aplicado y validado diferentes mÉtodos empíricos para estimar la evapotranspiraciÓn potencial en sus diferentes zonas pluviomÉtricas, dentro de los cuales se encuentran los siguientes:

Para la Zona pluviomÉtrica Caribe - Cesar (zona en donde se encuentra el departamento de La Guajira), el mÉtodo que presenta un mejor ajuste es el FAO Penman Monteith, seguido del Hargreaves (IDEAM, 2018).

Para este trabajo, se aplicÓ el mÉtodo Hargreaves para estimar la evapotranspiraciÓn potencial en el departamento de La Guajira, teniendo en cuenta que su ecuaciÓn de cÁlculo es mucho mÁs sencilla, veamos:

\(ETP= 0,0023 * (Tmed + 17,78)*Ro*(tmax - tmin)*0,5\)

Donde,

ETP= EvapotranspiraciÓn potencial Tmed= Temperatura media Ro= RadiaciÓn en el tope de la atmÓsfera Tmax= Temperatura mÁxima Tmin= Temperatura mínima

A continuaciÓn, se presenta el script utilizado para estimar la evapotranspiraciÓn potencial a partir de informaciÓn satelital.

3.1 Preparación de los datos

  • Los datos se organizaron por variables, en este caso, se crearon cuatro (4) archivos que contenian la información de las estaciones enviadas. Es importante mencionar que se seleccionaron solo aquellas estaciones que cumplieran la condicion de tener maximo 30% de datos faltantes.

  • Se creo un archivo .csv que contenia los códigos y coordenadas de las estaciones por cada variable

  • Del shape Departamentos (descargado desde el IGAC), se exportó unicamente el departamento de La Guajira

3.2 Interpolación

Para realizar la interpolación de los datos, se utilizó el procedimiento propuesto por Cristian Vargas el cual puede encotrarse en el siguiente link https://rpubs.com/CDVS/885607

  • Utilizamos la herramienta Fishnet desde ArcToolbox para generar la malla para el departamento de La Guajira (ver figura 2)

Figura 2. Malla para el departamento de La Guajira

INICIA EL CODIGO

Cargamos las librerias necesarias

 rm(list = ls()) #Limpiar el Environment
library(zoo)#Infraestructura S3 para series temporales regulares e irregulares
library(xts) #Series temporales
library(lubridate) #Manejo de fechas
library(terra) # Información Geográfica
library(sp)#Metodos espaciales 
library(raster) #Procesar información raster
library(phylin) #Metodo de interpolación IDW
library(rgdal) #Enlaces para la biblioteca de abstracción de datos 'Geoespacial'
library(spatial)
library(SPEI)
library(gtools)
library(future)
library(furrr)
#Directorio donde esta el archivo netcdf
setwd("D:/DOCUMENTOS/Parcial/")
#Temperatura mínima 
estacionesTMN <- vect("Descargas/XYestTMN_pl.shp") # Lee el shp de estacionesTMN
as.data.frame(estacionesTMN[,1:6]) # Muestra un extrato de la tabla de atributos
plot(estacionesTMN, 3) # Plotea la ubicación de las estacionesTMN

Datos_TMN=read.csv('datos/TMN_CON.csv',sep = ';')
Datos_TMN_Interpolados<-Datos_TMN
columnas<-length(Datos_TMN_Interpolados)
filas=nrow(Datos_TMN_Interpolados)

for (i in 1:1) {
  
  
  datos.sinteticosTMN <- as.numeric(as.vector(Datos_TMN_Interpolados[i,2:24])) # Datos sintéticos 
  
  stn <- geom(estacionesTMN)[,c("x", "y")] # Coordenadas de las estacionesTMN
  stn <- as.data.frame(stn)
  
  grid <- vect("Referencias/mallacorregidaOK.shp") # shp de la red de puntos
  grid <- geom(grid)[,c("x", "y")] # coordenadas de los puntos
  grid <- as.data.frame(grid)
  
  datos.interpoladosTMN <- idw(datos.sinteticosTMN, stn, grid, p = 1) # Aplicación IDW
  datos.interpoladosTMN<- round(datos.interpoladosTMN[,1],2)
  
  grid <- data.frame(grid, valores = datos.interpoladosTMN) 
  coordinates(grid) <- ~x+y # Se especifica que 'x' y 'y' son las coordenadas
  gridded(grid) <- T # Se indica que se trata de una red de puntos
  gridTMN <- raster(grid, "valores") # Se crea el raster
  projection(grid) <- crs("+init=epsg:3116") # Se le asigna un sistemas de coordenadas
  path='D:/DOCUMENTOS/Parcial/Interpolacion/TMN'
  numeracion=as.character(i)
  Nombre="Campo_TMN"
  Extension=".tif"
  writeRaster(gridTMN, paste(path,Nombre,numeracion,Extension), overwrite = T) # Se escribe el raster
  
}

plot(gridTMN, main = "TMN 01/01/1990")
#Temperatura maxima

estacionesTMX <- vect("D:/DOCUMENTOS/Parcial/Descargas/XYestTMX_pl.shp") # Lee el shp de estacionesTMX
as.data.frame(estacionesTMX[,1:6]) # Muestra un extrato de la tabla de atributos
plot(estacionesTMX, 3) # Plotea la ubicación de las estacionesTMX

Datos_TMX=read.csv('datos/TMX_CON.csv',sep = ';')
Datos_TMX_Interpolados<-Datos_TMX
columnas<-length(Datos_TMX_Interpolados)
filas=nrow(Datos_TMX_Interpolados)

for (i in 1:1) {
  
  
  datos.sinteticosTMX <- as.numeric(as.vector(Datos_TMX_Interpolados[i,2:18])) # Datos sintéticos 
  
  stn <- geom(estacionesTMX)[,c("x", "y")] # Coordenadas de las estacionesTMX
  stn <- as.data.frame(stn)
  
  grid <- vect("Referencias/mallacorregidaOK.shp") # shp de la red de puntos
  grid <- geom(grid)[,c("x", "y")] # coordenadas de los puntos
  grid <- as.data.frame(grid)
  
  datos.interpoladosTMX <- idw(datos.sinteticosTMX, stn, grid, p = 1) # Aplicación IDW
  datos.interpoladosTMX<- round(datos.interpoladosTMX[,1],2)
  
  grid <- data.frame(grid, valores = datos.interpoladosTMX) 
  coordinates(grid) <- ~x+y # Se especifica que 'x' y 'y' son las coordenadas
  gridded(grid) <- T # Se indica que se trata de una red de puntos
  gridTMX <- raster(grid, "valores") # Se crea el raster
  projection(grid) <- crs("+init=epsg:3116") # Se le asigna un sistemas de coordenadas
  path='D:/DOCUMENTOS/Parcial/Interpolacion/TMX'
  numeracion=as.character(i)
  Nombre="Campo_TMX"
  Extension=".tif"
  writeRaster(gridTMX, paste(path,Nombre,numeracion,Extension), overwrite = T) # Se escribe el raster
  
}

plot(gridTMX,main = "TMX 01/01/1990")
#Temperatura media

estacionesTSSM <- vect("D:/DOCUMENTOS/Parcial/Descargas/XYestTSSM_pl.shp") # Lee el shp de estacionesTSSM
as.data.frame(estacionesTSSM[,1:6]) # Muestra un extrato de la tabla de atributos
plot(estacionesTSSM, 3) # Plotea la ubicación de las estacionesTSSM

Datos_TSSM=read.csv('datos/TSSM.csv',sep = ';')
Datos_TSSM_Interpolados<-Datos_TSSM
columnas<-length(Datos_TSSM_Interpolados)
filas=nrow(Datos_TSSM_Interpolados)

for (i in 1:1) {
  
  datos.sinteticosTSSM <- as.numeric(as.vector(Datos_TSSM_Interpolados[1,2:6])) # Datos sintéticos 
  
  stn <- geom(estacionesTSSM)[,c("x", "y")] # Coordenadas de las estacionesTSSM
  stn <- as.data.frame(stn)
  
  grid <- vect("Referencias/mallacorregidaOK.shp") # shp de la red de puntos
  grid <- geom(grid)[,c("x", "y")] # coordenadas de los puntos
  grid <- as.data.frame(grid)
  
  datos.interpoladosTSSM <- idw(datos.sinteticosTSSM, stn, grid, p = 1) # Aplicación IDW
  datos.interpoladosTSSM<- round(datos.interpoladosTSSM[,1],2)
  
  grid <- data.frame(grid, valores = datos.interpoladosTSSM) 
  coordinates(grid) <- ~x+y # Se especifica que 'x' y 'y' son las coordenadas
  gridded(grid) <- T # Se indica que se trata de una red de puntos
  gridTSSM <- raster(grid, "valores") # Se crea el raster
  projection(grid) <- crs("+init=epsg:3116") # Se le asigna un sistemas de coordenadas
  path='D:/DOCUMENTOS/Parcial/Interpolacion/TSSM'
  numeracion=as.character(1)
  Nombre="Campo_TSSM"
  Extension=".tif"
  writeRaster(gridTSSM, paste(path,Nombre,numeracion,Extension), overwrite = T) # Se escribe el raster
  
}

plot(gridTSSM,main = "TSMM 01/01/1990")
#As raster 

stck <- addLayer(gridTMN, gridTMX, gridTSSM)

rasterToPoints(stck,spatial = FALSE)

#As a Table 


plan(cluster, workers=1,gc=TRUE)
etpt <- furrr:: future_map(.x=1:length(unique(stck[,3])), .f=1:(length(unique(stck[,3])))) {
       
  tbl <- stck[,3]
  lat<-unique(tbl)
  tbl<- (etp=as.numeric(hargreaves(Tmin=stck$valores.1, Tmax =stck$valores.2, lat= lat)))
  print("Done!")
 }                                             

plot(Evapotranspiración)