Implementación de una Métrica de Evaluación de Similitud entre Raster
Diego Cortés Ramos - diacortesra@unal.edu.co
2022-10-19
En el presente cuaderno se expone el proyecto final de la asignatura SIG en el Manejo de Recursos Hídricos de la Maestría en Ingeniería - Recursos Hidráulicos de la Universidad Nacional de Colombia.
1 Introducción
En el marco del proyecto de tesis titulado “Evaluación de la Estimación Espaciotemporal de la Recarga Mediante un Modelo Hidrológico Utilizando una Calibración Multiobjetivo que Incorpore Información de Sensores Remotos” es necesario la implementación de métricas que permitan establecer la similitud entre dos raster para poder llevar a cabo el proceso de calibración propuesto.
Para la comparación entre series hidrológicas existen múltiples métricas propuestas como el Nash-Sutcliffe Efficiency (NSE), error cuadrático medio (RMSE), Kling-Gupta Model Efficiency (KGE), sin embargo, no es posible evaluar la similitud espacial de series espaciales mediante estas métricas. Debido al auge del uso de sensores remotos y el aumento de la capacidad computacional de los ordenadores, se ha venido involucrando la información de sensores remotos en modelación hidrológica y autores como Dembélé et. al. (2020), Koch et. al. (2015, 2018) han empezado a proponer diversas métricas para evaluar el desempeño espacial de los modelos.
Mediante el presente proyecto se va realizar la implementaición de una métrica de comparación, así como la aplicación otras estrategias para analizar la similitud de dos series temporales de rasters de evapotranspiración potencial (ETP) en una cuenca del Valle Medio del Magdalena, la primera es estimada mediante el método de Hargreaves-Samani (George H. Hargreaves & Zohrab A. Samani, 1985) a partir de registros en superficie del IDEAM; la segunda es descargada del satélite MODIS.
En el cuaderno se tratan temas como descarga automática de datos de MODIS, interpolación de datos de temperatura, álgebra de mapas para el cálculo de la ETP, entre otros; los cuales pueden ser de consulta para cualquier persona que lo esté leyendo.
rm(list = ls())
# setwd("~/Diego/03.Academico/MAESTRIA/SIG RH/PROYECTO/")
wd <- "D:/Documentos/Diego/03.Academico/MAESTRIA/SIG RH/PROYECTO/"2 Área de Estudio
El río Lebrija es un afluente del río Magdalena en su cuenca media. Su desembocadura marca la división entre los departamentos de Cesar y Santander, y forma un aplio complejo cenagoso. A la altura del corregimiento del San Rafael de Lebrija, en Rionegro, Santander, se delimitó la cuenca tomando como referencia la estación hidrológica San Rafael (23197370). Esta se muestra en la Figura 2.1, así como los principales drenajes y la estación mencionada.
library(dplyr) # Funciones auxiliares
library(terra) # Datos geográficos
library(tidyterra) # Plot terra
library(raster) # Datos geográficos
library(leaflet) # Visor html de mapas
library(sf) # Datos geográficoscnca <- vect(paste0(wd, "GIS/Cuenca.shp")) # Lee la cuenca
drnjs <- vect(paste0(wd, "GIS/Drenajes.shp")) # Drenajes
stn.SR <- vect(paste0(wd, "GIS/EstacionCierre.shp")) # Estación San Rafael
xtn <- ext(project(cnca, "epsg:4326")) # Estensión de la cuenca en WGS84
# Solo para visualizar
cncaSpatial <- as(project(cnca, "epsg:4326"), "Spatial")
drnjsSpatial <- as(project(drnjs, "epsg:4326"), "Spatial")
stnSRSpatial <- as(project(stn.SR, "epsg:4326"), "Spatial")
map <- leaflet() %>% addTiles() %>%
setView(lat = mean(xtn[3:4]), lng = mean(xtn[1:2]), zoom = 9) %>%
addPolygons(data = cncaSpatial, color = "#830000", fillOpacity = 0.1, popup = "Cuenca") %>%
addPolygons(data = drnjsSpatial, color = "#104E8B", fillOpacity = 0.8, weight = 1.2,
popup = drnjs$NOMBRE_GEO) %>%
addCircleMarkers(data = stnSRSpatial, color = "#551A8B", opacity = 1, popup = stn.SR$nombre,
radius = 3)
mapFigure 2.1: Localización de la cuenca de estudio.
3 Información de Entrada
La información para el cálculo de la evapotranspiración provino de dos fuentes: de superficie mediante estaciones del IDEAM y de sensoramiento remoto del satélite MODIS. Adicionalmente, se utilizó un DEM de la zona de estudio en resolución de 30m proveniente de la misión ALOS PALSAR.
3.1 Información de Superficie
En el área de estudio y alrededores se recopilaron series de temperatura media, mínima y máxima de seis estaciones meteorológicas del IDEAM. Cabe destacar que no todas las estaciones cuentan con los registros de las tres variables, el caso más crítico corresponde con la temperatura máxima que solo cuenta con tres estaciones. Las estaciones seleccionadas se muestran en la Figura 3.1.
# Estaciones de temperatura
stn.temp <- vect(paste0(wd, "GIS/Estaciones.shp"))
# Para visualizar
stnTempSpatial <- as(project(stn.temp, "epsg:4326"), "Spatial")
mapstn <- addMarkers(map, data = stnTempSpatial, popup = stn.temp$nombre)
mapstnFigure 3.1: Estaciones con registros de temperatura
Los registros de temperatura coresponden al periodo 2000 a 2020 (ver Figura 3.2), previamente se les había realizado un análisis de datos anómalos. Sin embargo, debido a la poca cantidad de estaciones, para garantizar que en todos los días todas las estaciones contaran con registros para la interpolación se realizó un rellenado de datos con las librerías ‘RMAWGEN’ Y ‘RGENERATE’. En Anexo 8.1 se muestra la comparación de la caracterización estadística antes y después del llenado para verificar que este procedimiento no haya afectado considerablemente la serie; además en la Figura 3.3 muestra, a manera de ejemplo, la serie rellenada de temperatura media de la estación Aeropuerto Palonegro (23195502).
library(openxlsx) # Para leer y manejar archivos xlsx
library(ggplot2) # Para graficar
library(ggpubr) # Para graficar
library(reshape2) # Reorganizar data frame
library(RMAWGEN) # Ajustar modelo autoregresivo
library(RGENERATE) # Rellenar serie
# Hojas en archivo de excel con temperatura mínima, máxima y media, respectivamente.
shts <- c("TMN_CON", "TMX_CON", "TSSM_D")
temp <- tempF <- list() # Para guardar las temperaturas
plt <- list() # Para guardar plots
clrs <- c("#EEDC82", "#8B2500", "#CD853F") # colores por variable
names(clrs) <- shts
est0 <- est1 <- array(dim = c(5, nrow(stn.temp), 3),
dimnames = list(c("Media", "Mínima", "Máxima", "Mediana", "Desviación Estandar"),
stn.temp$CODIGO, shts)) # Estadsiticas
for(i in shts){
# Datos de temperatura de la hoja i
data <- read.xlsx(paste0(wd, "ARCHIVOS DE TRABAJO/Temperatura.xlsx"), i)
# Códigos de las estaciones en data
cod <- substring(colnames(data[,-1]), 2)
# Estadísticos en serie cruda
est0[,cod,i] <- apply(data[,-1], 2, function(x){
md <- mean(x, na.rm = T)
mn <- min(x, na.rm = T)
mx <- max(x, na.rm = T)
m50 <- median(x, na.rm = T)
sde <- sd(x, na.rm = T)
return(c(md, mn, mx, m50, sde))
})
# Para plot
dta.plot <- (t((data[, -1] + 1) / (data[, -1] + 1)) * 1:(ncol(data[, -1]))) %>% t
dta.plot <- cbind(as.Date(data[,1]), dta.plot) %>% as.data.frame()
dta.plot <- melt(dta.plot, id.vars = "V1", variable.name = "Temp")
dta.plot[,1] <- as.Date(dta.plot[,1], origin = "1970/01/01")
tlt <- switch (i,
TMN_CON = "Temperatura Mínima",
TMX_CON = "Temperatura Máxima",
TSSM_D = "Temperatura Media"
)
# Gráfica de registros
plt[[i]] <- ggplot(dta.plot, aes(V1, value)) + theme_bw() +
geom_point(color = clrs[i]) +
scale_y_continuous(breaks = 1:(ncol(data[, -1])), labels = cod) +
labs(x = "Fechas", y = NULL, title = tlt) +
theme(plot.title = element_text(hjust = 0.5))
# Rellenado
# Ajusta un modelo autoregresivo multisitio
md <- getVARmodel(data = data[,-1], suffix = NULL, n_GPCA_iteration=5,
n_GPCA_iteration_residuals=5)
# Aplica el modelo para rellenar la serie
dataFill <- generate(md, n = nrow(data), gap.filling = data[,-1], names = names(data)[-1])
# Estadísticos en serie rellenada
est1[,cod,i] <- apply(dataFill, 2, function(x){
md <- mean(x, na.rm = T)
mn <- min(x, na.rm = T)
mx <- max(x, na.rm = T)
m50 <- median(x, na.rm = T)
sde <- sd(x, na.rm = T)
return(c(md, mn, mx, m50, sde))
})
# Guarda la serie rellenada y la serie cruda
dataFill <- cbind(data[,1], dataFill)
temp[[i]] <- data
tempF[[i]] <- dataFill
}
# Formato estadísticas antes del rellenado
est0 <- melt(est0)
colnames(est0) <- c("Estadítica", "Estación", "Variable", "Valor")
est0 <- est0[!is.na(est0$Valor),]
# Formato estadísticas después del rellenado
est1 <- melt(est1)
colnames(est1) <- c("Estadítica", "Estación", "Variable", "Valor")
est1 <- est1[!is.na(est1$Valor),]
# Estadísticas resumen
est <- cbind(est0, est1$Valor)
est$Diferencia <- apply(est[,4:5], 1, function(x){
r <- x[1] - x[2]
})
colnames(est)[4:5] <- c("Original", "Rellenada")# Gráfica de registros
p <- ggarrange(plt[[1]], plt[[3]], plt[[2]], nrow = 3, ncol = 1, heights = c(4, 4, 3))
p
Figure 3.2: Registros de temperatura recopilados
# Gráfica ejemplo de serie rellenada
# data frame para ggplot
df <- cbind(tempF[["TSSM_D"]][,c("data[, 1]", "X23195502")], temp[["TSSM_D"]][,"X23195502"])
colnames(df) <- c("Fecha", "Cruda", "Rellenada")
df <- melt(df, id.vars = "Fecha", variable.name = "Serie")
df$Fecha <- as.Date(df$Fecha)
# Gráfica
p <- ggplot(df, aes(Fecha, value, color = Serie)) + geom_line() +
scale_color_manual(values = c("#1874CD", "#ADD8E6"), labels = c("Rellenada", "Original")) +
labs(x = NULL, y = "Temperatura (ºC)",
title = "Serie de Temperatura Media Diaria Rellenada\nEstación Aeropuerto Palonegro (23195502)") +
theme(plot.title = element_text(hjust = 0.5), legend.position = "bottom")
p
Figure 3.3: Ejemplo de serie rellenada
3.2 Información de sensores remotos
El satélite Moderate Resolution Imaging Spectroradiometer (MODIS) opera desde principios de sigo y cuenta con un producto de evapotranspiración (Running et al., 2019) aplicando un método basado en la ecuación del Penman-Monteith.
MODIS utiliza una modificación del algoritmo desarrollado por Cleugh et. al. (2007) el cual requiere de información de cobertura de la tierra y albedo, ambos calculados con la información satelital. Así mismo, utiliza información de temperatura, humedad y radiación proveniente de datos de superficie. El producto de ET se ofrece en resolución temporal de 8 días, mensual o anual y en resolución temporal de celdas de 500m.
Usando el paquete MODISstp (Busetto & Ranghetti, 2016) se descargaron datos el 2000 al 2022. Posteriormente se recortaron y proyectaron al sistema de coordenadas MAGNA/SRIGAS Origen Único Nacional (EPSG:9377). Así mismo, se le retiraron las celdas clasificadas como nubes, zonas urbanas y zonas cubiertas por agua; y se realizó una identificación de datos atípicos mediante boxplopt. Cuando en un raster, producto de la remoción de valores inválidos, se hubiera removido mas del 30% de los datos se descartó.
# Revisa si ya no hay datos descargados previamente
# NOTA: Larga duración de ejecución
if(!file.exists(paste0(wd, "/ET_MODIS/Net_ET_8Day_500m_v6/Time_Series/RData/Terra/PET_500m/MOD16A2_PET_500m_1_2001_1_2022_RData.RData"))){
library(MODIStsp)
# Crea directorio para guardar la información de MODIS
outfldr <- paste0(wd, "ET_MODIS")
dir.create(outfldr)
# Descarga la información del satélite MODIS
# Este código fue editado para ocultar el usuario y contraseña del autor
# remplace los **** por su propio usuario y contraseña.
# out_folder: directorio para guardar la información descargada
# out_folder_mod: directorio para guardar los archivos originales hdf
# selprod: producto de MODIS a descargar
# sensor: zona de interés en terrestre u oceanica
# bandsel: banda especifica dentro del producto
MODIStsp(gui = F,
out_folder = "ET_MODIS", out_folder_mod = "ET_MODIS",
user = "*****", password = "*****",
selprod = "Net_ET_8Day_500m (M*D16A2)",
sensor = "Terra", bandsel = "PET_500m",
start_date = "2000.01.01", end_date = "2022.01.01",
bbox = c(-8895610,0,-7783640,1111940),
spatmeth = 'bbox',
out_format = "GTiff", compress = "LZW",
delete_hdf = T, parallel = T
)
}
# Revisa si ya no hay datos procesados previamente
# NOTA: Larga duración de ejecución
if (!file.exists(paste0(wd, "ET_MODIS/ETP2319/ETP2319.RData"))){
# Los resultados se descargan como tif y en un RData por lo que se pueden usar
# Carga el stack descargado
load(paste0(wd, "/ET_MODIS/Net_ET_8Day_500m_v6/Time_Series/RData/Terra/PET_500m/MOD16A2_PET_500m_1_2001_1_2022_RData.RData"))
grd <- rast(raster_ts[[1]]) %>% project("epsg:9377") %>%
crop(cnca, mask = T) %>% crds(df = T)
# Proyecta a EPSG 9377 y recorta a la cuenca de interés
etpM <- rast(raster_ts) %>% project("epsg:9377") %>% crop(cnca, mask = T)
# Remueve atípicos de los rasters y descarga si se remueven mas del 50%
etpM2319 <- rast()
for(i in 1:dim(etpM)[3]){
ri <- etpM[[i]]
# Si no hay valores validos sigue al siguiente
n1 <- sum(!is.na(values(ri)))
if(n1 == 0) next
# valores > 32700 son para ciertas coberturas donde el algoritmo no es aplicable
ri[ri > 32700] <- NA
# Factor de escala de MODIS
ri <- 0.1 * ri
# Remueve iterativamente los atípicos mediante boxplots
while(T){
otl <- boxplot.stats(na.omit(values(ri)))$out
if(length(otl) == 0) break
ri[ri %in% otl] <- NA
}
# Valores validos después de la remoción
n2 <- sum(!is.na(values(ri)))
# Si se remueven mas del 30% se omite
if(n2/n1 < 0.7) next
# Lo guarda como stack
etpM2319 <- c(etpM2319, ri)
# Guarda el tif como respaldo
writeRaster(ri, paste0(wd, "ET_MODIS/ETP2319/", names(etpM)[i], ".tif"), overwrite = T)
}
# No se puede guardar un SpatRaster en un RData
etpM2319.Stack <- stack(etpM2319)
save(etpM2319.Stack, grd, file = paste0(wd, "ET_MODIS/ETP2319/ETP2319.RData"))
}else{
# Carga los datos procesados previamente
load(paste0(wd, "ET_MODIS/ETP2319/ETP2319.RData"))
etpM2319 <- rast(etpM2319.Stack)
}
yrs <- names(etpM2319) %>% substr(18,21) %>% as.numeric
etpM2319 <- etpM2319[[yrs <= 2020]]4 Métodos
En este punto ya se cuenta con la información espacial de ETP descargada del satélite MODIS, pero aún no se tiene la ETP con la información de superficie. Para tal fin se va a utilizar el método de Hargreaves - Samani (George H. Hargreaves & Zohrab A. Samani, 1985) mediante álgebra de mapas de los datos interpolados de temperatura. Con esto se van a aplicar dos técnicas de comparación, inicialmente se van a aplicar técnicas de comparación temática (Foody, 2002) y posteriormente se va a implementar la Métrica de Eficiencia de los Patrones Espaciales (\(E_{SP}\)).
4.1 Aplicación de Hargreaves - Samani para cálculo de la ETP
Para realizar la interpolación espacial de las temperaturas se va a utilizar el Inverso Ponderado de la Distancia (IDW, por sus siglas en ingles) de orden 1. Debido a la fuerte relación entre la temperatura y la elevación, se va a incluir el gradiente térmico de -6.5 ºC por cada mil metros de elevación. El procedimiento consiste en llevar todas las mediciones a una misma elevación (0 msnm) realizar la interpolación en esta condición y luego regresa las mediciones a la elevación real. Para facilitar las comparaciones se va a interpolar sobre la misma malla de coordenadas de la información descargada de MODIS.
# Revisa si ya hay valores interpolados
if(!all(file.exists(paste0(wd, "GIS/", names(tempF), ".tif")))){
library(phylin)
# Coordenadas las celdas de MODIS
crdsM <- grd
# DEM del área de estudio
dem <- rast(paste0(wd, "GIS/DEM.tif"))
# Elevaciones de las estaciones de acuerdo al DEM
stn.temp$elev <- extract(dem, stn.temp)[,1]
# Elevaciones de las celdas de MODIS
crdsM$elev <- extract(dem, crdsM)[,1]
# Interpolaciones para cada variable de temperatura
tidw <- list()
for(i in 1:length(tempF)){
# Datos de cada variable de temperatura
ti <- tempF[[i]]
cod <- colnames(ti)[-1] %>% substring(2)
# Estaciones y coordenadas
sti <- stn.temp[stn.temp$CODIGO %in% cod,]
sti.crds <- geom(sti, df = T)[,c("x", "y")]
# Asegura que el orden de las columnas coincida con el de las estaciones
ti <- ti[,paste0("X", sti$CODIGO)]
# Temperatura en cota 0
t0 <- ti + matrix(6.5 * sti$elev / 1000, ncol = ncol(ti), nrow = nrow(ti), byrow = T)
# Interpolación
tinterp <- apply(t0, 1, function(x, st, grd){
r <- idw(x, st, grd[,1:2], p = 1)[,1] - 6.5 * grd[,3] /1000
r <- cbind(grd[,1:2], r)
r <- rast(r, type = "xyz", crs = "epsg:9377")
}, st = sti.crds, grd = crdsM)
names(tinterp) <- tempF[[i]][,1]
tidw[[i]] <- tinterp
rm(tinterp)
}
tidw <- lapply(tidw, function(x) rast(x))
# Escribe in tif multicapa por cada variable
for(i in 1:length(tidw)){
writeRaster(tidw[[i]], paste0(wd, "GIS/", names(tempF)[i], ".tif"), overwrite = T)
}
}Posteriormente, se realiza el álgebra de mapas para calcular la ETP. Para esto se siguió el documento de Goméz-Blanco & Cadena (Gómez-Blanco & Cadena, 2018):
\[ETP=0.0023*R_{o}*\left(T_{med}+17.78\right)*\left(T_{max}-T_{min}\right)^{0.5}\]
Donde, \(T\) es la temperatura, en grados celsius, media (\(med\)), mínima (\(min\)) o máxima (\(max\)), según corresponda. Y \(R_{o}\) es la radiación en el tope de la atmósfera. Esta se expresa en milímetros de agua evaporada por día (\(mm/día\)), y se calcula con \(R_{o}=R_{a}/2.45\), donde \(R_{a}\) es la radiación extraterrestre, en \(MJ m^{-2}día^{-1}\). Típicamente, esta se obtiene de tablas debido a su complejidad en el cálculo; pero al usar un lenguaje de programación es incluso más sencillo programar la formula original que hacer la búsqueda en tablas.
\[R_{a}=\frac{24*60}{\pi}G_{sc}d_{r}\left[\omega_{s}\sin\varphi\sin\delta+\cos\varphi\cos\delta\sin\omega_{s}\right]\]
Donde, \(G_{sc}\) es una constante solar equivalente a \(0.082MJM^{-2}min^{-1}\), \(d_{r}\) es la distancia relativa al sol, \(\omega_{s}\) es el ángulo de radiación a la hora de puesta del sol, \(\delta\) es la declinación solar, y \(\varphi\) es la latitud, estos últimos tres en radianes. Y se calculan con las siguientes expresiones:
\[d_{r}=1+0.03\cos\left(\frac{2\pi}{365}J\right)\]
\[\delta=0.409\sin\left(\frac{2\pi}{365}J-1.39\right)\]
\[\omega_{s}=\cos^{-1}\left(-\tan\varphi\tan\delta\right)\]
Donde, \(J\) es el día del año, de 1 a 365. En resumen, el método de Hargreaves-Samani depende de las temperaturas, el día del año y la latitud.
# Indice de fechas donde hay información satelital
ind <- names(etpM2319) %>% substring(nchar(names(etpM2319)) - 7)
ind <- as.Date((substr(ind, 6, 8) %>% as.numeric) - 1,
origin = paste0(substr(ind, 1, 4), "-01-01"))
ind <- sapply(ind, function(x) seq(x - 7, x, by = "day"))
ind <- as.vector(ind) %>% as.Date(origin = "1970-01-01") %>% as.character
if(!all(file.exists(paste0(wd, "GIS/ETPS2319/", ind, ".tif")))){
library(parallel)
# Raster de latitudes
if(!file.exists(paste0(wd, "GIS/lat.tif"))){
lat <- grd
lat$z <- crds(project(tidw[[1]][[1]], "epsg:4686"), df = T)$y
lat <- rast(lat, type = "xyz", crs = "epsg:9377")
writeRaster(lat, paste0(wd, "GIS/lat.tif"))
}
t <- Sys.time()
cl <- makeCluster(8)
invisible(
parLapply(cl, ind, function(x, wd){
if(file.exists(paste0(wd, "GIS/ETPS2319/", x, ".tif"))) return(NULL)
library(dplyr)
library(lubridate)
library(terra)
tmn <- rast(paste0(wd, "GIS/TMN_CON.tif"))
tmx <- rast(paste0(wd, "GIS/TMX_CON.tif"))
tssm <- rast(paste0(wd, "GIS/TSSM_D.tif"))
lat <- rast(paste0(wd, "GIS/lat.tif"))
# Constante solar
Gr <- 0.082 #MJ*m^{-2}*min^{-1}
# Día del año
j <- as.Date(names(tmn)) %>% yday
names(j) <- names(tmn)
# dr
dr <- 1+0.033*cos(2*pi/365*j[x])
# delta
delta <- 23.45*pi/180*sin((360/365*(284+j[x]))*pi/180)
# Omega s
ws <- acos(-tan(lat*pi/180)*tan(delta))
# Radiación solar Extra terrestre
ra <- 24*60/pi*Gr*dr*(ws*sin(lat*pi/180)*sin(delta)+cos(lat*pi/180)*cos(delta)*sin(ws))
# ETP día x
eti <- 0.0023 * (ra / 2.45) * (tssm[[x]] + 17.78) * (tmx[[x]] - tmn[[x]]) ^ 0.5
writeRaster(eti, paste0(wd, "GIS/ETPS2319/", x, ".tif"), overwrite = T)
}, wd = wd)
)
stopCluster(cl)
}
# Agrega 8 días para hacerlo equivalente a MODIS
# Revisa si ya hay proceso previo
if(!file.exists(paste0(wd, "GIS/ETPS2319/0_ETPS2319.tif"))){
# Fechas para cada raster de MODIS
ind <- names(etpM2319) %>% substring(nchar(names(etpM2319)) - 7)
ind <- as.Date((substr(ind, 6, 8) %>% as.numeric) - 1,
origin = paste0(substr(ind, 1, 4), "-01-01"))
# Agrega las fechas
etpS2319 <- list()
for(i in 1:nlyr(etpM2319)){
indi <- seq(ind[i] - 7, ind[i], by = "day")
lyrsi <- paste0(wd, "GIS/ETPS2319/", indi, ".tif") %>% rast
etpS2319[[i]] <- sum(lyrsi)
}
etpS2319 <- rast(etpS2319)
names(etpS2319) <- ind
writeRaster(etpS2319, paste0(wd, "GIS/ETPS2319/0_ETPS2319.tif"))
}else{
# Lee si existe proceso previo
etpS2319 <- rast(paste0(wd, "GIS/ETPS2319/0_ETPS2319.tif"))
}4.2 Matriz de confusión
Para la evaluación de la simulutud de las series de rasters se va a hacer un muestreo aleatorio de mil puntos en la cuenca.
Imagen 4.1. Matriz de confusión. Fuente: Foody
La matriz de confusión (Imagen 4.1) es eje de la evaluación de exactitud de mapas temáticos (Foody, 2002). Consiste en tabular la categoría de las predicciones contra la categoría de las observaciones. A partir de este es posible calcular el porcentaje de exactitud global (\(I_{G}\)) o el indice kappa (\(\kappa\)) de acuerdo con las siguientes expresiones:
\[I_{G}=\frac{\sum_{n=1}^{q}n_{kk}}{n}*100\]
\[\kappa=\frac{\sum_{n=1}^{q}n_{kk}-\sum_{n=1}^{q}n_{k+}n_{+k}}{n^{2}-\sum_{n=1}^{q}n_{k+}n_{+k}}\]
Donde, \(\sum_{n=1}^{q}n_{kk}\) es la suma de la diagonal de la matrix, y \(\sum_{n=1}^{q}n_{k+}n_{+k}\) es el producto de la fila \(+k\) con la columna \(k+\), es decir la fila y columna de sumatoria.
# Seis intervalos
n <- 6
# Intervalos uniformes
brks <- seq(min(values(etpS2319), na.rm = T), max(values(etpS2319), na.rm= T),
length.out = n+1) %>% round(1)
# MODIS reclasificado
etpMR2319 <- etpM2319
etpMR2319[!is.na(etpMR2319)] <- 6
# ETP superficie reclasificada
etpSR2319 <- etpS2319
etpSR2319[!is.na(etpSR2319)] <- 6
# Reclasificación
for (i in (n-1):1){
etpMR2319[etpM2319 < brks[i+1]] <- i
etpSR2319[etpS2319 < brks[i+1]] <- i
}
# Puntos aleatorios para muestreo
set.seed(77170)
spl <- spatSample(cnca, 1000)
# Extraer valores de los raster de ETP
spl.val <- cbind(extract(etpM2319, spl) %>% as.matrix %>% as.vector,
extract(etpS2319, spl) %>% as.matrix %>% as.vector) %>%
as.data.frame
# Extraer valores de los raster de ETP reclasificados
spl.rnk <- cbind(extract(etpMR2319, spl) %>% as.matrix %>% as.vector %>% as.factor,
extract(etpSR2319, spl) %>% as.matrix %>% as.vector %>% as.factor) %>%
as.data.frame
# Matriz de confusión
cm <- table(spl.rnk$V2, spl.rnk$V1)
# Numero total de muestras
n = sum(cm)
# Diagonal
diag = diag(cm)
# Exactitud global
OA = sum(diag) / n
# Casos observados por clase
rowsums = apply(cm, 1, sum)
p = rowsums / n
# Casos estimados por clase
colsums = apply(cm, 2, sum)
q = colsums / n
# Indice capa
expAccuracy = sum(p*q)
kappa = (OA - expAccuracy) / (1 - expAccuracy)4.3 Métrica de eficiencia de los patrones espaciales
La Métrica de Eficiencia de los Patrones Espaciales (\(E_{SP}\)) evalúa el ajuste de dos series de rasters sin tener en cuenta las magnitudes sino los patrones espaciales. Se basa en métrica de eficiencia de Kling-Gupta (KGE) y se calcula de la siguiente forma:
\[E_{SP}=1-\sqrt{\left(r_{s}-1\right)^{2}+\left(\gamma-1\right)^{2}+\left(\alpha-1\right)^{2}}\]
Donde, \(R_{s}\) es el coeficiente de correlación por rangos de Spearman, \(\gamma\) es el relación de variabilidad (razón entre los coeficientes variación) y \(\alpha\) es el termino de coincidencia espacial que se calculan de la siguiente manera:
\[r_{s}=1-\frac{6\sum_{1}^{n}d^{2}}{n(n^{2}-1)}\]
\[\gamma=\frac{\frac{\sigma_{mod}}{\mu_{mod}}}{\frac{\sigma_{obs}}{\mu_{obs}}}\]
\[\alpha=1-RMSE(Z_{X_{mod}},Z_{X_{obs}})\]
Donde, \(d\) es la distancia entre rangos, \(\sigma\) es la desviación estandar, \(\mu\) es la media, \(RMSE(Z_{X_{mod}},Z_{X_{obs}})\) es la raíz del error cuadrático medio entre los valores estandarizados (\(Z\)) modelados (\(X_{mod}\)) y observados (\(X_{obs}\)).
# Estandarización de rasters
zetpM2319 <- lapply(etpM2319, function(x){
xv <- values(x)
r <- (x - mean(xv, na.rm = T)) / sd(xv, na.rm = T)
})
zetpS2319 <- lapply(etpS2319, function(x){
xv <- values(x)
r <- (x - mean(xv, na.rm = T)) / sd(xv, na.rm = T)
})
zetpM2319 <- rast(zetpM2319)
zetpS2319 <- rast(zetpS2319)
# Error promedio en los valores estandarizados en cada punto de muestreo
spl.pz <- array(dim = c(nrow(spl), nlyr(etpM2319), 2))
spl.pz[,,1] <- extract(zetpM2319, spl) %>% as.matrix
spl.pz[,,2] <- extract(zetpS2319, spl) %>% as.matrix
spl.pz <- apply(spl.pz, c(1, 2), function(x) (x[1] - x[2])^2)
spl.pz <- rowMeans(spl.pz, na.rm = T) %>% sqrt
# Error
spl$Zerror <- spl.pz
writeVector(spl, paste0(wd, "GIS/SPL.shp"), overwrite = T)
# Valores de los rasters estandarizados
spl.z <- cbind(extract(zetpM2319, spl) %>% as.matrix %>% as.vector,
extract(zetpS2319, spl) %>% as.matrix %>% as.vector) %>%
as.data.frame
# Métrica de eficiencia de los patrones espaciales
rs <- cor.test(spl.val$V1, spl.val$V2, method = "spearman")$estimate
gamma <- (sd(spl.val$V2, na.rm = T) / mean(spl.val$V2, na.rm = T)) /
(sd(spl.val$V1, na.rm = T) / mean(spl.val$V1, na.rm = T))
alpha <- 1 - (apply(spl.z, 1, function(x) (x[1] - x[2])^2 ) %>% mean(na.rm = T) %>% sqrt)
esp <- 1 - sqrt((rs - 1) ^ 2 + (gamma - 1) ^ 2 + (alpha - 1) ^ 2)5 Resultados
En la Figura 5.1 se muestra la ETP calculada por ambos métodos. Se puede ver como el satelite MODIS, pese a que se le hayan removido valores atípicos, presenta valores mucho mas altos de lo que se puede esperar. Estos valores se presentan especialmente alrededor del área metropolitana de Bucaramanga. En la ETP estimada mediante Hargreaves-Samani se ve como hay una clara influencia de la elevación debido a la técnica de interpolación utilizada en las temperaturas.
etpMean <- c(MODIS = mean(etpM2319, na.rm = T), Hargreaves = mean(etpS2319, na.rm = T))
etpMean <- rast(etpMean)
pM <- ggplot() + theme_bw() +
geom_spatraster(data = etpMean[["MODIS"]]) +
scale_fill_whitebox_c(palette = "muted") +
labs(title = "MODIS") +
theme(plot.title = element_text(hjust = 0.5))
pS <- ggplot() + theme_bw() +
geom_spatraster(data = etpMean[["Hargreaves"]]) +
scale_fill_whitebox_c(palette = "muted") +
labs(title = "Hargreaves") +
theme(plot.title = element_text(hjust = 0.5))
ggarrange(pM, pS, ncol = 2)
Figure 5.1: Evapotranspiración Acumulada en 8 días promedio entre 2001 y 2020
Por lo tanto, es de esperarse que la matriz de confusión (Tabla 5.1) no de los mejores resultados, como efectivamente sucede. El índice de exactitud global y el indice kappa arrojaron valores de 0.064 y 0.001, respectivamente. Este pobre resultado en la comparación de las magnitudes es un motivo para pensar que el planteamiento del \(E_{SP}\) tiene sentido al comparar únicamente los patrones espaciales.
lbs <- sprintf("%.1f - %.1f", brks[1:6], brks[2:7])
cm <- as.data.frame.matrix(cm)
rownames(cm) <- lbs
colnames(cm) <- lbs
knitr::kable(cm, digits = 3, caption = 'Matriz de confusión')| 12.9 - 18.5 | 18.5 - 24.1 | 24.1 - 29.7 | 29.7 - 35.3 | 35.3 - 40.8 | 40.8 - 46.4 | |
|---|---|---|---|---|---|---|
| 12.9 - 18.5 | 18 | 165 | 321 | 338 | 409 | 2470 |
| 18.5 - 24.1 | 313 | 813 | 1876 | 2448 | 2093 | 24749 |
| 24.1 - 29.7 | 756 | 1376 | 3764 | 5770 | 5116 | 74666 |
| 29.7 - 35.3 | 1226 | 645 | 3287 | 6328 | 6831 | 112144 |
| 35.3 - 40.8 | 864 | 99 | 1284 | 2195 | 2625 | 60480 |
| 40.8 - 46.4 | 165 | 2 | 180 | 318 | 258 | 7882 |
En la Figura 5.2 se ve la ETP estandarizada. Si bien, visualmente no se presentan mayores diferencias respecto a la Figura 5.1, pueden ser mas fácilmente comparables al estar estandarizados. El \(E_{SP}\) dió un resultado de -0.87, los tres componentes mostraron los siguientes resultados \(r_s=\) 0.07, \(\gamma=\) 0.16 y \(\alpha=\) -0.38.
etpMean <- c(MODIS = mean(zetpM2319, na.rm = T), Hargreaves = mean(zetpS2319, na.rm = T))
etpMean <- rast(etpMean)
pM <- ggplot() + theme_bw() +
geom_spatraster(data = etpMean[["MODIS"]]) +
scale_fill_whitebox_c(palette = "muted") +
labs(title = "MODIS") +
theme(plot.title = element_text(hjust = 0.5))
pS <- ggplot() + theme_bw() +
geom_spatraster(data = etpMean[["Hargreaves"]]) +
scale_fill_whitebox_c(palette = "muted") +
labs(title = "Hargreaves-Samani") +
theme(plot.title = element_text(hjust = 0.5))
ggarrange(pM, pS, ncol = 2)
Figure 5.2: Promedio de la ETP estandarizada.
En la Figura 5.3 se muestra el raíz del error cuadrático medio en cada uno de los mil puntos de muestreo. Se ve como en media montaña se producen los menores errores mientras que en las partes altas o bajas tienden a aumentar.
p <- ggplot() + theme_bw() +
geom_spatvector(data = cnca, fill = NA, color = "#A80084") +
geom_spatvector(data = spl, aes(color = Zerror)) +
scale_color_whitebox_c(palette = "viridi")
p
Figure 5.3: RMSE en los puntos de muestreo
6 Discusión
La implementación de ambas comparaciones de las series de rasters se pudo hacer satisfacotoriamente, desde el punto de vista del procedimiento. Existen algunos aspectos que vale la pena tener en cuenta en futuras aplicaciones de estas técnicas de comparación y los cuales se analizan a continuación.
Comparar los valores directamente como se hizo con la matriz de confusión no resultó adecuado y una posible razón de esto es que es difícil que dos medidas o estimaciones sean totalmente comparables cuando se originaron de maneras tan diferentes. Comparar los patrones espaciales mostró mejores resultados y puede ser una mejor estrategia para aprovechar la información de sensores remotos. El \(E_{SP}\) mostró mejores resultados que los indices de exactitud temática. Aún así, el valor fue negativo, lo cual suele tomarse como un ajuste pobre cuando se analiza el KGE o el NSE. Al ser una aplicación inicial se puede esperar que mejore cuando se haga la comparación de evapotranspiración real (ETR) durante una modelación hidrológica y se calibre el modelo; en todo caso, es posible que no se lleguen a los valores que se suelen esperar (>0.6) de una modelación hidrológica típica en las dos métricas ya mencionadas.
En el análisis espacial de los puntos de muestreo, se vio que los mejores resultados se dieron en la media montaña. Es posible que en elevaciones bajas, las altas temperaturas que hay en la zona hagan que el método de Hargreaves-Samani sobrestime la ETP. Así mismo, en elevaciones altas, las dinámicas de un ecosistema de páramo son muy propias de territorios como el colombiano y puede que no se reproduzcan facilmente por el método.
7 Conclusiones
Las técnicas tradicionales de evaluación temática mediante la matriz de confusión pueden no ser las mejores si se desea usar la información un producto de sensoramiento remoto en una calibración hidrológica. En este sentido la métrica de eficiencia espacial (\(E_{SP}\)) propuesta por Dembélé et. al. (2020) puede ofrecer mejores resultados al comparar los patrones espaciales en lugar de las magnitudes.
Los altos valores en mostrados por MODIS aún después de la remoción de atípicos muestra que el pre-procesamiento de esa información debe refinarse pues no es posible que en una zona se presenten mas de 3000 mm de evapotranspiración en 8 días.
La distribución espacial de los errores también muestra que la estimación de ETP sigue siendo un tema abierto, teniendo en cuenta que es un componente tan importante del ciclo hidrológico y suele ser una entrada dentro de la modelación hidrológica.
8 Anexos
knitr::kable(est, digits = 3, caption = 'Diferencias estadísticas en las series despúes de rellenarlas.')| Estadítica | Estación | Variable | Original | Rellenada | Diferencia | |
|---|---|---|---|---|---|---|
| 1 | Media | 23195180 | TMN_CON | 12.720 | 12.740 | -0.020 |
| 2 | Mínima | 23195180 | TMN_CON | 5.900 | 5.900 | 0.000 |
| 3 | Máxima | 23195180 | TMN_CON | 17.400 | 17.400 | 0.000 |
| 4 | Mediana | 23195180 | TMN_CON | 12.800 | 12.800 | 0.000 |
| 5 | Desviación Estandar | 23195180 | TMN_CON | 1.544 | 1.521 | 0.023 |
| 11 | Media | 23195090 | TMN_CON | 13.921 | 13.945 | -0.023 |
| 12 | Mínima | 23195090 | TMN_CON | 9.000 | 9.000 | 0.000 |
| 13 | Máxima | 23195090 | TMN_CON | 17.800 | 17.800 | 0.000 |
| 14 | Mediana | 23195090 | TMN_CON | 14.000 | 14.000 | 0.000 |
| 15 | Desviación Estandar | 23195090 | TMN_CON | 1.246 | 1.239 | 0.007 |
| 16 | Media | 23195110 | TMN_CON | 19.661 | 19.680 | -0.019 |
| 17 | Mínima | 23195110 | TMN_CON | 14.000 | 14.000 | 0.000 |
| 18 | Máxima | 23195110 | TMN_CON | 23.600 | 23.600 | 0.000 |
| 19 | Mediana | 23195110 | TMN_CON | 19.800 | 19.800 | 0.000 |
| 20 | Desviación Estandar | 23195110 | TMN_CON | 1.040 | 1.045 | -0.005 |
| 21 | Media | 37015020 | TMN_CON | 4.720 | 4.734 | -0.014 |
| 22 | Mínima | 37015020 | TMN_CON | -9.800 | -9.800 | 0.000 |
| 23 | Máxima | 37015020 | TMN_CON | 9.000 | 9.000 | 0.000 |
| 24 | Mediana | 37015020 | TMN_CON | 5.600 | 5.600 | 0.000 |
| 25 | Desviación Estandar | 37015020 | TMN_CON | 2.902 | 2.895 | 0.007 |
| 26 | Media | 23185010 | TMN_CON | 23.701 | 23.719 | -0.018 |
| 27 | Mínima | 23185010 | TMN_CON | 19.000 | 19.000 | 0.000 |
| 28 | Máxima | 23185010 | TMN_CON | 27.800 | 27.800 | 0.000 |
| 29 | Mediana | 23185010 | TMN_CON | 23.800 | 23.800 | 0.000 |
| 30 | Desviación Estandar | 23185010 | TMN_CON | 1.104 | 1.097 | 0.006 |
| 31 | Media | 23195502 | TMN_CON | 18.635 | 18.644 | -0.009 |
| 32 | Mínima | 23195502 | TMN_CON | 15.800 | 15.800 | 0.000 |
| 33 | Máxima | 23195502 | TMN_CON | 21.600 | 21.600 | 0.000 |
| 34 | Mediana | 23195502 | TMN_CON | 18.600 | 18.600 | 0.000 |
| 35 | Desviación Estandar | 23195502 | TMN_CON | 0.776 | 0.781 | -0.004 |
| 36 | Media | 23195180 | TMX_CON | 22.627 | 22.643 | -0.016 |
| 37 | Mínima | 23195180 | TMX_CON | 16.000 | 16.000 | 0.000 |
| 38 | Máxima | 23195180 | TMX_CON | 28.800 | 28.800 | 0.000 |
| 39 | Mediana | 23195180 | TMX_CON | 22.600 | 22.600 | 0.000 |
| 40 | Desviación Estandar | 23195180 | TMX_CON | 1.585 | 1.579 | 0.006 |
| 56 | Media | 37015020 | TMX_CON | 13.962 | 13.967 | -0.005 |
| 57 | Mínima | 37015020 | TMX_CON | 7.200 | 7.200 | 0.000 |
| 58 | Máxima | 37015020 | TMX_CON | 20.000 | 20.000 | 0.000 |
| 59 | Mediana | 37015020 | TMX_CON | 14.000 | 14.000 | 0.000 |
| 60 | Desviación Estandar | 37015020 | TMX_CON | 1.651 | 1.654 | -0.003 |
| 66 | Media | 23195502 | TMX_CON | 25.898 | 25.906 | -0.008 |
| 67 | Mínima | 23195502 | TMX_CON | 19.600 | 19.600 | 0.000 |
| 68 | Máxima | 23195502 | TMX_CON | 31.800 | 31.800 | 0.000 |
| 69 | Mediana | 23195502 | TMX_CON | 26.000 | 26.000 | 0.000 |
| 70 | Desviación Estandar | 23195502 | TMX_CON | 1.510 | 1.505 | 0.005 |
| 71 | Media | 23195180 | TSSM_D | 17.621 | 17.610 | 0.011 |
| 72 | Mínima | 23195180 | TSSM_D | 13.933 | 13.933 | 0.000 |
| 73 | Máxima | 23195180 | TSSM_D | 22.867 | 22.867 | 0.000 |
| 74 | Mediana | 23195180 | TSSM_D | 17.567 | 17.533 | 0.033 |
| 75 | Desviación Estandar | 23195180 | TSSM_D | 1.068 | 1.063 | 0.005 |
| 76 | Media | 24055030 | TSSM_D | 19.032 | 19.029 | 0.003 |
| 77 | Mínima | 24055030 | TSSM_D | 14.450 | 14.450 | 0.000 |
| 78 | Máxima | 24055030 | TSSM_D | 23.950 | 23.950 | 0.000 |
| 79 | Mediana | 24055030 | TSSM_D | 19.000 | 19.000 | 0.000 |
| 80 | Desviación Estandar | 24055030 | TSSM_D | 1.062 | 1.063 | -0.001 |
| 81 | Media | 23195090 | TSSM_D | 18.946 | 18.934 | 0.012 |
| 82 | Mínima | 23195090 | TSSM_D | 15.000 | 15.000 | 0.000 |
| 83 | Máxima | 23195090 | TSSM_D | 23.867 | 23.867 | 0.000 |
| 84 | Mediana | 23195090 | TSSM_D | 18.933 | 18.933 | 0.000 |
| 85 | Desviación Estandar | 23195090 | TSSM_D | 1.086 | 1.090 | -0.003 |
| 91 | Media | 37015020 | TSSM_D | 9.117 | 9.119 | -0.002 |
| 92 | Mínima | 37015020 | TSSM_D | 3.400 | 3.400 | 0.000 |
| 93 | Máxima | 37015020 | TSSM_D | 12.500 | 12.500 | 0.000 |
| 94 | Mediana | 37015020 | TSSM_D | 9.150 | 9.150 | 0.000 |
| 95 | Desviación Estandar | 37015020 | TSSM_D | 0.949 | 0.947 | 0.002 |
| 96 | Media | 23185010 | TSSM_D | 27.899 | 27.872 | 0.027 |
| 97 | Mínima | 23185010 | TSSM_D | 23.300 | 23.300 | 0.000 |
| 98 | Máxima | 23185010 | TSSM_D | 33.067 | 33.067 | 0.000 |
| 99 | Mediana | 23185010 | TSSM_D | 27.900 | 27.867 | 0.033 |
| 100 | Desviación Estandar | 23185010 | TSSM_D | 1.243 | 1.233 | 0.010 |
| 101 | Media | 23195502 | TSSM_D | 21.692 | 21.676 | 0.016 |
| 102 | Mínima | 23195502 | TSSM_D | 17.400 | 17.400 | 0.000 |
| 103 | Máxima | 23195502 | TSSM_D | 25.750 | 25.750 | 0.000 |
| 104 | Mediana | 23195502 | TSSM_D | 21.700 | 21.675 | 0.025 |
| 105 | Desviación Estandar | 23195502 | TSSM_D | 1.060 | 1.056 | 0.003 |
9 Replicabilidad
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] RGENERATE_1.3.7 magrittr_2.0.3 RMAWGEN_1.3.7 vars_1.5-6
## [5] lmtest_0.9-40 urca_1.3-3 strucchange_1.5-3 sandwich_3.0-2
## [9] zoo_1.8-11 MASS_7.3-53 date_1.2-40 chron_2.3-58
## [13] reshape2_1.4.4 ggpubr_0.4.0 ggplot2_3.3.6 openxlsx_4.2.5
## [17] sf_1.0-8 leaflet_2.1.1 raster_3.6-3 sp_1.5-0
## [21] tidyterra_0.3.1 terra_1.6-17 dplyr_1.0.10
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-149 tools_4.0.3 backports_1.4.1 bslib_0.4.0
## [5] utf8_1.2.2 rgdal_1.5-32 R6_2.5.1 KernSmooth_2.23-17
## [9] rgeos_0.5-9 DBI_1.1.3 colorspace_2.0-3 withr_2.5.0
## [13] tidyselect_1.2.0 compiler_4.0.3 cli_3.2.0 labeling_0.4.2
## [17] bookdown_0.29 sass_0.4.2 scales_1.2.1 classInt_0.4-8
## [21] proxy_0.4-27 stringr_1.4.1 digest_0.6.30 rmarkdown_2.17
## [25] pkgconfig_2.0.3 htmltools_0.5.3 fastmap_1.1.0 highr_0.9
## [29] htmlwidgets_1.5.4 rlang_1.0.6 rstudioapi_0.14 jquerylib_0.1.4
## [33] generics_0.1.3 farver_2.1.1 jsonlite_1.8.2 crosstalk_1.2.0
## [37] zip_2.2.1 car_3.1-0 s2_1.1.0 Rcpp_1.0.9
## [41] munsell_0.5.0 fansi_1.0.3 abind_1.4-5 lifecycle_1.0.3
## [45] stringi_1.7.8 yaml_2.3.6 carData_3.0-5 plyr_1.8.7
## [49] grid_4.0.3 lattice_0.20-41 cowplot_1.1.1 knitr_1.40
## [53] pillar_1.8.1 ggsignif_0.6.4 codetools_0.2-16 wk_0.7.0
## [57] glue_1.6.2 evaluate_0.17 data.table_1.14.4 vctrs_0.4.2
## [61] gtable_0.3.1 purrr_0.3.5 tidyr_1.2.1 assertthat_0.2.1
## [65] cachem_1.0.6 xfun_0.34 broom_1.0.1 e1071_1.7-11
## [69] rstatix_0.7.0 class_7.3-17 tibble_3.1.8 units_0.8-0
## [73] ellipsis_0.3.2