Uso de imágenes MODIS para estimar el estado de masas forestales (I)

Los índices de vegetación calculados a partir de la información espectral contenida en las imágenes MODIS están diseñados para ofrecer comparaciones espaciales y temporales de las condiciones de la vegetación. La utilización de este tipo de herramientas puede ofrecer, sin embargo, una serie de problemas técnicos en su procesamiento y manejo. En primer lugar las imágenes vienen en un formato especial (*.hdf) que no leen la mayoría de los SIG y softwares de teledetección. En segundo lugar, la información contenida en las distintas bandas de las imágenes MODIS viene en un tipo de proyección (Sinusoidal) bastante atípico, por lo que suele ser necesaria la reproyección de esta información. En el pasado, nuestro grupo de investigacion había utilizado el MODIS import tool y una rutina generada en ArcGIS para procesar un número grande de estas imágenes. El proceso podía llevar uno o dos días y requería cierto grado de supervisión. El uso de software libre, en particular si se integra bajo Linux, ofrece una alternativa mucho más rápida y potente (además de gratuita) para el manejo de información espacial.

Esta guía pretende dar unas pautas para la extracción de información de índices de vegetación para un conjunto de puntos (inventarios de campo). Los pasos descritos implican fundamentalmente el uso del software de libre distribución R, que funciona bien bajo Windows o Linux. El uso de otros paquetes, como GDAL (Geospatial Data Abstraction Library), pueden ser también de gran ayuda a la hora de leer la información de las bandas contenidas en los archivos *.hdf en el que las imágenes MODIS son distribuidas.

Paso 1. Descarga de imágenes MODIS

Las imágenes MODIS pueden solicitarse en la página https://wist.echo.nasa.gov/api/. Es necesario crearse un usuario y una clave. Después habrá que seleccionar en el panel de búsqueda (Search) la región geográfica en la que nos interese trabajar, las fechas en las que queremos buscar imágenes, y el tipo de imágenes requeridas. Para trabajar con índices de vegetación, lo más conveniente es trabajar con las imágenes MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid V005. Es importante que seleccionemos bien la versión V005 ya que la versión V004 no contiene exactamente las mismas bandas, sobretodo en lo que respecta al índice de calidad de la información espectral. Además, la V005 está parcialmente validada en campo, mientras que la V004 no lo está.

Una vez hecha la solicitud, nos llegará a nuestro e-mail una confirmación de que la solicitud se ha procesado correctamente y una dirección de ftp en donde nos podemos descargar las imágenes (una a una o compiladas en un archivo *.zip).

Paso 2. Lectura de la información de cabecera de las imágenes MODIS (opcional)

Se puede mirar la información contenida en las imágenes utilizando GDAL desde Python. Si trabajamos en Linux, se puede acceder a Python desde la terminal usando el comando:

$ python

Luego debemos importar los siguientes paquetes:

$ import gdal
$ import os

Para acceder a la carpeta en donde están las imágenes MODIS en formato *.hdf debemos de cambiar el directorio de trabajo usando los siguientes comandos:

$ os.getcwd()
$ os.chdir("path del directorio de trabajo")

Ahora podemos utilizar los siguientes comandos para leer la información de las bandas contenidas en alguno de los archivos *.hdf:

$ ds=gdal.Open("MOD13Q1.A2000049.h17v05.005.
2006269163246.hdf")
$ sds_md=ds.GetMetadata('SUBDATASETS')
$ sds_md['SUBDATASET_1_NAME]
'HDF4_EOS:EOS_GRID:"MOD13Q1.A2000049.h17v05.005.2006269163
246.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days NDVI'
$ sds_md['SUBDATASET_2_NAME']
'HDF4_EOS:EOS_GRID:"MOD13Q1.A2000049.h17v05.005.2006269163
246.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI'

Y así sucesivamente. Esto nos da los nombres completos de las bandas que tendremos que utilizar cuando usemos GDAL desde R (paquete rgdal) para importarlas como objetos espaciales (SpatialGridDataFrame). Se puede cerrar la sesión en Python con Ctrl + D.

Paso 3. Importar imágenes MODIS en R y asignar valores de índices de vegetación (NDVI, EVI) y calidad (VI Quality) a una serie de puntos tomados en campo

Para llevar a cabo este paso tenemos que abrir la consola de R (que asumimos está instalado). Si trabajamos en Linux podemos hacerlo directamente desde la terminal escribiendo el comando:

$ R

Una vez en R se puede utilizar el siguiente código (las explicaciones intercaladas en el código están en inglés; en mayúsculas se han señalado las líneas que ha de modificar cada usuario del código):

##################################################################################################
 LOAD LIBRARIES
# install.packages(c("sp", "maps", "maptools", "rgdal", "gpclib", "spatstat")) # This has to be done only the first time the code is run in your computer

library(sp)
library(maps)
library(maptools)
library(rgdal)
library(gpclib)
library(spatstat)

##################################################################################################
SET WORKING DIRECTORY AND READ FIELD POINTS FROM WHICH VI VALUES WILL BE EXTRACTED

setwd("NAME OF THE FOLDER WHERE YOUR *.hdf FILES ARE FOUND") # Only *.hdf files can be placed in this folder 
id <- dir(pattern = ".hdf") # Create an object with the names of the files contained in the working directory folder
id <- id[nchar(id) <- max(nchar(id))]   # Select only *.hdf files (excluding *.xlm metadata files) based on the length of character string

# If reading field data from a shapefile, then go to Chunk 1
# If reading field data from a database (*.txt file), then go to Chunk 2

####### Chunk 1

points <- readShapePoints("PATH TO THE *.shp FILE WITH THE FIELD DATA") # Read the shapefile with the field points 
proj4string(points) <- CRS("+proj=utm +zone=30 +ellps=WGS84 +datum=WGS84 +units=m +no_defs") # Assign the projection of the SpatialPointsDataFrame 
points <- spTransform(points, CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")) # Transform the projection 
# Redefine and/or rename the variables contained in the data slot. In this particular case, it is expected to have information on the site id (Id), the level of tree damage 
# (Damage) and the planted species (Species)
slot(points, "data") <- data.frame(Id = points@data$VARIABLE NAME, Damage = points@data$VARIABLE NAME, Species = points@data$VARIABLE NAME) 

####### Chunk 2

points <- read.table("PATH TO *.txt FILE WITH THE FIELD DATA", header = T, sep = "\t") # Read the *.txt file with the field points
coordinates(points) <- data.frame(points$x, points$y) # This converts the data frame in a SpatialPointsDataFrame
proj4string(points) <- CRS("+proj=utm +zone=30 +ellps=WGS84 +datum=WGS84 +units=m +no_defs") # Assign the projection of the SpatialPointsDataFrame
points <- spTransform(points, CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")) # Transform the projection
# Redefine and/or rename the variables contained in the data slot. In this particular case, it is expected to have information on the site id (Id), the level of tree damage 
# (Damage) and the planted species (Species)
slot(points, "data") <- data.frame(Id = points@data$VARIABLE NAME, Damage = points@data$VARIABLE NAME, Species = points@data$VARIABLE NAME)


##################################################################################################
LOAD MODIS *.HDF4 FILES AND EXTRACT VALUES FOR FIELD POINTS (NDVI, EVI, QA)

# Create a list with three identical spatial objects where extracted values from *.hdf images will be assigned to (NDVI, EVI, and pixel reliability respectively)
list_vi <- list(points_ndvi = points, points_evi = points, points_rel = points)  index <- c("NDVI", "EVI", "pixel reliability") # Character vector referring to the layers in the *.hdf files from which the information will be extracted

for (j in 1:length(vi_index)) {
 for (i in 1:length(id)) {
  NDVI_i <- readGDAL(paste("HDF4_EOS:EOS_GRID:", id[i], ":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days ", index[j], sep = ""),     offset=c(1300,3400),region.dim=c(300,500)) # Read the *.hdf file [i] band[j]                                      overlay_i <- overlay(NDVI_i, list_vi[[j]]) # Extract values from VI grid[i,j] using the SpatialPointsDataFrame
  ntemp <- names(slot(list_vi[[j]], "data")) # Assign temporal column names
  slot(list_vi[[j]], "data") <- cbind(slot(list_vi[[j]], "data"), slot(overlay_i, "data")$band1)   # Bind assigned extracted values to existing SpatialPointsDataFrame[[j]]
  ntemp_i <- substr(id[i], 10, 16) # Create column name for extracted values from image[i,j]   names(slot(list_vi[[j]], "data")) <- c(ntemp, ntemp_i) # Assign final column names to SpatialPointsDataFrame[[j]]
 }
}

attach(list_vi)
proj_back <- CRS("+proj=utm +zone=30 +ellps=WGS84 +datum=WGS84 +units=m +no_defs") # Define the original CRS of field points
points_ndvi <- spTranform(points_ndvi, proj_back) # Change projection of points_ndvi
points_evi <- spTranform(points_evi, proj_back) # Change projection of points_evi
points_rel <- spTransform(points_rel, proj_back) # Change projection of points_rel
save(points_ndvi, points_evi, points_rel, file = "PATH TO THE FOLDER WHERE THESE R OBJECTS WILL BE SAVED/points.Rob") # Save SpatialPointsDataFrames