Uso de imágenes MODIS para estimar el estado de masas forestales (II y III)

Las imágenes MODIS MOD13Q01 a una resolución de 250 m y con una periodicidad de 16 días contienen varias capas de información relacionadas con los índices de vegetación. Anteriormente se explicó cómo automatizar el proceso de lectura de estas imágenes y extracción de datos de las capas subyacentes para una serie de puntos georreferenciados en campo. En este apartado se explica cómo: (1) filtrar los datos utilizando la información contenida en la capa de fiabilidad (pixel reliability); y (2) interpolar los datos que faltan (datos poco fiables) en las series temporales.

1. Filtrado de los datos

En la versión 5 de las imágenes MODIS se incorpora una capa (pixel reliability) que contiene información resumida sobre la fiabilidad de los datos, en dónde:

0 = Fiable (sin nubes, sin nieve, sin carga atmosférica de aerosoles);
1 = Semi-fiable (sin nubes, sin nieve, el objetivo es visible y las correcciones atmosféricas han mejorado la calidad del dato);
2 = Nieve (superficie cubierta total o parcialmente por nieve);
3 = Nubes (presencia total o parcial de nubes);
4 = Sin datos (información no procesada por razones diversas).

Es necesario que se use esta información para eliminar datos que pudieran estar afectados por la cobertura de nieve o la presencia de nubes antes de proceder a su análisis. Para ello podemos utilizar el siguiente código de R que nos permitirá hacer esto de manera fácil y rápida. En la sesión anterior, se generaron tres objetos en R (points_ndvi, points_evi, points_rel) que contenían la información de los índices de vegetación NDVI y EVI y el grado de fiabilidad de esta información, respectivamente, para 76 puntos georreferenciados en campo y todas las imágenes MODIS disponibles (201) desde comienzos del año 2000 hasta finales del 2008.

Antes de comenzar, se instalan y cargan los paquetes que van a ser necesarios para el procesamiento de los datos.

install.packages(c("pamr", "cluster", "survival"))
library(pamr)

A continuación se importa el archivo que salvamos en la sesión anterior con los objetos de R que contienen la información de interés (points_ndvi, points_evi, points_rel). Este archivo se ha subido a una ubicación pública por lo que todo el mundo puede acceder a ellos utilizando la siguiente línea de código.

load(url("http://archivos-para-subir.googlegroups.com/web/points.Rob?hl=es&gda=mttqMzwAAACkOfwqHAVd4YqgfIB09GDRafzhFXYXrkFb8NzylP1tVw7qabgqw0xDKbwB-h3MnSf9Wm-ajmzVoAFUlE7c_fAt"))
ls()

[1] "points_evi" "points_ndvi" "points_rel"

Estos tres objetos son objetos espaciales del tipos SpatialPointsDataFrame. La información que nos interesa está contenida en un apartado específico de estos objetos al que se accede nombrando a dicho objeto y postponiéndole @data.

ndvi <- points_ndvi@data
evi <- points_evi@data
rel <- points_rel@data

Esto es, se han creado tres objetos nuevos llamados ndvi, evi y rel, respectivamente, que contienen solamente la información que nos interesa sobre los índices de vegetación NDVI y EVI y su fiabilidad, excluyendo toda referencia espacial. Creamos ahora una función que nos permita resumir la capa de información referente a la fiabilidad de los datos (rel).

table2 <- function(x) {
x1 <- as.data.frame(table(x))
x2 <- data.frame(x = as.factor(c(0:4)))
x3 <- merge(x2, x1, by.x = "x", all = T)[,2]
ifelse(is.na(x3) == T, 0, x3)
}

Y la aplicamos.

relsum <- as.data.frame(t(apply(rel[ , -c(1:4)], 2, table2)))
colnames(relsum) <- c("Reliable", "Semireliable", "Snow", "Cloud", "Faulty data")
apply(relsum, 2, sum)

Reliable Semireliable Snow Cloud Faulty data
12490 2471 288 27 0

En resumen, del total de datos (76 parcelas de trabajo x 201 fechas = 15276), la mayor parte (12490) son totalmente fiables y una proporción importante (2471) son semi-fiables. El problema está en datos con nieve (288) o presencia de nubes (27). Vamos a ver gráficamente en qué fechas se concentran más los datos agrupados por categoría de fiabilidad.

ts.rel <- ts(relsum, start = c(2000, 4), frequency = 23)
#png("Reliability index.png", width = 1000, height = 400)
plot(ts.rel[,3], xlab = "Time", ylab = "Number of field points", col = "blue", type = "l", ylim = c(0,80), lwd = 2, cex.axis = 1.3, cex.lab = 1.5)
par(new = T)
plot(ts.rel[,4], xlab = "", ylab = "", col = "red", type = "l", ylim = c(0,80), axes = F, lwd = 2)
par(new = T)
plot(apply(ts.rel[,c(1,2)], 1, sum), xlab = "", ylab = "", col = "grey", type = "l", ylim = c(0,80), axes = F, lwd = 2)
#dev.off()

Figura 1. Representación del grado de fiabilidad de la información obtenida de las imágenes MODIS para los índices de vegetación en las distintas fechas (gris = fiable o semifiable; azul = nieve; rojo = nubes). Nota: si se quita el símbolo # del código anterior la figura se salvará en el directorio de trabajo establecido como un archivo gráfico en formato *.png.

Por lo que se observa en la Figura 1 la mayoría de los datos con cobertura de nieve se concentran en unas pocas fechas, repartidas entre los meses de invierno (noviembre a febrero) de los distintos años. El último paso sería eliminar estos datos poco fiables (con nieve o nubes) de las capas que contienen la información sobre los índices de vegetación NDVI o EVI. Estos valores se transformarán a 'NA' (non available).

for(i in 5:length(ndvi[1,])) {
ndvi[, i] <- ifelse(rel[, i] > 1, NA, ndvi[, i])
}

for(i in 5:length(evi[1,])) {
evi[, i] <- ifelse(rel[, i] > 1, NA, evi[, i])
}

2. Interpolación de datos

¿Qué hacer ahora con los datos que faltan en las series de índices de vegetación? Una opción razonable es interpolar estos datos. Para interpolar los datos que faltan podemos utilizar la información de otras series que se comportan de modo parecido o utilizar los datos próximos en el tiempo dentro de esa misma serie. Cualquiera de estas opciones es buena.

Vamos a aplicar el uso de estas dos metodologías con los datos de NDVI (ndvi). Lo primero que intentaremos será interpolar los datos que faltan (NA) utilizando información de otras series que se comportan de modo parecido a lo largo del tiempo. Para ello usaremos la función pamr.knnimpute() del paquete pamr. Esta función calcula el valor de un dato que falta promediando los valores para la misma fecha de las series más cercanas. Está función sólo permite interpolar un número máximo de datos para cada fecha. Si faltan muchos datos a la vez para una única fecha entonces el algoritmo no funciona y da un error. Este número, ncol, está definido por defecto en 0.8 (80%). Si tenemos 76 series, esto quiere decir que cuando tengamos aproximadamente 60 NA para una misma fecha, no podremos interpolar los valores utilizando está metodología.

x1 <- ndvi[,-c(1:4)]
f1 <- function(x) sum(is.na(x))
t1 <- apply(x1, 2, f1)
x.imp <- x1[ , t1<60]
str(x.imp)

[1] 76 200

Así que para 200 de las 201 fechas podremos interpolar los NA usando los vecinos más próximos, mientras que para una de las fechas para las que tenemos más del 80% de NA no podremos utilizar esta metodología (ver Figura 1).

El argumento principal que se necesita definir en esta función es una lista con dos componentes: (1) una matriz con las series en las columnas y los sitios en las filas; y (2) un vector de clases para los sitios. Para este último tomaremos los datos de daño (tres clases) medidos en campo, asumiendo que el comportamiento de las series para cada uno de los grupos de daño ha de ser más o menos parecido. Otro argumento, k, permitirá especificar el número de datos que se utilizarán para calcular el valor que falta en una determinada serie.

class <- as.factor(ndvi[, 2])
ndvi.list <- list(x = as.matrix(x.imp), y = class) i
mputed.data <- pamr.knnimpute(ndvi.list, k=5)
x.imp <- imputed.data$x
for (i in 1:length(t1[t1<60])) {
x1[ , t1<60][, i] <- x.imp[, i]
}

El resultado es un arreglo de datos (data.frame) con todos los datos interpolados usando los vecinos más próximos dentro de grupos de sitios definidos por el vector class para todos los datos faltantes, menos para una de las fechas, en donde teníamos más del 80% de datos faltantes y, por tanto, no podíamos utilizar esta metodología. Para interpolar estos datos vamos a utilizar la segunda de las metodologías propuestas, que consiste en interpolar los datos utilizando los valores de la propia serie a lo largo del tiempo. Para ello usaremos la función spline() del paquete stats, que viene instalado por defecto en la instalación básica de R (por lo que no es necesario cargarlo al comienzo de cada sesión).

series.ndvi <- ts(x2, start=c(2000,4), frequency=23)
dates <- time(series.ndvi)
for (i in 1:dim(series.ndvi)[2]) {
if (length(which(is.na(series.ndvi[,i])) > 0 ) ) {
interp <- spline(x=dates, y=series.ndvi[,i], xout=dates[which(is.na(series.ndvi[,i])) ])
series.ndvi[which(is.na(series.ndvi[,i])) , i] <- interp$y
}
}

Por lo que ya tendríamos nuestros datos de NDVI interpolados usando dos técnicas distintas (vecino más próximo y cubic splines). El resultado final es un objeto llamado series.ndvi del tipo ts (time series) con las fechas en filas y los sitios en columnas.