Cómo determinar la matriz de varianza-covarianza entre bandas ráster mediante lenguaje R

La varianza es una medida de variabilidad o dispersión de los datos. Matemáticamente, corresponde al promedio de las desviaciones al cuadrado de los valores de una serie de datos con relación a la media. Por otra parte, la covarianza es una medida de la extensión en la cual los elementos constituyentes de dos series diferentes de datos se mueven en la misma dirección. Las varianzas y covarianzas frecuentemente se despliegan juntas en una matriz de varianza-covarianza (matriz de covarianza) donde las varianzas aparecen a lo largo de la diagonal principal y las covarianzas fuera de ésta. No obstante, esta matriz es simétrica por lo que los elementos aij son iguales a los elementos aji.

Los detalles algebráicos los pueden encontrar en:

Variance-Covariance Matrix

pero aquí sólo vamos a señalar el procedimiento de adecuación de los datos ráster para determinar este tipo de matrices. Las bandas Landsat 5, por ejemplo, son 7 y puede interesarnos saber si existe asociación o no entre ellas (como en este trabajo). Los valores en cada ráster se encuentran dispuestos a su vez en un arreglo matricial pero nos interesa, para el análisis, desplegarlos como si fuesen vectores columna cuya amplitud será m x n (número de filas por columna). Técnicamente, el procedimiento operativo que se encuentra delineado en el link de arriba consiste en construir esas matrices de 2, 3, 4, 5, 6 o 7 columnas por m x n filas, restarles la media por columna a cada valor en la primera matriz y luego, a la matriz resultado, tomar su transpuesta, multiplicarla por la matriz resultado y, finalmente, dividir cada elemento entre m x n. Si se tomaron 2 bandas el resultado es una matriz simétrica de 2 x 2, para 3 bandas de 3 x 3 y así sucesivamente.

Hacer lo anterior con R es muy sencillo, rápido y eficiente porque este lenguaje, a diferencia de otros, ya tiene implementados por defecto los métodos de álgebra vectorial sin necesidad de usar ningún bucle o estructuras repetitivas (for, while); las cuales por cierto son muy ineficientes en R por lo que se desaconseja su uso (especialmente cuando hay que lidiar con gran volumen de datos como en el caso de rásteres). Para probar el procedimiento voy a usar sólo 2 rásteres, de 2181 filas y 2668 columnas (casi 6 millones de píxeles cada uno), que corresponden a las bandas 1 (blue) y 6 (térmica) de Landsat 5 para cierta zona de Venezuela. El compendio de pasos para obtener la matriz de covarianza es el siguiente:

>setwd('proyectoR') #establece el directorio de trabajo
>list.files(pattern='b') #se verifica la existencia de los rásteres
>library(rgdal) #carga la libreria rgdal
>#se leen ambos raster
>b1<-readGDAL('b1_clip.tif')
b1_clip.tif has GDAL driver GTiff 
and has 2181 rows and 2668 columns
>b6<-readGDAL('b6_clip.tif')
b6_clip.tif has GDAL driver GTiff 
and has 2181 rows and 2668 columns
>#se asignan sus values a dos vectores
>z1<-b1$band1
>z6<-b6$band1
>length(z1)
[1] 5818908
>length(z6)
[1] 5818908
> range(z1)
[1]  34 255
> 255-34
[1] 221
> range(z6)
[1] 106 152
> 152-106
[1] 46
>a<-cbind(z1,z6) #se combinan los values para formar una matriz
>#se determinan las medias de los values
> meanz1<-mean(z1)
> meanz6<-mean(z6)
> meanz1
[1] 64.7216
> meanz6
[1] 128.1068
>v1<-matrix(meanz1,5818908,1)
>v6<-matrix(meanz6,5818908,1)
>b<-cbind(v1,v6) #se crea la matriz de medias de 5818908 x 2
>c<- a-b #matriz de valores menos matriz de medias
>d<- t(c)%*%c 
>d
           z1         z6
z1 2089897502 -234345982
z6 -234345982   65189791
e<-d/5818908 #matriz de covarianza
e
          z1        z6
z1 359.15631 -40.27319
z6 -40.27319  11.20310

A continuación tenemos los histogramas de las variables z1 y z6:

hist(z1, col='green', border='blue')
hist(z6, col='blue', border='green')

hist

y previo al análisis de componentes principales, con estos comandos:

> auto<-eigen(e)
> auto
$values
[1] 363.756824   6.602579

$vectors
           [,1]       [,2]
[1,] -0.9935386 -0.1134947
[2,]  0.1134947 -0.9935386

ya se tienen los autovalores y autovectores de la matriz de covarianza.

Cada raster considerado tiene casi 6 millones de píxeles y las operaciones matriciales con R son prácticamente instantáneas. Por cierto, como R también es un lenguaje de scripts todo lo anterior puede automatizarse de esa manera.

Esta entrada fue publicada en Lenguaje R, SIG, Software Libre. Guarda el enlace permanente.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s