1 Introducción

Este documento calcula la evapotranspiración potencial diaria (ETP) con la ecuación de Hargreaves-Samani:

\[ETP = 0.0023 \cdot \frac{R_a}{2.45} \cdot (T_{med} + 17.78) \cdot (T_{max} - T_{min})^{0.5}\]

donde \(R_a\) es la radiación extraterrestre (MJ m⁻² día⁻¹) y las temperaturas están en °C. El procedimiento consiste en:

  1. Interpolar la temperatura (media, mínima y máxima) a partir de las estaciones, corrigiendo por elevación con un gradiente de 6.5 °C/km.
  2. Calcular la radiación extraterrestre para cada celda según su latitud y el día del año.
  3. Aplicar la ecuación de Hargreaves-Samani.

2 Preparación

2.1 Librerías y datos de entrada

rm(list = ls())
library(terra)

dem <- rast("DEM_ETP.tif")

gauges <- vect("ETP.gpkg", "Estaciones")

2.2 Elevación de las estaciones

Se extrae del DEM la elevación en cada estación.

gauges$elev <- extract(dem, gauges)[, 2]

3 Interpolación de la temperatura

3.1 Ejemplo con la temperatura media

Primero se lleva la temperatura de cada estación al nivel del mar, se interpola con IDW y luego se vuelve a la elevación de cada celda del DEM.

gauges$t_med_0 <- gauges$TSSM_MEDIA_D + 6.5 * gauges$elev / 1000

t_med_0 <- interpIDW(dem, gauges, "t_med_0", 1e6)

plot(t_med_0, main = "Temperatura media reducida al nivel del mar (°C)")

t_med <- t_med_0 - 6.5 * dem / 1000

plot(t_med, main = "Temperatura media (°C)")

3.2 Proceso en un ciclo para las tres variables

Se repite el procedimiento para la temperatura media, mínima y máxima.

fields <- c("TSSM_MEDIA_D", "TSSM_MN_D", "TSSM_MX_D")

temp <- list()
for (field in fields) {
  gauges[, paste0(field, "_0")] <- gauges[[field]] + 6.5 * gauges$elev / 1000
  temp[[field]] <- interpIDW(dem, gauges, paste0(field, "_0"), 1e6) - 6.5 * dem / 1000
}
par(mfrow = c(1, 3))
for (field in fields) plot(temp[[field]], main = field)

par(mfrow = c(1, 1))

4 Radiación extraterrestre

4.1 Raster de latitud

Se construye un raster con la latitud (en grados) de cada celda, proyectando de EPSG:9377 (MAGNA-SIRGAS / Origen Nacional) a EPSG:4686 (MAGNA-SIRGAS geográficas).

lat <- as.data.frame(dem, xy = TRUE) |>
  (\(x) {
    crds <- project(as.matrix(x)[, 1:2], from = "EPSG:9377", to = "EPSG:4686")
    cbind(x[, 1:2], crds[, 2])
  })() |>
  rast(crs = "EPSG:9377")

plot(lat, main = "Latitud (°)")

4.2 Cálculo de \(R_a\)

# Constante solar
Gr <- 0.082 # MJ*m^{-2}*min^{-1}

# Día del año
j <- as.Date("2026-10-27") |> lubridate::yday()

# Distancia relativa Tierra-Sol
dr <- 1 + 0.033 * cos(2 * pi / 365 * j)

# Declinación solar
delta <- 23.45 * pi / 180 * sin((360 / 365 * (284 + j)) * pi / 180)

# Ángulo horario al atardecer (omega s)
ws <- acos(-tan(lat * pi / 180) * tan(delta))

# Radiación solar extraterrestre
ra <- 24 * 60 / pi * Gr * dr *
  (ws * sin(lat * pi / 180) * sin(delta) +
     cos(lat * pi / 180) * cos(delta) * sin(ws))

plot(ra, main = "Radiación extraterrestre Ra (MJ m⁻² día⁻¹)")

5 ETP diaria

etp <- 0.0023 * (ra / 2.45) *
  (temp[["TSSM_MEDIA_D"]] + 17.78) *
  (temp[["TSSM_MX_D"]] - temp[["TSSM_MN_D"]])^0.5

plot(etp, main = "ETP diaria (mm/día)")