jueves, 18 de diciembre de 2008

Análisis del decaimiento de masas forestales en la Sierra de los Filabres

Una de las líneas de investigación en las que estoy más involucrado actualmente es el análisis de los factores causales y la detección de indicadores que permitan predecir el decaimiento de masas forestales de pinos. Este proyecto se desarrolla fundamentalmente en la Sierra de los Filabres (Almería), en donde desde el año 2008 se viene detectando un proceso gradual de decaimiento de los pinares de Pinus sylvestris y Pinus nigra.

Dentro de esta línea de investigación nos hemos planteado las siguientes preguntas:

  • ¿Existen indicadores que nos permitan predecir el decaimiento? En concreto ¿se puede utilizar información espectral de imágenes satelitales para este propósito?
  • ¿Existen factores de incitación (bióticos y abióticos) que favorecen el decaimiento de las masas forestales?
  • ¿Genera el cambio climático una predisposición de estas masas al decaimiento?

Para abordarlas, hemos diseñado el siguiente esquema de trabajo.

Esquema conceptual del estudio del decaimiento en la Sierra de los Filabres.

Para obtener información de campo, en noviembre de 2008 se llevó a cabo un muestreo estratificado por especie (Pinus sylvestris, Pinus nigra) y grado de afectación de las masas utilizando una cartografía de daños generada en años previos a partir de imágenes ASTER de alta resolución (15 x 15 m). En 59 rodales se estimó visualmente el grado de afectación de las masas utilizando la guía Ferreti (1994) para la estimación de daños en especies forestales. Además, se midieron distintas variables relacionadas con la estructura del rodal, como por ejemplo el área basimétrica, la altura media o la altura de copa viva. Este muestreó se complementó con información obtenida a partir de fotografías aéreas de otros 17 rodales con un grado de afectación bajo, que no pudieron ser visitados en campo debido a la cobertura de nieve.

En relación a la detección de indicadores, nos planteamos trabajar con imágenes MODIS con una resolución de 250 x 250 m. Aunque estas imágenes tienen mucha menor resolución que las imágenes Landsat o ASTER, tienen la ventaja de ser gratuitas, de fácil acceso y están disponibles con una periodicidad de 16 días. Estas imágenes contienen información sobre índices de vegetación (NDVI y EVI) que nos pueden dar una idea de la evolución temporal de la actividad fotosintética (NDVI) y la estructura del dosel arbóreo (EVI) de los rodales estudiados. Para obtener dicha información seguimos el siguiente protocolo de trabajo:

  1. Adquisición de las imágenes y obtención de la información de los índices de vegetación para los rodales estudiados.
  2. Filtrado de los datos utilizando la información contenida en la capa de fiabilidad (pixel reliability) de las imágenes MODIS.
  3. Interpolación de los datos que faltan en las series temporales.
  4. Estandarización de los datos.

Con esta información, se pretende generar una serie de indicadores que nos permitan predecir el decaimiento en fases tempranas. Para ello se generarán distintos indicadores a partir de las series temporales de índices de vegetación y se relacionarán con los datos de daños observados en campo mediante técnicas de análisis discriminante y análisis correlacionales.

Con respecto a los factores causales del decaimiento, P.D. Manion (1991) identifica tres grandes tipos de factores: (1) Factores de predisposición, que son aquellos, mayoritariamente abióticos (clima), que debilitan al árbol y/o reducen sus defensas a largo plazo haciéndolo más vulnerable; (2) Factores de incitación, que son factores bióticos y abióticos (p.e. heladas) que aceleran su declive a corto plazo; y (3) Factores de contribución, generalmente patógenos oportunistas (bacterias, hongos, insectos, plantas hemiparásitas como el muérdago), que rematan al árbol previamente debilitado. En este trabajo abordaremos el análisis de factores de predisposición e incitación.

Para abordar los factores de predisposición se pretende utilizar técnicas de regresión de series temporales para relacionar las series de índices de vegetación con series climáticas obtenidas a partir de datos de 16 estaciones meteorológicas próximas al área de estudio.

Finalmente, para el estudio de los factores de incitación, se utilizarán análisis correlacionales del daño observado en campo (variable y) y distintas variables topográficas, espaciales, estructurales y microclimáticas (variables x).

Este estudio se está realizando en el marco del proyecto GESBOME (Gestión Sostenible del Bosque Mediterráneo en un Escenario de Cambio Globla) en colaboración con la Universidad de Córdoba, la Universidad de Granada, y el INIA-CIFOR.

lunes, 8 de diciembre de 2008

¿Contribuyen las abejas a la conservación del medio natural?

La desaparición de la apicultura tradicional de montaña y la aparición de distintas enfermedades que atacan a la abeja desde hace años han desencadenado una disminución muy grande de las poblaciones de abeja melífera en el medio rural. ¿Cuál es el papel ecológico que juega la abeja en los ecosistemas de montaña? Si la abeja desaparece ¿suplen otros polinizadores el papel de la abeja en la polinización de la flora de montaña? ¿Hasta qué punto su ausencia en zonas rurales puede repercutir sobre la conservación de otras especies? Para responder a estas preguntas, se inició a comienzos del 2008 un estudio para evaluar el papel de la abeja melífera en la polinización de distintas especies características de la flora atlántica de montaña en la Cordillera Cantábrica.

Para llevar a cabo este estudio se buscaron colmenares aislados en valles de montaña (para evitar efectos solapados entre colmenares). En estas localidades se seleccionaron especies representativas de la vegetación de montaña y que fueran importantes como recurso alimenticio de otras especies como el oso pardo. Las especies seleccionadas fueron el arándano (Vaccinium myrtillus), el cerezo (Prunus avium), el castaño (Castanea sativa), el endrino (Prunus spinosa) y el majuelo (Crataegus monogyna). Se marcaron distintos individuos de cada una de estas especies a distintas distancias de las colmenas (de 0 a 1500 m). En cada planta (o pie de planta) se seleccionaron varias ramas, la mitad de estas se cubrieron con una malla y la otra mitad no. Con esto se pretendía cuantificar el papel de la abeja frente a otros polinizadores. Desde marzo a octubre, se fué varias veces a campo y se contaron flores y frutos en cada una de las ramas marcadas.

El modelo general que se pretende comprobar es el siguiente:

% Frutos = Distancia * Malla + Error(Localidad)

en donde:

% Frutos = (Frutos / Flores)*100
Distancia = Distancia a la colmena (m)
Malla = Factor con dos niveles (con malla, sin malla)
Localidad = Factor con siete niveles correspondientes a las localidades de muestreo. Este factor se ha de tratar como un factor aleatorio para controlar la variabilidad del experimento.

Una representación de los datos obtenidos en campo (ver figura) parece apuntar a un efecto importante de la abeja en la polinización del arándano (bilberry), el endrino (blackthorn) y el majuelo (hawthorn). En breve se espera poder ajustar modelos estadísticos adecuados para la comprobación de las hipótesis subyacentes.


Este estudio está siendo realizado por la Fundación FIRE y el FAPAS con financiación de la Fundación Biodiversidad.

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

Buscar entradas