martes, 2 de diciembre de 2008

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

7 comentarios:

Germán Mauricio Valencia Hernández dijo...

Saludos Luis,
Gracias por publicar esta guia. Voy a tratar de aplicarla y te comento que tal me fue.

Luis Cayuela dijo...

Gracias Germán y perdona por tardar tanto en contestar. Si encuentras cualquier problema cuando uses este código no dudes en consultarme.

Unknown dijo...

Hola Luis,
Me parece un trabajo muy interesante el que haces, sobretodo el de la sierra de Filabres. Yo estoy trabajando con imágenes MODIS (MOD13Q1 y MOD09Q1)y me gustaría transformar los valores de la banda de VI_Quality en número binarios con la finalidad de extraer los píxeles de NDVI que contenga una mayor calidad durante el periodo considerado. Has aplicado tu este preceso a las imágenes que has utilizado?. Estoy intentando escribir un script en IDL pero no soy capaz de que funciones.
cualquier sugerencia será de gran ayuda.

Un saludo

Belen

Luis Cayuela dijo...

Hola Ana Belén,

la respuesta a tu pregunta es sí. Hay una manera de hacerlo en R y esto está detallado en otra entrada de mi blog (http://luiscayuela.blogspot.com/2009/01/uso-de-imgenes-modis-para-estimar-el.html). El procedimiento es básicamente filtrar los datos con la capa de VI_Quality y luego interpolar los datos faltantes con técnicas estadísticas de interpolación. Todo esto está detallado en el blog. El problema que puedes encontrar al utilizar R en windows (yo trabajo con Linux) es que hay que configurar una serie de herramientas espaciales para que funcionen algunos de los paquetes espaciales de R, como rgdal.

Espero que esta guía te sea de utilidad. Si tienes cualquier pregunta no dudes en escribirme.

Un saludo,

Luis

Fevegal dijo...

Consulta: Necesito extraer los datos de calidad MODIS sds VI_Quality bit:2-5. Luego reproyectar y transformar a geotif. El problema consiste en que se extraen estos datos usando LDOPE tool, pero ya no se puede reproyectar con MRT tool o con R, debido a que ya no tiene ningun archivo de cabecera. Como puedo trabajar con estos archivos?.

Luis Cayuela dijo...

Creo que tu problema no hace relación a esta entrada. La metodología propuesta es para extraer valores de las imágenes MODIS para una serie de puntos (en forma, por ejemplo, de capa shapefile o puntos con coordenadas geográficas). Tu pregunta hace relación al procesamiento de imágenes completas, lo cual es difícil de hacer en R, aunque posiblemente en Grass puedas hacerlo.

Anónimo dijo...

Hola Luis:

Gracias por tu blog. Recién comienzo en el mundo del análisis de imágenes satelitales y quería hacerte una pregunta muy básica. Estoy bajando imágenes MODIS MOD13Q1 de los últimos años para realizar un análisis multitemporal sobre el Valle de Guadalhorce. Pordrías decirme por favor que aplicación es la mejor usar para realizar las correcciones geométricas y radiométricas?? En lo posible si conoces de aplicaciones freeware. Gracias!

Buscar entradas