Más sobre fumigaciones aéreas y procesionaria

Nuestro artículo publicado recientemente en Forest Ecology and Management sobre la inutilidad de las fumigaciones aéreas para el control de la procesionaria ha suscitado algunas críticas. La principal crítica (reiterada por el revisor de otro artículo que tenemos actualmente en revisión en Climatic Change utilizando esta misma base de datos) se refiere al hecho de que nuestro artículo sólo se refiere a rodales que sufren un grado de infestación de procesionaria alta (nivel 3 o superior sobre una escala de 5). En estos casos es lógico pensar que la plaga ya está al borde del colapso poblacional y, que por tanto, los rodales fumigados van a tener una respuesta similar a los no fumigados: esto es, un colapso da las poblaciones de procesionaria al siguiente año (se puede encontrar una explicación más detallada en esta otra entrada). Sin embargo, puede haber rodales que tengan un grado de infestación medio (por ejemplo 2) y que también sean fumigados. Se podría pensar que en estos rodales la respuesta no es la misma, y que mientras que en los fumigados se rompería el ciclo poblacional, los no fumigados tendrían más probabilidades de sufrir una superpoblación al siguiente año.

Explicaré a continuación qué parte de esta crítica tiene sustento y qué parte no la tiene. Cuando comenzamos este trabajo, tomamos la información de la Consejería de Medio Ambiente de la Junta de Andalucía. En las propias directrices de actuación de la Junta se establecía que sólo los rodales con un nivel de infestación de 3 o más son tratados, mientras que los que sufrían un nivel 2 sólo eran tratados si estaban próximos a un rodal con un grado de infestación de 3 o más. Esta información está disponible en el siguiente enlace (ver páginas 58 a 64). Como esta era información oficial, la dimos por válida sin cuestionarla. Más tarde, cuando surgieron las críticas, procedimos a verificar que, efectivamente, se fumigaban mayoritariamente rodales con un grado de infestación de 3 o más. La siguiente tabla muestra, para el período 2002-2005 (que es aquel para el que disponemos de la información más completa), el número de rodales fumigados (Trat.) y sin fumigar.

	>=3		2		1		0		Total
	Trat.	Sin	Trat.	Sin	Trat.	Sin	Trat.	Sin	Trat.	Sin
2002	47	271	23	414	44	985	29	1550	149	4240
2003	29	371	30	629	38	1074	36	1404	135	4254
2004	29	409	22	482	23	1010	39	1666	118	4271
2005	34	184	16	361	28	919	31	2132	111	4278

A partir de estos valores podemos calcular la odds ratio de la prevalencia de rodales tratados frente a no tratados para cada año entre las categorías de daño 3 o más y el resto. Para ello, haríamos el siguiente cálculo (como ejemplo tomaremos las categorías 3 o más y 2):

A = Rodales daño >= 3 Trat. / Rodales daño >= 3 Sin

B = Rodales daño = 2 Trat. / Rodales daño = 2 Sin

Odds ratio = A / B

Un valor positivo de la odds ratio indicaría que la presencia de tratamientos en rodales con un nivel daño de 3 o más ocurre con ese valor más frecuentemente que en rodales con un nivel de daño de 2. Por ejemplo, una odds ratio de 5 indicaría que los rodales con daño 3 o más se tratan 5 veces más frecuentemente que los rodales con daño 2. Hacemos estos cálculos comparando el nivel de daño de 3 o más con el resto. Los resultados se muestran a continuación.

	OR >=3 / 2	OR >=3 / 1	OR >=3 / 0
2002	3.12	3.88	9.26
2003	1.63	2.21	3.04
2004	1.55	3.11	3.03
2005	4.17	6.06	12.71

De todo esto se deduce lo siguiente. Primero, no todos los rodales que sufren un nivel de daño alto son tratados, como decíamos en nuestro trabajo. Sólo un porcentaje relativamente pequeño como vemos en la primera tabla. Esto es tranquilizador en parte, porque como hemos demostrado en nuestro trabajo, las fumigaciones producen exactamente los mismos resultados que el no hacer nada (la razón es que el propio insecto es controlado por la escasez de recursos alimenticios y el aumento de predadores y parasitoides). Segundo, es cierto que los rodales con un nivel de daño intenso son tratados con mayor frecuencia que los rodales que tienen un nivel de daño medio o incluso bajo, pero resulta curioso ver qué todavía hay muchos rodales con niveles de daño bajo que son tratados ¿qué criterios -más allá de la proximidad a rodales con un grado de infestación alto- utilizan los técnicos para decidir que estos rodales merecen ser tratados? Por último, alguien podría pensar que esto podría invalidar los resultados de nuestro trabajo, pero hemos repetido los análisis publicados en el artículo de Forest Ecology and Management, utilizando sólamente los rodales con un nivel de daño 2. Los resultados, que no muestro aquí por brevedad, muestran claramente que tampoco hay diferencias en la respuesta cuando comparamos rodales fumigados y no fumigados con un nivel de daño moderado. Por tanto, las conclusiones de nuestro estudio siguen siendo válidas.

Un festín de orugas

La procesionaria del pino (Thaumetopoea pityocampa) es un lepidóptero típico de la región Mediterránea. La mariposa de la procesionaria se aparea en verano. La hembra pone sus huevos sobre las copas de los árboles y 30 ó 40 días después nacen las orugas, generalmente en los meses de agosto y septiembre, que construyen sus nidos sobre las ramas y pasan el invierno en ellos. Entre febrero y abril descienden al suelo, forman las características filas indias –de ahí su nombre común de “procesionaria”– y se entierran finalmente en el suelo, donde pasan a la fase de crisálida. En verano las crisálidas eclosionan y surgen las mariposas, que se aparean y reinician de nuevo el ciclo. Durante el invierno, las orugas se alimentan de las hojas de los pinos en los que construyen sus nidos y esa es precisamente la causa de la defoliación.


A pesar de su toxicidad, existen varias especies de aves que han desarrollado mecanismos y estrategias para poder alimentarse de las larvas de procesionaria. Así, por ejemplo, el críalo europeo (Clamator glandorius) y el cuco (Cuculus canorus) son capaces de regurgitar los pelos urticantes de las orugas. El herrerillo capuchino (Lophophanes cristatus) y el carbonero común (Parus major) y garrapino (P. ater) no son capaces de ingerir la oruga entera, sino que las pelan como pipas, quitándoles la cabeza y el tegumento para alimentarse solamente de la parte carnosa de la larva. Estos últimos son los principales responsables de los agujeros que vemos en los nidos de la procesionaria. Existen otras especies que se alimentan de la procesionaria, pero no durante la fase de oruga, como la abubilla (Upupa epops; ver foto más abajo) que desentierra las crisálidas del suelo, o el chotacabras gris (Caprimulgus europaeus), que es capaz de cazar a la mariposa al vuelo durante su efímera existencia (generalmente no más de 24 horas).


Todas estas especies actúan como agentes efectivos para el control natural de este insecto, evitando así con su presencia la aparición de superpoblaciones de procesionaria en los pinares mediterráneos.

Barbaro, L., & Battisti, A. (2011). Birds as predators of the pine processionary moth (Lepidoptera: Notodontidae) Biological Control, 56, 107-114 : doi:10.1016/j.biocontrol.2010.10.009

Gestión de plagas en Andalucia: ¿Es efectiva la fumigación de rodales para combatir la procesionaria?

Anualmente la Consejería de Medio Ambiente de la Junta de Andalucía invierte cerca de un millón o millón y medio de euros en la fumigación de rodales atacados por procesionaria. Sin embargo, cabe preguntarse si dicha medida es efectiva y si la disminución de los brotes de procesionaria responde a una respuesta a la fumigación o, por el contrario, forma parte de los ciclos naturales de esta plaga. Responder a esta pregunta es, por tanto, de gran importancia para la gestión si se quieren minimizar los costes de tratamientos que son, a veces, tan caros como ineficaces.

Para ello, nuestro grupo de investigación ha trabajado con una base de datos recopilada por la Junta de Andalucía (una de las aplicaciones desarrolladas para la Red de Daños y Equilibrios) que contiene información sobre el grado de afectación de los rodales forestales y de los tratamientos aplicados sobre dichos rodales desde 1992 hasta la fecha. Anteriormente se analizó el potencial de una versión preliminar de esta base de datos para el análisis de patrones espacio-temporales de la procesionaria en Andalucía. En este estudio, nos centraremos en la respuesta de la procesionaria en rodales sometidos a fumigación y no tratados. La metodología utilizada es sencilla y consta de los siguientes pasos:

1. Selección de rodales con un grado de afectación de 3 o más (la escala ordinal utilizada va de 0 a 5, en dónde 3 = defoliación fuerte y bolsones numerosos en los bordes del rodal y algo de defoliación en el centro del rodal) que hayan sido sometidos en el otoño de ese mismo año (septiembre-octubre) a fumigación aérea con objeto de controlar la plaga (los datos más completos se tienen sólo para el período 2002-2005 y es con estos datos con los que se han realizado los sucesivos análisis). Hay que puntualizar que, en principio, sólo se fumigan rodales con un grado de afección de 3 o más, y por eso tomamos este criterio a la hora de seleccionar los rodales para el análisis.
2. Cálculo de un "índice de recuperación", que se calcularía como la diferencia entre el grado de afectación de un año y el siguiente. Si un año el rodal ha sido asignado a un nivel de daño 3 y en el siguiente año, el grado de afección disminuye, el índice tendría un valor negativo e indicaría una buena recuperación del rodal.
3. Selección de rodales con un grado de afectación de 3 o más que no hayan sido fumigados. Las muestras se aparean de tal forma que, en cada año, cada rodal con un grado de afectación de 3 o más fumigado se "empareja" con el rodal más próximo geográficamente que también haya tenido un grado de afectación de 3 o más pero que no haya sido fumigado.
4. Se utiliza un test de la t pareado para comparar los índices de recuperación de rodales tratados y no tratados en cada una de las cuatro especies principales de pino (Pinus halepensis, P. nigra, P. pinaster y P. pinea). Para P. sylvestris no hubo una muestra suficientemente representativa cómo para realizar un test estadístico.

Los resultados mostraron que no hay diferencias significativas en el índice de recuperación de rodales fumigados y no fumigados para las cuatro especies: (Figura 1) Pinus halepensis (p-value = 0.7738), P. nigra (p-value = 0.6987), P. pinaster (p-value = 0.2939) y P. pinea (p-value = 0.4793). Esto indica claramente que, cuando menos, no hay evidencias de que los rodales sometidos a fumigación se recuperen antes que los rodales no fumigados y por tanto cuestiona la validez de esta práctica tan habitualmente usada en la gestión forestal.

Figura 1. Distribuciones del índice de recuperación de rodales forestales fumigados (línea contínua) y no fumigados (línea discontínua). Las barras verticales indican la media de las distribuciones de este índice para cada grupo en cada una de las cuatro especies analizadas: P. halepensis (n = 19 pares de rodales), P. nigra (n = 26 pares de rodales), P. pinaster (n = 14 pares de rodales), P. pinea (n = 71 pares de rodales).

Por otro lado, la Comisión Europea quiere prohibir con carácter general, o como mínimo restringir al máximo, el uso de la fumigación aérea en la UE por los daños que este método puede provocar sobre el medioambiente y la salud de las personas.

Queda claro el mensaje, si no es efectivo y encima causa riesgos innecesarios ¿para qué se siguen fumigando miles de hectáreas cada año?

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.

Presentación al 5º Congreso Forestal Español

Algunos de los resultados obtenidos hasta el momento en relación al uso de imágenes MODIS (250 x 250 m) para detectar el decaimiento en repoblaciones de pino silvestre y pino salgareño en la Sierra de los Filabres (Almería) van a ser presentados al 5º Congreso Forestal Español, que tendrá lugar en Ávila, del 21 al 25 de septiembre de 2009. A continuación se muestra el resumen y el texto completo en pdf.

Detección del decaimiento en masas de Pinus sylvestris L. y P. nigra Arnold. procedentes de repoblación en la Sierra de los Filabres (Almería) mediante el uso del sensor MODIS

Durante las últimas décadas se han detectado en Europa y Norte América numerosos casos de decaimiento en masas forestales asociados, de forma general, con la contaminación, la presencia de plagas forestales y el cambio climático. En España, se han descrito procesos de decaimiento en abetales, encinares y pinares. Se prevé que el cambio climático resulte en una intensificación de estos procesos, por lo que la detección temprana de este fenómeno es un paso esencial para la gestión forestal sostenible. El diseño de modelos con base biológica para llevar a cabo dicha detección en Andalucía es uno de los objetivos prioritarios del proyecto GESBOME.

En este estudio se muestran resultados preliminares sobre la aplicación de imágenes MODIS (resolución espacial ~ 250 m) como herramienta para el seguimiento del estado vegetativo de masas forestales en la Sierra de los Filabres (Almería). En esta localidad, el problema del decaimiento de pinares de Pinus nigra Arnold. y Pinus sylvestris L. se viene observando desde el año 2001.

Para ello, se seleccionaron 76 puntos de control en campo y se asignaron a tres niveles de daño: sin afectar, moderado e intenso. En cada parcela se extrajeron los valores del índice de vegetación NDVI (Normalized Difference Vegetation Index) de las imágenes MODIS cada 16 días desde comienzos del año 2000 hasta agosto de 2008. Las series temporales de los tres niveles de daños se compararon entre sí utilizando un índice de referencia basado en la mediana del total de las series temporales para todas las parcelas. Las señales espectrales permitieron discriminar correctamente el nivel de daño extremos de los daños intermedios y moderado, especialmente en pinares de Pinus sylvestris, por lo que pueden utilizarse como indicadores para la detección del decaimiento. El uso de series temporales basadas en índices espectrales puede contribuir además a explorar procesos causales del decaimiento forestal mediante la comparación con variables edafoclimáticas y/o fisiográficas.

Análisis de los factores de incitación del decaimiento en masas forestales

En la estructura inicial planteada para analizar el proceso de decaimiento en la Sierra de los Filabres, Almería, se propuso la exploración de los factores de incitación mediante el análisis correlacional de distintas variables (topográficas, estructurales, espaciales) con el porcentaje de decaimiento estimado de visu en campo.

Uno de los problemas técnicos que se planteó previo a la realización de dichos análisis fue el cálculo de las superficies de insolación, problema que se ha resuelto mediante la aplicación de modelos específicos de insolacióimplementados en GRASS. En una entrada reciente en este blog se proveen los detalles técnicos para la preparación de la bases de datos points. Esta base contiene la siguiente información:

1. Porcentaje de decaimiento observado in situ en campo.
2. Especie (Pinus nigra y P. sylvestris).
3. Área basimétrica del rodal (var. estructural).
4. Densidad de pies del rodal (var. estructural).
5. Coordenadas geográficas x e y (var. espaciales).
6. Elevación (var. topográfica).
7. Pendiente (var. topográfica).
8. Orientación (var. topográfica).
9. Insolación en el solsticio de invierno (var. topográfica).
10. Insolación en el solsticio de verano (var. topográfica).
11. Índice de humedad del suelo (var. topográfica).
    

La información referente a variables topográficas (6-10) ha sido recalculada a tres resoluciones distintas (10x10 m, 30x30 m, 50x50 m). El modelo digital de elevaciones, del que se deriva el resto de información, está a una resolución original de 10x10 m. Sin embargo, dado el grano de las imágenes MODIS utilizadas (250x250 m) y el posible error de localización de los puntos puede ser conveniente re-escalar la información a escalas más groseras como 30x30 m o 50x50 m.

¿Está esta información proyectada a distintas resoluciones muy correlacionada entre sí? Si es así (p.e. r = 0.95), no importará mucho a qué resolución midamos las variables. De lo contrario, habrá que investigar si el cambio de resolución tiene un efecto sobre los resultados de los análisis de correlación con la variable respuesta (i.e. decaimiento).

Las tres bases de datos (points10, points30, points50) con la información referente a variables topográficas calculada a tres resoluciones distintas (10x10 m, 30x30 m y 50x50 m respectivamente) se puede descargar en R con el siguiente código (sólo la base de datos de 30x30 contiene además la información referente a la variable índice de humedad del suelo).

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

[1] "points10" "points30" "points50"

Representamos ahora la matriz de correlaciones.


Cómo podemos observar, el cambio de resolución en la elevación no tiene ningún efecto perceptible sobre los valores de esta variable. Sin embargo, en la pendiente y orientación sí que se observa un cambio importante en los valores de estas variables según se cambia la resolución, sobretodo cuando se cambia la resolución de 10x10 m o 30x30 m a 50x50 m.

Vamos a seleccionar una resolución de 30x30 m para las variables topográficas, ya que es más representativa de los procesos a escala de rodal que la de 10x10 m y, al mismo tiempo, no se pierde mucha información con respecto a una escala de más detalle, cosa que sí pasa con la resolución de 50x50 m. Representamos ahora gráficamente los datos.
Parece que las variables que influyen más como factores de incitación en el proceso de decaimiento son la elevación, sobretodo para el caso de Pinus sylvestris, y la insolación invernal. Relaciones con las variables topográficas reescaladas a otra resolución (10x10 m, 50x50 m) no cambian apenas los resultados obtenidos. Es necesario, no obstante, analizar en mayor profundidad estos datos.

Como utilizar el potencial de GRASS desde R

R tiene varios paquetes que pueden utilizarse para el análisis y visualización de datos espaciales (sp, maptools, rgdal, etc.). Sin embargo, R no es un SIG y, por tanto, tiene ciertas limitaciones a la hora de trabajar con información espacial. Por ejemplo, para algo tan simple como extraer la pendiente y la orientación a partir de un modelo digital de elevaciones, R necesita del paquete RSAGA, que a su vez necesita del software gratuito SAGA diseñado por Alexander Brenning, que funciona únicamente bajo Windows. Esto supone ciertas limitaciones al uso de R como SIG. Estas limitaciones no existen cuando se utiliza GRASS en combinación con R. GRASS ofrece un complemento perfecto a R para la implementación de operaciones típicamente realizadas en SIG y R permite potenciar las capacidades de GRASS para la implementación de análisis geoespaciales.

Para poder hacer uso de ambos programas en necesario tener instalados GRASS y una versión de R posterior a 2.1.0. Desde el 'shell' se puede abrir una sesión de R. Lo primero será instalar los paquetes necesarios para poder utilizar R en combinación con GRASS (sp, rgdal, maptools, spgrass6, spGDAL, spmaptools).

GRASS 6.3.0 (Filabres):~ > R

install.packages(c("sp", "rgdal", "maptools"), dependencies= TRUE)

rS <- "http://r-spatial.sourceforge.net/R" install.packages(c("spgrass6", "spGDAL", "spmaptools"), repos=rS, dependencies=TRUE)

Sin salir de la sesión de R, debemos de cargar el paquete spgrass6 y la 'location' de GRASS donde se encuentran todas las capas de información raster o vectorial.

library(spgrass6)
G <- gmeta6()
str(G)

Si ahora queremos acceder a alguna de las capas raster del 'location' en el que estamos trabajando dentro de GRASS tendremos que utilizar la función de R readCELL6sp() ò readRAST6(). Esto genera un objeto del tipo SpatialGridDataFrame.

mde <- readCELL6sp("mde")

Sin embargo, es importante tener en cuenta que si cargamos muchas capas simultáneamente, especialmente si son capas con un gran número de celdas, se puede agotar la memoria virtual. En este caso R dará el mensaje 'vector memory exhausted'. Cuando esto me ha ocurrido a mi, R se ha quedado colgado y he tenido que reiniciar la sesión. Así que lo que aconsejo es ir eliminando los objetos que ya no necesitemos una vez que extraigamos de ellos la información necesaria. Para conocer el uso que se va haciendo de la memoria se puede utilizar la función gc().

gc()

used (Mb) gc trigger (Mb) max used (Mb)
Ncells 279590 7.5 531268 14.2 361111 9.7
Vcells 3160949 24.2 15608917 119.1 30379099 231.8

Es decir que se están utilizando alrededor de 31 MB de memoria virtual (hay que mirar la primera columna). Ahora eliminamos la capa mde.

rm(mde)
gc()

used (Mb) gc trigger (Mb) max used (Mb)
Ncells 279513 7.5 531268 14.2 361111 9.7
Vcells 137524 1.1 12487133 95.3 30379099 231.8

Y comprobamos que disminuye considerablemente el uso de la memoria virtual.

Antes de acceder a las capas raster, vamos a leer una capa de puntos en formato *.txt (parcelas muestreadas en campo) para la cuál extraeremos la información de elevación, pendiente, altitud e insolación de las capas raster.

points <- read.table(url("http://archivos-para-subir.googlegroups.com/web/Analisis+correlacionales.txt?gda=0ywgME4AAACkOfwqHAVd4YqgfIB09GDRN3No94DJejqv6LWrffOcUtQTec6xGbbz9B_IWZmp37itXp5Ud9d8afGj09bAbTSQ47Cl1bPl-23V2XOW7kn5sQ"), header = T, sep = "\t", dec = ",")
str(points)

'data.frame': 76 obs. of 7 variables:
$ x : num 543726 543486 542903 535211 539086 ...
$ y : num 4124164 4124224 4122417 4118071 4121908 ...
$ Id : Factor w/ 76 levels "1_12","1_13",..: 1 2 4 6 7 8 9 10 11 12 ...
$ Especie : Factor w/ 2 levels "Pinus nigra",..: 1 1 1 1 1 1 1 2 1 1 ...
$ Defoliacion: num 55 55 30 55 30 55 55 60 45 55 ...
$ Area_basim : num 25.8 28.1 24.2 22.0 18.0 ...
$ Densidad : num 286 286 1592 923 1974 ...

El siguiente paso es convertir esta capa de puntos en un objeto espacial del tipo SpatialPointsDataFrame y asignarle una proyección (ha de ser la misma que la de las capas raster y sino habría que reproyectar esta capa).

coordinates(points) <- data.frame(points$x, points$y)
proj4string(points) <- CRS("+proj=utm +zone=30 +ellps=intl +units=m +no_defs")

Y ahora sí, importamos cada una de las capas raster de GRASS a R, extraemos los valores de las capas raster para los puntos incluidos en la capa vectorial points y, por último, borramos las capas raster para no agotar la memoria virtual.

mde <- readCELL6sp("mde")
extract.mde <- overlay(mde, points)
rm(mde)

slope <- readCELL6sp("slope")
extract.slope <- overlay(slope, points)
rm(slope)

aspect <- readCELL6sp("aspect")
extract.aspect <- overlay(aspect, points)
rm(aspect)

rad.summer <- readCELL6sp("beam_radiation_172")
extract.rad.sum <- overlay(rad.summer, points)
rm(rad.summer)

rad.winter <- readCELL6sp("beam_radiation_354")

extract.rad.win <- overlay(rad.winter, points)
rm(rad.winter)

wetness <- readCELL6sp("wetness_index")
extract.wetness <- overlay(wetness, points)
rm(wetness)

Y ahora juntamos toda la información extraida en un único data.frame.

points@data <- cbind(points@data, mde = extract.mde@data, slope = extract.slope@data, aspect = extract.aspect@data, rad.summer = extract.rad.sum@data, rad.winter = extract.rad.win@data, wetness = extract.wetness)

El resultado final, es por tanto, un data.frame con información sobre valores de decaimiento (incluidos en el archivo points) y toda una serie de variables físicas obtenidas a partir de capas raster (pendiente, altitud, etc) que podemos analizar estadísticamente en R, cosa que en principio no podríamos hacer en GRASS. En la siguiente entrada se muestra un análisis preliminar de estos datos en R.

Evaluación del grado de afección de pinares por medio de series temporales MODIS

En el presente trabajo se reproduce la metodología propuesta por Verbesselt et al. (en prensa) para evaluar el potencial de series temporales de índices de vegetación NDVI (Normalized Difference Vegetation Index) y EVI (Enhanced Vegetation Index) derivados de imágenes MODIS 13Q01 v.05, obtenidas con una periodicidad de 16 días. En sesiones anteriores se explicó: (1) como automatizar la lectura de las imágenes MODIS y la extracción de los valores de los índices de vegetación NDVI y EVI a partir de una serie de puntos georreferenciados y visitados en campo; (2) cómo filtrar los datos utilizando la capa de fiabilidad provista en las imágenes MODIS y cómo interpolar estos datos utilizando distintas técnicas de interpolación. En esta sesión se analizan los datos ya filtrados e interpolados. El código en R que genera los resultados que se muestran en esta sesión se puede descargar aquí.

La metodología propuesta en el citado trabajo utiliza el índice relativo propuesto por Díaz-Delgado & Pons (2000). Dicho índice se basa en una serie temporal de referencia que representa el estado de salud promedio de las plantaciones. Esta serie temporal se obtiene calculando la mediana de los valores de los índices de vegetación para todas las parcelas de estudio en cada una de las fechas. La fórmula es la siguiente:

RIt = (VIt/VIt reference) - 1

Donde RIt es el valor de referencia en tiempo t, VIt es la mediana de los índices de vegetación de las parcelas de estudio de una clase de daño determinada, y VIt reference es la mediana de los índices de vegetación de todas las parcelas de estudio.

La estimación de los daños en campo se hizo utilizando la guía Ferreti (1994) en el interior de rodales más o menos homogéneos, a por lo menos 100 m del camino más próximo, y haciendo un recorrido en un radio de aproximadamente 50 m alrededor del punto seleccionado. Se seleccionaron 76 parcelas de campo, 35 en rodales de Pinus nigra y 41 en rodales de Pinus sylvestris. Las clases de daño se agruparon de la siguiente forma: (1) de 0-30% de daño (daño leve); (2) de 30-60% de daño (daño moderado); (3) >60% daño (daño severo). Las parcelas de Pinus nigra tuvieron la siguiente representación por clase de daño: (1) 16 parcelas con daño entre 0-30%; (2) 18 parcelas con daño entre 30-60%; (3) 1 única parcela con daño >60%. Las parcelas de Pinus sylvestris tuvieron la siguiente representación por clase de daño: (1) 12 parcelas con daño entre 0-30%; (2) 12 parcelas con daño entre 30-60%; (3) 17 parcelas con daño >60%. La clasificación en tres categorías de daño es ciertamente arbitraria y se corresponde en parte con la realidad del área de estudio, Filabres, en dónde prácticamente no existen ya rodales exentos de daño. Estas categorías contrastan notablemente con las categorías definidas por Verbesselt et al. (en prensa), en donde la categoría de daño severo se define >30% de daño.

La figura 1(a) muestra la serie temporal de las medianas del NDVI para cada una de las clases de daño agrupadas por especie. El NDVI es un índice que indica la actividad fotosintética de las plantas, o en este caso, de las masas forestales.

Figura 1a. Series temporales del NDVI obtenidas a partir de la mediana de los valores del NDVI de todas las parcelas agrupadas en cada clase de daño.

La figura 1(b) muestra la serie temporal de las medianas del EVI para cada una de las clases de daño agrupadas por especie. El EVI es un índice que refleja bien la estructura del dosel.

Figura 1b. Series temporales del EVI obtenidas a partir de la mediana de los valores del EVI de todas las parcelas agrupadas en cada clase de daño.

De ambas figuras se pueden sacar las siguientes conclusiones: (1) en Pinus sylvestris es posible detectar cierto alejamiento de los series temporales de índices de vegetación para la clase de daño más severa con respecto a las otras dos clases, sobretodo a partir de los años 2004 y 2005 en adelante; (2) el NDVI capta mejor estas diferencias que el EVI; (3) en Pinus nigra no se ve una tendencia clara. En particular para la clase de daño más severa, el comportamiento es bastante caótico, posiblemente como consecuencia de que sólo hay una parcela que representa esta clase de daño. Contrario a lo que cabría esperar, la clase de daño leve tiene en promedio menor actividad fotosintética que la clase de daño moderado.

La figura 2a muestra la serie temporal del índice relativo (RIt) de NDVI para cada una de las clases de daño agrupadas por especie.

Figura 2a. Series temporales del índice relativo de NDVI (RI).

La figura 2b muestra la serie temporal del índice relativo (RIt) de EVI para cada una de las clases de daño agrupadas por especie.

Figura 2b. Series temporales del índice relativo de EVI (RI).

A continuación, se han agrupado los valores relativos de los índices de vegetación (NDVI y EVI) en cada una de las clases de daño por mes y año, y se han representado gráficos de cajas para ver si las distintas clases de daño se diferencian mejor en determinadas épocas del año y/o en determinados años para cada una de las especies estudiadas .

La figura 3a muestra los gráficos de cajas de los valores del índice relativo del NDVI para Pinus nigra y Pinus sylvestris.

Figura 3a. Gráficos de cajas del índice relativo de NDVI (RI) agrupado por meses.

La figura 3b muestra los gráficos de cajas de los valores del índice relativo del EVI para Pinus nigra y Pinus sylvestris.

Figura 3b. Gráficos de cajas del índice relativo de EVI (RI) agrupado por meses.

Finalmente, las figuras 4a y 4b muestran los gráficos de cajas de los valores del índice relativo del NDVI y EVI, respectivamente, agrupados por años.

Figura 4a. Gráficos de cajas del índice relativo de NDVI (RI) agrupado por años.

Figura 4b. Gráficos de cajas del índice relativo de EVI (RI) agrupado por años.

En general, nuestros resultados son coincidentes con los de Verbesselt et al. (en prensa) en varios aspectos. En particular, para Pinus sylvestris, nuestros resultados muestran que las imágenes MODIS permiten discriminar perfectamente las masas con un grado de decaimiento más severo del resto. El índice NDVI parece detectar mejor estas diferencias que el EVI. Temporalmente, parece que las masas de Pinus sylvestris que actualmente tienen tienen un grado de afección severo y aquellas con un grado de afección leve se distinguen muy bien en cualquier época del año y en todos los años de estudio. Esto puede indicar que: (1) la afección se ha producido muchos años atrás, lo cual es poco probable, ya que los técnicos de campo empezaron a detectar el daño en Filabres en 2002; (2) que las masas tienen características estructurales distintas y por tanto sus condiciones de partida en términos de NDVI son también distintas. Esto se puede investigar a través de análisis correlacionales del grado de afección de las parcelas con distintas variables físicas y estructurales.

Por último, son muy extraños los resultados obtenidos para Pinus nigra. Si bien la respuesta caótica de los índices de vegetación para el grado de afección mayor se puede explicar por la baja representatividad de parcelas en esta clase de daño (n = 1), no encontramos de momento una explicación lógica para la respuesta inversa de las otras dos clases de daño. Es necesario investigar por qué se ha obtenido esta señal. ¿Tal vez la complejidad estructural de los rodales de Pinus nigra es más compleja o heterogénea que los rodales de Pinus sylvestris? ¿Existen factores físicos (p.e. la pendiente) que puedan estar condicionando la respuesta de los rodales al decaimiento y a la vez la respuesta espectral de estas parcelas? Las respuestas a éstas y otras preguntas son críticas para poder determinar el papel de las imágenes MODIS como herramienta de detección de los procesos de decaimiento en plantaciones forestales.