Modelos binomiales en R

Muchos problemas estadísticos implican el análisis de una respuesta binaria. Por ejemplo, muchas veces se clasifican los individuos muestrales como:

• Vivos o muertos,
• Presencia o ausencia de especies,
• Hombre o mujer,
• Sano o enfermo,
• Deforestado o no deforestado.

En estos y otros casos, los modelos binarios suponen una herramienta de análisis muy interesante. Para ello es conveniente que, al menos, una de las variables explicativas sea continua. En el caso de la deforestación (o cualquier otro cambio de uso binario), los datos consistirán, por ejemplo, en una lista de variables (ej. distancia a caminos, elevación, pendiente) medidas en píxeles que han sido deforestados en un periodo determinado, y una lista similar para píxeles que no han sido deforestados en ese mismo periodo. La cuestión es tratar de averiguar cuál de estas variables, si es que hay alguna, aumentan la probabilidad de que un píxel individual sea deforestado para ese periodo.

La variable respuesta contiene sólo 0s y 1s. Por ejemplo, 0 representaría píxeles que permanecen como bosque y 1 píxeles que han sido deforestados. Por tanto, habría una única columna para la variable respuesta. Los pasos a seguir para modelar una variable respuesta de tipo binomial son los siguientes:

1. Crear un vector de 0s y 1s como variable respuesta.
2. Ajustar un modelo lineal generalizado (funcion glm) con una familia de distribución de errores de tipo binomial (argumento family = binomial).
3. Reducir el modelo eliminando los términos que no sean significativos y comparando el cambio en devianza con respecto al modelo total.
4. Interpretación de los resultados.
   

La función de vínculo en los glm hace referencia a la transformación de la(s) variable(s) antes de ajustar el modelo. En el caso de los modelos binomiales, la función de vínculo que se usa habitualmente es la logística (logit), que transforma con logaritmos naturales la variable respuesta. Por defecto, no hace falta especificar ésto en los modelos.

A la hora de comparar modelos se utiliza la devianza como un indicador de la bonanza del modelo. Para aquellos que no estén muy familiarizados con la jerga estadística, la devianza es una medida del ajuste del modelo en relación al número de parámetros del mismo. Si un modelo explica lo mismo que otro, pero es más sencillo (menor número de párametros) es, por regla general, mejor modelo. Y esto es lo que va a indicar la devianza.

Un ejemplo: El caso de la deforestación en Chile Central

En este ejemplo, la variable respuesta va a ser la ocurrencia de deforestación en Chile Central para el periodo 1975-1985, en dónde 1 = deforestación y 0 = no deforestación. Para obtener estos datos, se llevó a cabo un muestreo sistemático cada 1 km en el área de estudio y se registró la cobertura en 1975 y 1985 para todos aquellos píxeles que eran bosque en 1975. Se midieron también toda una serie de variables explicativas que se detallan a continuación:

• Elevación (m),
• Pendiente (grados),
• Orientación (grados),
• Distancia a ríos (m),
• Distancia a ciudades grandes (m),
• Distancia a pueblos (m),
• Distancia a carreteras principales (m),
• Distancia a carreteras secundarias (m),
• Distancia a suelo agrícola (m),
  
• Distancia al borde del bosque (m).

Los datos originales están en Excel. Para leerlos directamente en R suele ser más rápido pasarlos a *.txt, aunque R puede leer también datos en Excel. Los datos están accesibles en una página web, por lo que se pueden leer desde R utilizando el comando url(). Basta con copiar el código que se ve en azul y pegarlo en R. Funcionará.

fnf75_85 <- read.table(url("http://archivos-para-subir.googlegroups.com/web/F-NF1975-85.txt?gda=f500bkEAAACkOfwqHAVd4YqgfIB09GDRCcOUO-z8uvHRarwVOPM5xOEDwyLfg5d1mmtcx5AUJ3JTCT_pCLcFTwcI3Sro5jAzlXFeCn-cdYleF-vtiGpWAA"), header = T, dec = ".", sep = "\t")

Ya se ha creado un objeto en R llamado fnf75_85, que contiene los datos de las distintas variables que se utilizarán en el análisis. Con el argumento header = T se especifica que la primera línea contiene los nombres de las variables. Con el argumento dec = "." se dice cómo vienen especificados los decimales (alternativamente podría ser ","). Por último, el argumento sep = "\t" especifica que los campos están separados por tabulaciones. Con la función str() se obtiene la estructura del objeto fnf75_85 e información sobre la naturaleza (integer, factor, character, etc) de cada una de las variables.

str(fnf75_85)

'data.frame': 1980 obs. of 15 variables:
$ ID_ALL : int 10 20 77 80 81 118 158 208 255 262 ...
$ ID_FOREST : int 1 2 3 4 5 6 7 8 9 10 ...
$ Class1975 : int 1 1 1 1 1 1 1 1 1 1 ...
$ Class1985 : Factor w/ 8 levels "1","2","3","5",..: 4 3 2 2 2 3 1 1 2 1 ...
$ Change : Factor w/ 3 levels "0","1","NoData": 2 2 2 2 2 2 1 1 2 1 ...
$ Elevation : num 435 136 299 300 300 ...
$ Slope : num 19.6 0.4 0.2 0 0 0 0.2 12.3 0.3 3.8 ...
$ Aspect : num 149 290 357 -1 -1 ...
$ Dist_river : num 4536 1025 751 2005 2672 ...
$ Dist_City.20T: num 27750 52254 28529 25994 25180 ...
$ Dist_village : num 2271 2620 1142 1084 742 ...
$ Dist_road_P : num 4601 2313 4621 1973 1442 ...
$ Dist_road_S : num 752 395 323 552 0 ...
$ Dist_agri : num 94.9 60.0 30.0 30.0 30.0 ...
$ Dist_f_forest: num 84.9 524.8 67.1 161.6 67.1 ...

La variable respuesta es $Change. Se observa que es un factor con 3 niveles: "0", "1" y "NoData". No interesa tener muestras con "NoData", por lo que se deben de eliminar estos datos. Se puede hacer sobreescribiendo el archivo fnf75_85 (! indica exclusión de casos).

fnf75_85 <- fnf75_85[!fnf75_85$Change == "NoData", ]
fnf75_85$Change <- as.integer(fnf75_85$Change)
fnf75_85$Change <- ifelse(fnf75_85$Change == 1, 0, 1)

Y ahora un resumen de la variable respuesta.

table(fnf75_85$Change)

0 1
774 1184

Así que hay 774 píxeles que continúan siendo bosque en 1985 y 1184 píxeles que han sido deforestados entre 1975 y 1985. El diseño es claramente no balanceado (distinto número de casos en la variable respuesta), aunque esto no tiene por qué influir necesariamente en el ajuste del modelo. El ajuste del modelo completo (excluyendo la variable Dist_f_forest, que ha sido mal calculada) se haría de la siguiente forma:

model.full <- glm(Change ~ Elevation + Slope +Aspect + Dist_river + Dist_City.20T + Dist_village + Dist_road_P + Dist_road_S + Dist_agri, data = fnf75_85, family = binomial)

Para explorar los resultados se puede usar la función summary().

summary(model.full)

Call:
glm(formula = Change ~ Elevation + Slope + Aspect + Dist_river +
Dist_City.20T + Dist_village + Dist_road_P + Dist_road_S +
Dist_agri, family = binomial, data = fnf75_85)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.5265 -1.3419 0.9214 1.0039 1.3176

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 6.177e-01 1.650e-01 3.742 0.000182 ***
Elevation 3.529e-04 2.209e-04 1.597 0.110201
Slope 3.913e-03 4.902e-03 0.798 0.424748
Aspect -8.098e-04 4.641e-04 -1.745 0.081044 .
Dist_river 1.350e-05 6.302e-05 0.214 0.830348
Dist_City.20T -3.136e-06 4.025e-06 -0.779 0.435960
Dist_village 2.353e-05 3.188e-05 0.738 0.460615
Dist_road_P -2.579e-05 9.662e-06 -2.669 0.007608 **
Dist_road_S -8.984e-05 6.592e-05 -1.363 0.172909
Dist_agri -4.340e-05 5.029e-05 -0.863 0.388227
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 2627.9 on 1957 degrees of freedom
Residual deviance: 2607.7 on 1948 degrees of freedom
AIC: 2627.7

Number of Fisher Scoring iterations: 4

Cómo se observa, las variables que influyen significativamente sobre la deforestación en el periodo 1975-1985 son la distancia a carreteras principales (Dist_road_P) y, de manera más marginal (p = 0.081), la orientación (Aspect). La primera variable influye negativamente sobre la probabilidad de deforestacion (Estimate = -2.579 e-05), es decir, que a medida que nos alejamos de las carreteras principales la probabilidad de deforestación disminuye. La segunda es más difícil de intepretar, puesto que la orientación va de 0º (Norte) a 360º (Norte). En este caso, sería más interpretable utilizar una variable como insolación, que está directamente relacionada con la orientación y la pendiente.

Puesto que el resto de variables no son significativas se puede reducir el modelo.

model.red1 <- glm(Change ~ Aspect + Dist_road_P, data = fnf75_85, family = binomial)
summary(model.red1)

Call: glm(formula = Change ~ Aspect + Dist_road_P, family = binomial, data = fnf75_85) Deviance Residuals: Min 1Q Median 3Q Max -1.5055 -1.3468 0.9385 1.0064 1.2628 Coefficients: Estimate Std. Error z value Pr(>|z|)

(Intercept) 7.514e-01 1.092e-01 6.881 5.92e-12 ***
Aspect -7.502e-04 4.593e-04 -1.633 0.10236
Dist_road_P -2.496e-05 8.126e-06 -3.072 0.00213 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 2627.9 on 1957 degrees of freedom
Residual deviance: 2615.3 on 1955 degrees of freedom
AIC: 2621.3

Number of Fisher Scoring iterations: 4

Y se comparan los dos modelos con la función anova().

anova(model.full, model.red1, test = "Chi")

Analysis of Deviance Table

Model 1: Change ~ Elevation + Slope + Aspect + Dist_river + Dist_City.20T +
Dist_village + Dist_road_P + Dist_road_S + Dist_agri
Model 2: Change ~ Aspect + Dist_road_P
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 1948 2607.66
2 1955 2615.32 -7 -7.66 0.36

Lo que indica que el modelo más simple no es significativamente peor que el modelo completo (p = 0.36). Así que aceptamos este modelo por el momento. Sin embargo, cuando analizamos los resultados se observa que la variable orientación no es significativa (p = 0.102), por lo que probamos a reducir el modelo un poco más.

model.red2 <- glm(Change ~ Dist_road_P, data = fnf75_85, family = binomial)
summary(model.red2)

Call: glm(formula = Change ~ Dist_road_P, family = binomial, data = fnf75_85) Deviance Residuals: Min 1Q Median 3Q Max -1.4540 -1.3500 0.9432 1.0034 1.2042 Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) 6.333e-01 8.129e-02 7.791 6.66e-15 ***
Dist_road_P -2.551e-05 8.115e-06 -3.144 0.00167 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 2627.9 on 1957 degrees of freedom
Residual deviance: 2618.0 on 1956 degrees of freedom
AIC: 2622

Number of Fisher Scoring iterations: 4

Y volvemos a comparar los dos modelos:

anova(model.red1, model.red2, test = "Chi")

Analysis of Deviance Table

Model 1: Change ~ Aspect + Dist_road_P
Model 2: Change ~ Dist_road_P
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 1955 2615.32
2 1956 2617.99 -1 -2.67 0.10

Así que de nuevo nos encontramos que el modelo reducido (model.red2) no es significativamente peor que el modelo anterior (model.red1). Aceptamos este último como válido, por lo que concluimos que la probabilidad de deforestación en el periodo 1975-1985 estuvo influida fundamentalmente por la presencia de carreteras principales (¿se debió la deforestación a un uso maderero por parte de empresas forestales?).

Por último, se representan gráficamente los resultados.

xv<-seq(0,30000,10)
yv<-predict(model.red2, list(Dist_road_P = xv), type = "response")
plot(as.integer(fnf75_85$Change) ~ fnf75_85$Dist_road_P, type = "p", ylab = "Probability of forest loss", xlab = "Distance to main roads (m)")
lines(xv,yv)

¿Qué nos dice este gráfico? Básicamente, que la probabilidad de deforestación disminuye con la distancia a carretera principales, pero disminuye poco. Así, cuando uno se aleja 25 kms de las carreteras principales (una distancia bastante grande), la probabilidad de deforestación disminuye del 65% a apenas el 50%.

Cómo hay muchos puntos en el gráfico, no se aprecia realmente si hay más 0s que 1s en distancias próximas a las carreteras principales y a la inversa... sería bueno intentar ajustar funciones de densidad de los 0s y 1s en el mismo gráfico.

¿Qué ocurre si la variable respuesta es una proporción?

Puede ser que la variable respuesta no sean 0s y 1s sino que sea una proporción o porcentaje, ej. porcentaje de semillas que sobreviven en un experimento, sex ratios, proporción de píxeles deforestados en un determinado área, etc. En este caso, tambińe podemos llevar a cabo un análisis glm con distribución de errrores de tipo binomial. Las únicas diferencias con respecto a los pasos descritos anteriormente serían:

• Para definir la variable respuesta creamos un objeto con dos columnas, utilizando la función cbind() para juntar dos vectores que contengan el número de éxitos y fracasos en cada uno de los casos.

n.fracasos <- denominador - n.exitos
var.respuesta <- cbind(n.exitos, n.fracasos)

• Comprobar que no hay sobredispersión (la devianza residual no es mayor que los grados de libertad residuales) y corregirlo en caso necesario utilizando una familia de tipo cuasibinomial (family = quasibinomial).
  

Bibliografía recomendada:

- Crawley, M.J. 2007. The R Book. Wiley, London.

Comparación de la pendiente de dos modelos (normal y poisson) de riqueza de especies

En un estudio que está llevando a cabo Fabio Suzart de Alburquerque (Departamento de Ecología, Universidad de Alcalá) se quiere comparar la pendiente de la recta de regresión de la riqueza de especies nativas y exóticas en el Reino Unido frente a una variable explicativa relacionada con la entrada de energía en el sistema (el primer eje de un análisis de componente principales), correlacionado fundamentalmente con la temperatura media anual, la temperatura media del mes de enero, y la evapotranspiración potencial). El principal problema con el que se enfrenta es que la distribución de la riqueza de nativas sigue una distribución normal, mientras que la distribución de especies exóticas sigue una distribución más de tipo Poisson o log-normal (ver Figura 1).

El siguiente código de R permite la lectura de los datos (accesibles con conexión a internet), el ajuste de dos modelos "normales" para la riqueza de nativas y exóticas, y la creación de: (1) histogramas de los residuos de los modelos de riqueza de especies nativas y exóticas; y (2) gráficos de correlación de estas variables con el primer eje del análisis de componentes principales (PCA).

data <- read.table(url("http://archivos-para-subir.googlegroups.com/web/datos+regresion.txt?gda=x8KpMUUAAACkOfwqHAVd4YqgfIB09GDRRDpDrIn9MKnJm3pWH7Uz3yEZ0eBGq9HKqr7NTMdj8uts7Hm76flVp_vYI1lbK8kaGu1iLHeqhw4ZZRj3RjJ_-A"), header = T, sep = "\t")

glm.native <- glm(data$NAT_RICH ~ data$Factor1.energy) glm.exotic <- glm(data$ALIEN_RICH ~ data$Factor1.energy)
par(mfcol= c(2, 2))
hist(residuals(glm.native), main = "Histogram of native SR residuals", xlab = "Number of species")
par(new = T)
plot(density(residuals(glm.native)), xlab = "", ylab = "", main = "", axes = F)
hist(residuals(glm.exotic), main = "Histogram of exotic SR residuals", xlab = "Number of species")
par(new = T)
plot(density(residuals(glm.exotic)), xlab = "", ylab = "", main = "", axes = F)
plot(data$NAT_RICH ~ data$Factor1.energy, main = "Native species richness vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "Number of species")
plot(data$ALIEN_RICH ~ data$Factor1.energy, main = "Exotic species richness vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "Number of species")

Figura 1. Histogramas de los residuos de los modelos de especies nativas (arriba izquierda) y exóticas (abajo izquierda) y gráficos de correlación de la riqueza de especies nativas (arriba derecha) y exóticas (abajo derecha) con el primer eje del análisis de componentes principales (PCA).

Una posible solución a este problema es utilizar modelos lineales generalizados (GLM) con una distribución de errores de tipo Poisson y una función de vínculo de tipo logarítmica para linealizar la relación entre la riqueza de especies nativas y la variable independiente. Esto impide comparar la pendiente de ambas rectas de regresión de manera directa pero nos permite, en cambio, establecer hipótesis sobre los procesos que generan la respuesta de las variables dependientes (riqueza de nativas frente a riqueza de exóticas).

glm.native <- glm(data$NAT_RICH ~ data$Factor1.energy)
glm.exotic <- glm(data$ALIEN_RICH ~ data$Factor1.energy, family = poisson)

El primer modelo es un modelo normal. Hay que considerar que los GLM no son más que una extensión de los modelos lineales, por lo que los modelos lineales que comúnmente se utilizan (ANOVA, regresión) no son más que un caso particular de GLM en dónde la distribución de errores es de tipo gaussiana (normal). Con la función summary() se obtiene un resumen de los modelos ajustados.

summary(glm.native)

Call:
glm(formula = data$NAT_RICH ~ data$Factor1.energy)

Deviance Residuals:
Min 1Q Median 3Q Max
-428.470 -45.604 5.102 47.498 218.190

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 421.489 1.515 278.12 <2e-16>
data$Factor1.energy 59.039 1.516 38.95 <2e-16>
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 5172.153)

Null deviance: 19483365 on 2251 degrees of freedom
Residual deviance: 11637345 on 2250 degrees of freedom
AIC: 25652

Number of Fisher Scoring iterations: 2

summary(glm.exotic)

Call:
glm(formula = data$ALIEN_RICH ~ data$Factor1.energy, family = poisson)

Deviance Residuals:
Min 1Q Median 3Q Max
-18.7702 -4.5256 -0.3740 3.1827 23.9153

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.900846 0.001945 2519.6 <2e-16 ***
data$Factor1.energy 0.543437 0.002224 244.3 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 158010 on 2251 degrees of freedom
Residual deviance: 87121 on 2250 degrees of freedom
AIC: 101806

Number of Fisher Scoring iterations: 5

Como se puede observar, las pendientes de la recta de regresión de ambos modelos no son comparables, entre otras cosas, porque en el modelo Poisson la variable respuesta (riqueza de exóticas) está transformada logarítmicamente. Esto ya nos da una idea sobre el comportamiento de ambos tipos de riqueza de especies. Así, mientras que la riqueza de especies nativas varía linealmente con el primer eje del PCA (energía), la riqueza de exóticas cambia exponencialmente, es decir, que para valores bajos de la variable explicativa, el incremento de la riqueza será relativamente pequeño, pero para valores altos de la variable explicativa, el incremento de la riqueza de especies será muy alto.

Esto se puede observar en el modelo ajustado a los puntos para cada uno de los dos tipos de riqueza.

par(mfrow=c(1,2))
yv.native <- predict(glm.native, type= "response")
plot(data$NAT_RICH ~ data$Factor1.energy, main = "Native species richness vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "Number of species")
lines(sort(data$Factor1.energy), sort(yv.native))
yv.exotic <- predict(glm.exotic, type= "response")
plot(data$ALIEN_RICH ~ data$Factor1.energy, main = "Exotic species richness (response) vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "Number of species")
lines(sort(data$Factor1.energy), sort(yv.exotic))


Y alternativamente podemos dibujar las predicciones del modelo en escala logarítmica. Esta relación es lineal, pero hay que tener en cuenta que la variable respuesta (riqueza de especies) está en escala logarítmica, por lo que la comparación de las pendientes de ambas rectas no es posible.

par(mfrow=c(1,2))
yv.native <- predict(glm.native, type= "response")
plot(data$NAT_RICH ~ data$Factor1.energy, main = "Native species richness vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "Number of species")
lines(sort(data$Factor1.energy), sort(yv.native))
yv.exotic2 <- predict(glm.exotic, type= "link")
plot(log(data$ALIEN_RICH) ~ data$Factor1.energy, main = "Exotic species richness (link) vs. PCA 1", xlab = "PCA factor 1 (energy)", ylab = "log(Number of species)")
lines(sort(data$Factor1.energy), sort(yv.exotic2))

Este estudio forma parte del siguiente trabajo: Suzart de Alburquerque, F., Castro, P. Cayuela, L. & Rodríguez, M.A. Assessing the influence of environmental and human-related factors on native and exotic species richness in Great Britain. In prep.

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.

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.

Análisis espacio-temporal de plagas de procesionaria en Andalucía

La procesionaria del pino es una de las plagas más importantes de los pinares mediterráneos. Debe su nombre de "Procesionaria" a que se desplaza en grupo de forma alineada, a modo de procesión. En Andalucía, la procesionaria es un problema importante. Durante los últimos años, la Consejería de Medio Ambiente de la Junta de Andalucía ha estado recopilando información sobre el grado de afección de plantaciones y pinares autóctonos en toda Andalucía.

Aunque los datos que se muestran aquí son datos todavía incompletos, el objetivo es mostrar el potencial de R para el manejo, depuración y análisis de bases de datos extensas. R no sólo permite todo ésto, sino que además ofrece la posibilidad de trabajar con datos espaciales simulando las capacidades y potencialidades de SIG como ArcGIS o Idrisi.

Los datos están accesibles a través de un ftp, por lo que todo el código que se muestra a continuación se puede reproducir en R desde cualquier ordenador. El código se muestra en azul y letra Courier.

Empecemos cargando (e instalando si fuera esta la primera sesión) los paquetes con los que vamos a trabajar:

install.packages(c("sp", "maps", "maptools", "rgdal", "gpclib", "spatstat"))

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

Leemos ahora los datos de daños (en formato *.txt) por procesionaria en rodales forestales de Andalucía. Los datos, como ya dije, están en un ftp, por lo que podemos leerlos directamente desde R y salvarlos como un objeto de R al que llamaremos daño. Si los puertos de acceso no están habilitados puede ser que este código no funcione y no se pueda leer el archivo.

daño <- read.table(url( "ftp://fundacionfire.org:Fire1234@fundacionfire.org/data/Procesionaria/daños.txt"), header = T, sep = "\t")
str(daño)

'data.frame': 30550 obs. of 3 variables:
$ Rodal: Factor w/ 3382 levels "AL-000-001","AL-000-002",..
$ Fecha: Factor w/ 675 levels "10/20/2003","1/11/2003",..
$ Grado: int 0 0 0 0 0 0 0 0 0 0 ...

Es decir, tenemos un arreglo de datos (data.frame) con 30.550 observaciones que contienen información sobre: (1) el rodal (identificador del rodal en dónde se tomo el dato de campo), (2) la fecha (dd/mm/aa), y (3) el grado de afección (de 0 a 5, siendo 0 = ningún daño y 5 = daño severo). El primer paso que vamos a seguir es poner todas las fechas en año y eliminar toda la información referida al mes y al día en el que se hizo la observación. De paso, haremos algunas correcciones de años que están mal escritos (p.e. 0002 en vez de 2002).

daño$Fecha <- do.call("rbind", strsplit(as.character(daño$Fecha), split = "/"))[,3]
daño$Fecha <- sub("0002", "2002", daño$Fecha)
daño$Fecha <- sub("0202", "2002", daño$Fecha)
daño$Fecha <- as.factor(daño$Fecha)

Ahora exploramos visualmente los datos.

plot(table(table(daño$Rodal)), xlab = "Número de datos por rodal", ylab = "Número de rodales", cex.lab = 1.3, cex.axis = 1.3)

Esta gráfica muestra el número de datos que tenemos por rodal. Como sólo hay una observación por año en cada rodal y tenemos información para el periodo 1992-2003 (ambos años inclusive), supuestamente deberíamos tener un máximo de 12 datos por rodal. Sin embargo vemos dos cosas en esta gráfica: (1) que muchos rodales no tienen información para el periodo completo, y (2) que para algunos rodales tenemos más datos de los esperados, llegando incluso a 81 datos por rodal. Esto indica que o bien los datos están repetidos o hay más información de la que cabría esperar. Vamos a ver si hay elementos duplicados y, en caso de que los haya, vamos a eliminarlos.

rodal.sel <- table(daño$Rodal)[table(daño$Rodal)>12]
daño <- data.frame(daño, dup = paste(daño$Rodal, daño$Fecha, sep = ""))
daño <- daño[!duplicated(daño$dup), ]
daño <- daño[ , -4]
str(daño)

'data.frame': 30243 obs. of 3 variables:
$ Rodal: Factor w/ 3382 levels "AL-000-001","AL-000-002",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Fecha: Factor w/ 12 levels "1992","1993",..: 12 11 10 9 8 7 6 5 4 3 ...
$ Grado: int 0 0 0 0 0 0 0 0 0 0 ...

Como se puede observar, ahora tenemos 30243 datos, frente a los 30550 que teníamos al principio. La diferencia entre estas dos cifras es la cantidad de elementos que estaban repetidos en la base de datos.

Lo que vamos a hacer ahora es separar los datos por años, de tal forma que tengamos una columna con el identificador del rodal y una columna por cada uno de los años para los cuales disponemos de datos. Esto podría ser tremendamente complicado con Excel ya que la información es altamente heterogénea y en cada año tenemos información sólo para unos rodales determinados. En R implica el uso de la función merge.

pt <- split(daño, daño$Fecha)
a <- merge(pt[[1]], pt[[2]], by = "Rodal", all = T)
daño92_93 <- data.frame(Rodal = a$Rodal, Daño.1992 = a$Grado.x, Daño.1993 = a$Grado.y)
b <- merge(pt[[3]], pt[[4]], by = "Rodal", all = T)
daño94_95 <- data.frame(Rodal = b$Rodal, Daño.1994 = b$Grado.x, Daño.1995 = b$Grado.y)
c <- merge(pt[[5]], pt[[6]], by = "Rodal", all = T)
daño96_97 <- data.frame(Rodal = c$Rodal, Daño.1996 = c$Grado.x, Daño.1997 = c$Grado.y)
d <- merge(pt[[7]], pt[[8]], by = "Rodal", all = T)
daño98_99 <- data.frame(Rodal = d$Rodal, Daño.1998 = d$Grado.x, Daño.1999 = d$Grado.y)
e <- merge(pt[[9]], pt[[10]], by = "Rodal", all = T)
daño00_01 <- data.frame(Rodal = e$Rodal, Daño.2000 = e$Grado.x, Daño.2001 = e$Grado.y)
f <- merge(pt[[11]], pt[[12]], by = "Rodal", all = T)
daño02_03 <- data.frame(Rodal = f$Rodal, Daño.2002 = f$Grado.x, Daño.2003 = f$Grado.y)
daño <- merge(daño92_93, daño94_95, by = "Rodal", all = T) daño <- merge(daño, daño96_97, by = "Rodal", all = T)
daño <- merge(daño, daño98_99, by = "Rodal", all = T)
daño <- merge(daño, daño00_01, by = "Rodal", all = T)
daño <- merge(daño, daño02_03, by = "Rodal", all = T) str(daño)

'data.frame': 3382 obs. of 13 variables:
$ Rodal : Factor w/ 3382 levels "AL-000-001","AL-000-002"
$ Daño.1992: int 0 0 NA NA NA NA NA NA NA NA ...
$ Daño.1993: int 0 0 NA NA NA NA NA NA NA NA ...
$ Daño.1994: int 0 0 0 0 0 3 0 NA 2 0 ...
$ Daño.1995: int 0 0 0 0 0 1 0 NA 1 0 ...
$ Daño.1996: int 0 0 0 0 1 0 0 NA 0 0 ...
$ Daño.1997: int 0 0 0 0 2 1 0 NA 1 0 ...
$ Daño.1998: int 0 0 1 0 1 3 3 NA 2 1 ...
$ Daño.1999: int 0 0 0 0 0 1 0 NA 0 1 ...
$ Daño.2000: int 0 0 NA NA 0 1 0 NA 0 3 ...
$ Daño.2001: int 0 0 0 0 0 0 0 NA 1 0 ...
$ Daño.2002: int 0 0 0 0 0 0 0 0 0 0 ...
$ Daño.2003: int 0 0 0 0 0 0 0 0 0 0 ...

Lo que tenemos ahora es un arreglo de datos con 13 variables, la primera correspondiente al identificador y las doce siguientes haciendo referencia al año para el cual se tienen datos de daños. Hay un total de 3382 rodales para los cuales tenemos información. Sin embargo, como vimos antes, no para todos los rodales tenemos información de cada uno de los doce años del período 1992-2003. Vamos a ver cuántos rodales tienen valores faltantes (NA).

f1 <- function(x) sum(is.na(x))
table(apply(daño, 1, f1))

0 1 2 3 4 5 6 7 8 9 10 11
313 216 1607 528 207 72 35 20 24 148 67 145

Es decir, que sólo tenemos información completa para todo el periodo para un total de 313 rodales. En 216 sólo falta un dato. En 1607 falta información de dos fechas. En 528 falta información de tres fechas, y así sucesivamente. Una vez más podemos observar que la información es bastante heterogénea.

El siguiente paso es asociar esta información a una capa vectorial de polígonos (formato shapefile *.shp) que representan los rodales de pinar en Andalucía. De esta manera podremos representar los datos en mapas. Para simplificar la lectura de datos, he incorporado esta capa y otra capa con las provincias de Andalucía en un objeto de R (*.Rob) que he subido al ftp.

load(url("ftp://fundacionfire.org:Fire1234@fundacionfire.org/data/Procesionaria/mapas.Rob"))

Esto importará a nuestra sesión de R dos objetos espaciales (SpatialPolygonsDataFrame) llamados mapa.rodales y andalucia. Lo primero que tenemos que ver es si los identificadores de los nombres de los rodales son iguales en la capa vectorial (mapa.rodales) y en nuestra base de datos (daño).

match(mapa.rodales@data$ID, daño$Rodal)

Esto indica que no hay ninguna coincidencia. Si miramos más en detalle los nombres de los rodales en la base de datos, vemos que hay diferencias con respecto a la forma en la que se llaman en la capa vectorial. Así un rodal llamado AL-000-001 en la base de datos se denomina AL000001 en la capa vectorial. Vamos a solucionar ésto:

daño$Rodal <- gsub("-", "", daño$Rodal)

Y ahora vamos a ver cuántos polígonos en la capa vectorial no tienen correspondencia con ninguno de los rodales en la base de datos y vice-versa.

sum(is.na(match(mapa.rodales@data$ID, daño$Rodal)))
[1] 625

sum(is.na(match(daño$Rodal, mapa.rodales@data$ID)))
[1] 59

Así que hay 625 polígonos para los cuáles no tenemos ni un sólo dato. ¿A qué es debido ésto? ¿Corresponden a otro tipo de vegetación? Igualmente hay 59 rodales en la base de datos que no tienen correspondencia con ninguno de los polígonos en la capa vectorial. Estos últimos datos no podremos usarlos porque no tienen representación espacial.

Lo siguiente es incorporar la información de la base de datos a la capa vectorial.

colnames(mapa.rodales@data)[1] <- "Rodal"
mapa.rodales@data <- merge(mapa.rodales@data, daño, by = "Rodal", all.x = T)

Y ahora ya sólo queda representar todos estos datos en mapas que generaremos para cada uno de los años del período de estudio. Es importante que antes de copiar este código establezcamos el directorio de trabajo. En Windows se puede hacer directamente desde Archivo (File). En Linux se puede hacer mediante el comando setwd().

mapa.rodales1 <- mapa.rodales
mapa.rodales1@data[, c(9:length(mapa.rodales1@data[1, ]))] <- mapa.rodales1@data[, c(9:length(mapa.rodales1@data[1, ]))] + 1
cols <- c("yellow", "orange1", "red1", "red2", "red3", "red4")
png(filename = "animation%02d.png", width = 800, height = 500)
for (i in 9:length(mapa.rodales@data[1, ])) {
plot(andalucia, axes = F, xlim = c(50000, 610000))
rrt <- mapa.rodales1@data[, i] plot(mapa.rodales, col = cols[rrt], border="transparent", las=1, add = T)
plot(andalucia, add = T)
axis(1, at = c(100000 + 0:4 * 125000), cex.axis = 1.2)
axis(2, at = c(3900000 + 0:4*100000), cex.axis = 1.2)
box()
legend("topleft", legend = c("Sin daño", "Daño bajo", "Daño medio bajo", "Daño medio alto", "Daño alto", "Daño severo"), fill = c("yellow", "orange1", "red1", "red2", "red3", "red4"), bty = "n", title = "Nivel de daño")
title(main = paste("Año", strsplit(colnames(mapa.rodales@data[i]), "\\.")[[1]][2]), cex.main = 2)
}
dev.off()

Esto va a crear doce mapas como el que se muestra a continuación.

¡Y ahora la parte más interesante! Vamos a crear una animación con estas imágenes utilizando una aplicación de Java creada por Barry Rowlingson (http://www.maths.lancs.ac.uk/~rowlings/Chicas/DIY/), de la University of Lancaster.

El resultado final se puede observar en http://rsgis.ugr.es/sacacorchos/Animation/sample.html. Se puede regular la velocidad con la que se suceden las imágenes ... ¡Fantástico! ¿no?