Estima de los atributos funcionales de ecosistemas forestales a partir del NDVI

En 2006, Domingo Alcaraz publicaba un trabajo en Global Ecology and Biogeography sobre la caracterización de los tipos funcionales de ecosistemas de la Península Ibérica a partir del uso de índices de vegetación obtenidos de imágenes de satélite. Este trabajo consiste, básicamente, en definir unos atributos que permitan distinguir desde un punto funcional distintos tipos de ecosistemas. En su trabajo, Alcaraz et al. (2006) definen con estos fines tres atributos: la integral del índice de verdor NDVI (NDVI-I), la diferencia entre el máximo y el mínimo NDVI (RREL) y el mes de máximo NDVI (MMAX) (Figura 1). El primero se relaciona funcionalmente con la productividad de los ecosistemas, mientras que el segundo y el tercero están más vinculados a la estacionalidad. Utilizando estos tres atributos, calculados para una serie temporal larga obtenida a partir de imágenes NOAA/AVHRR a una resolución de 1 x 1 km, Alcaraz et al. (2006) clasificaron todos los tipos de ecosistemas en la Península Ibérica a lo largo de gradientes de productividad y estacionalidad.

Figura 1. Los tres atributos del NDVI empleados en la caracterización de los tipos funcionales de ecosistemas por Alcaraz et al. (2006): la integral del NDVI (NDVI.I), la diferencia entre el máximo y el mínimo NDVI (RREL) y el mes de máximo NDVI (MMAX). Ejemplo tomado de un pinar de Pinus nigra en la Sierra de los Filabres, Almería.

Partiendo de esta idea, pensamos en utilizar estos atributos, no para distinguir distintos tipos funcionales de ecosistemas, sino para intentar detectar el decaimiento en plantaciones de Pinus sylvestris y P. nigra en la Sierra de los Filabres. Para ello seguimos los siguientes pasos:

1. Se promediaron en cada fecha los valores de las series temporales de NDVI obtenidos a partir de las imágenes MODIS (250 x 250 m). Esto supone unos 8 a 9 datos por fecha (2000-2008). Cómo estamos trabajando con valores promedios no hace falta interpolar los datos faltantes como hicimos anteriormente.
2. Se calculó la integral total anual (NDVI-I) como la suma del área bajo la curva en cada intervalo.
3. Se calculó la diferencia entre el mínimo y el máximo NDVI (RREL) y el mes del máximo NDVI (MMAX). Estos cálculos son bastante intuitivos así que no hace falta explicarlos (ver Figura 1).
   

Para los datos de Filabres, usamos el siguiente código en R, que permite calcular estos atributos a partir de los datos de las series temporales de NDVI obtenidos previamente de las imágenes MODIS, y obtener gráficos de dichos atributos para cada una de las 76 parcelas de estudio. Como en ocasiones anteriores, el código permite acceder directamente a los datos y ejecutar las distintas funciones diseñadas para tal fin.

load(url("http://archivos-para-subir.googlegroups.com/web/R+step+1.Rob?gda=GQGm0T4AAACkOfwqHAVd4YqgfIB09GDRsFEbW00qADF89i9HmLDg8UMvt5QpYAf3GGSwB4Eu5X3jsKXVs-X7bdXZc5buSfmx"))
ndvi <- points_ndvi@data[,-c(1:4)]
julday <- as.numeric(substr(labels(ndvi)[[2]], 5, 7)) year <- as.numeric(substr(labels(ndvi)[[2]], 1, 4))

NDVI.I <- rep(0,dim(ndvi)[1])
RREL <- rep(0, dim(ndvi)[1])
MMAX <- rep(0, dim(ndvi)[1])
for (i in 1:dim(ndvi)[1]) {
ndvi.mean <- tapply(as.numeric(ndvi[i, ]), as.factor(julday), mean)/10000
NDVI.i <- abs(ndvi.mean[1] - ndvi.mean[length(ndvi.mean)])/2 + min(c(ndvi.mean[1], ndvi.mean[length(ndvi.mean)]))
for (j in 1:(length(ndvi.mean)-1)) {
NDVI.i <- NDVI.i + abs(ndvi.mean[j+1] - ndvi.mean[j])/2 + min(c(ndvi.mean[j+1], ndvi.mean[j]))
}
NDVI.I[i] <- NDVI.i
RREL[i] <- max(ndvi.mean) - min(ndvi.mean)
MMAX[i] <- dimnames(ndvi.mean)[[1]][ndvi.mean == max(ndvi.mean)]
png(as.character(paste("NDVI.I", i, ".png", sep = ""), width = 550))
dimnames(ndvi.mean)[[1]] <- c("Jan", "", "Feb", "", "Mar", "", "Apr", "", "May", "", "Jun", "", "Jul", "", "Aug", "", "Sep", "", "Oct", "", "Nov", "", "Dec")
plot(c(0, ndvi.mean), type = "n", xlab = "", ylab = "NDVI", axes = F, ylim = c(0.1,0.9))
axis(1, at = c(2:24), labels = dimnames(ndvi.mean)[[1]], las = 2)
axis(2)
polygon(x = c(2:24, 24, 2), y = c(ndvi.mean, 0.1, 0.1), col = "grey90")
segments(x0 = 1, y0 = min(ndvi.mean), x1 = 1, y1 = max(ndvi.mean), col = "red", lwd = 3)
pos <- pmatch(ndvi.mean[ndvi.mean == max(ndvi.mean)], ndvi.mean)
segments(x0 = pos +1, y0 = 0.1, x1= pos+1, y1 = as.numeric(ndvi.mean[ndvi.mean == max(ndvi.mean)]), col = "blue", lwd = 3, lty = 3)
box()
title(main = paste(points_ndvi@data$Species[i], points_ndvi@data$Id[i]), sub =paste("Level of damage =", points_ndvi@data$Damage[i]))
text(x = 12, y = mean(c(max(ndvi.mean), 0.1)), labels = "NDVI.I", cex = 1.5)
text(x = 1.5, y =max(ndvi.mean) + 0.05, labels = "RREL", cex = 1.2)
text(x = ifelse(pos < y ="ifelse(pos" labels = "MMAX" cex =" 1.2)
dev.off()
}

El resultado de aplicar este código son tres vectores con los valores de estos atributos para cada una de las 76 parcelas (NDVI.I, RREL, MMAX) y 76 gráficas con la representación de dichos valores.

Estos valores se utilizarán junto a otras variables (p.e. insolación, altitud, etc) para intentar discriminar masas forestales con distintos grados de decaimiento. Pronto esperamos poder mostrar avances en esta línea, en la que llevamos trabajando ya varios meses.

Referencias:

Alcaraz, D., Paruelo, J.M. & Cabello, J. 2006. Identification of current ecosystem functional types in the Iberian Peninsula. Global Ecology and Biogeography 15: 200-212.