Álgebra de mapas con lenguaje R sin rgdal

El formato ráster se fundamenta en la división del área de estudio en una matriz de celdillas, generalmente cuadradas. Cada una de estas celdillas recibe un único valor que se considera representativo para toda la superficie abarcada por la misma. El nombre genérico de álgebra de mapas se refiere a un conjunto de herramientas de cálculo con matrices de datos. El álgebra de mapas incluye un amplio conjunto de operadores que se ejecutan sobre una o varias capas ráster de entrada para producir una o varias capas ráster de salida. Por operador se entiende un algoritmo que realiza una misma operación en todas las celdillas de una capa raster. Cada capa ráster es una matriz de números y la operación se realiza para todos los números de la matriz, por tanto, para todas las celdillas de la capa raster. En consecuencia, el álgebra de mapas es coincidente con el concepto matemático de álgebra matricial.

Prescindiendo en primera instancia de la complejidad inherente al formato ráster, el cual puede ser cargado directamente al ambiente de R en el formato SpatialGridDataFrame mediante rgdal, vamos a considerar pequeñas matrices de datos para recrear las operaciones de álgebra matricial a pequeña escala y pueda corroborarse fácilmente sus resultados. Por otra parte, los rásteres pueden ser transformados en un SIG al formato Arc/Info ASCII Grid (AAIGrid) y observar, al abrirlo con un editor de texto, que constan de un encabezado de texto seguido del arreglo matricial de filas y columnas. Eliminando el encabezado, se pueden cargar al intérprete de R como data frames y grabar nuevamente como archivos de texto prescindiendo de rgdal. Sin embargo, a estos últimos habría que colocarle manualmente el encabezado para adecuarlo al formato ASCII Grid.

Supongamos que tenemos las siguientes dos matrices en archivos de texto separados, m_sim y m_nosim respectivamente, en el directorio proyectoR:

1    2    3           1    2    3     
2    1    4           4    5    6
3    4    1           7    8    9

Una vez en el intérprete, ejecutamos:

> getwd() #Obtencion del directorio actual de trabajo
[1] "/home/zeito"
> dir() #Listar archivos y directorios del directorio actual de trabajo
 [1] "C++"                            "Clases_T..."               
 [3] "C:\\nppdf32Log\\debuglog.txt"   "datos.out"                     
 [5] "Descargas"                      "Documentos"                    
.
.
.
[25] "proyectoR"                      "Público"                       
[27] "qgis.mp4"                       "qgis.ogv"                      
[29] "R"                              "rad.in"                        
[31] "resultado.txt"                  "salida.out"                    
[33] "sesion1"                        "sesiones_de_R.doc"             
[35] "sesiones_de_R.odt"              "skype_Ubuntu-current_amd64.deb"
[37] "Temario de R.doc"               "Temario de R.odt"              
[39] "Vídeos"                         "Wallpapers"                    
[41] "wordpress.com"                  "wordpress.com~" 

donde se visualiza la existencia de “proyectoR”. Con la siguiente serie de comandos:

> file.info("proyectoR")
          size isdir mode               mtime               ctime
proyectoR 1024  TRUE  775 2014-02-23 00:32:12 2014-02-23 00:32:12
                        atime  uid  gid uname grname
proyectoR 2014-02-23 09:48:43 1000 1000 zeito  zeito
> setwd("proyectoR")
> dir()
 [1] "julio2.in" "julio.in"  "m_nosim"   "m_nosim~"  "m_sim"     "m_sim~"   
 [7] "m_sim_c1"  "prueba"    "sesion1"   "sesion2"  
> file.show('m_sim'); file.show('m_nosim')

se verifica que “proyectoR” es un directorio, lo establecemos como directorio de trabajo, listamos sus archivos y visualizamos el contenido de m_sim y m_nosim antes de cargar los archivos al intérprete.

Con la serie de comandos siguiente cargamos las matrices m_sim y m_nosim y las visualizamos en el intérprete:

> m_sim<-read.table('m_sim')
> m_nosim<-read.table('m_nosim')
> m_sim
  V1 V2 V3
1  1  2  3
2  2  1  4
3  3  4  1
> m_nosim
  V1 V2 V3
1  1  2  3
2  4  5  6
3  7  8  9

Recordemos que el Índice de Vegetación de Diferencias Normalizadas (NDVI por sus siglas en inglés) se puede determinar a partir de las bandas 3 (b3) y 4 (b4) de Landsat (5, 7) mediante la siguiente expresión:

    NDVI = (b4 – b3)/(b4+b3)

Asumamos arbitrariamente que b4 es 2*m_nosim y b3 es m_sim. Entonces:

> b4<- 2*m_nosim
> b3<- m_sim
> NDVI<- (b4-b3)/(b4+b3)
> b4
  V1 V2 V3
1  2  4  6
2  8 10 12
3 14 16 18
> b3
  V1 V2 V3
1  1  2  3
2  2  1  4
3  3  4  1
> NDVI
         V1        V2        V3
1 0.3333333 0.3333333 0.3333333
2 0.6000000 0.8181818 0.5000000
3 0.6470588 0.6000000 0.8947368

donde se pueden corroborar fácilmente los cálculos. Por otra parte, si se quiere determinar los promedios por columnas, por filas y total de los elementos de la matriz, para m_nosim se tiene entonces que:

> m_nosim
  V1 V2 V3
1  1  2  3
2  4  5  6
3  7  8  9
> colMeans(m_nosim)
V1 V2 V3 
 4  5  6 
> rowMeans(m_nosim)
[1] 2 5 8
> sum(m_nosim)/(ncol(m_nosim)*nrow(m_nosim))
[1] 5

Ahora tratemos de cargar al intérprete dos bandas b4 y b3 reales que he transformado en AAIGrid (ascii) mediante QGIS. Aquí tenemos una imagen del contenido de la b3 abierta mediante leafpad:

b3

Guardándolas sin encabezado como b3_se.asc y b4_se.asc, respectivamente, se cargan al intérprete con la siguiente serie de comandos para también determinar el NDVI y grabarlo en ascii como ndvi.asc:

> b3<-read.table('b3_se.asc')
> b4<-read.table('b4_se.asc')
> NDVI<- (b4-b3)/(b4+b3)
> write.table(NDVI, 'ndvi.asc', row.names=FALSE, col.names=FALSE)
> range(NDVI)
[1] -0.6482759  0.9632406

Después de colocarle el encabezado al ndvi.asc, lo abrí en QGIS, le asigné una tabla de colores automática con 5 clases y he aquí el resultado:

ndvi

Es posible entonces hacer álgebra de mapas con R sin necesidad del raster calculator del software SIG. No obstante, con la librería de R rgdal se simplifican las cosas porque no haría falta el concurso del software SIG como intermediario para convertir el ráster en AAIGrid; aunque con gdal_translate de gdal podríamos haber hecho lo mismo en una consola de bash.

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

4 respuestas a Álgebra de mapas con lenguaje R sin rgdal

  1. geoyons dijo:

    Muy bueno. Espero con mucho entusiasmo el curso de R.

  2. Pingback: Álgebra de mapas con lenguaje R sin rgda...

  3. Pingback: Álgebra de mapas con lenguaje R usando la librería rgdal | El Blog de José Guerrero

  4. Pingback: Cómo hacer álgebra de mapas con la librería raster de lenguaje R | El Blog de José Guerrero

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