R + SAGA Geoprocesamiento con RSAGA Package

En la presente entrada trataré un tema un poco diferente a lo acostumbrado, sobre todo porque será la primera vez que me refiero al Software R, al que de un tiempo muy reciente lo estoy practicando; por otro lado, aunque tampoco me acuerdo haber discutido sobre SAGA GIS específicamente, es un software del cual ya tengo muchos años trabajándolo, y en donde me siento más cómodo, sobre todo porque se complementa muy bien con el QGIS.

Bueno, en esta aventura voy a mostrar mi experiencia probando el package "RSAGA", por el que me animé luego de leer este enlace, a partir de ello, estuve revisando otras fuentes y decidí intentarlo empleando solamente una capa raster de elevación, el cual será nuestro punto de partida.

Fuente de Datos:

Para esta oportunidad como lo indicaba se empleará una capa raster tipo DEM, el cual está disponible a través del Modelo Digital de Superficie ALOS World 3D -30m (AW3D30), solamente se tienen que registrar y podrán descargar los datos, para el presente ejemplo se seleccionó la zona que comprende el Departamento de Lima-Perú.

Figura 1: DEM (30 m) de la cuenca del río Rimac 

Pasos Previos

Debemos tener en cuenta que el package "RSAGA", al parecer por temas de actualización no funciona con las últimas versiones de SAGA, por lo tanto, es recomendable instalar las versiones 2.14 o  2.2. En relación a R, se recomienda que instalen el RStudio, sobre todo para tener un entorno que nos facilite nuestro trabajo, como por ejemplo la instalación de otros paquetes adicionales requeridos.

Manos a la Obra

Paso 1: Como primer paso, una vez dentro de RStudio, debemos definir el directorio de trabajo, para ello, lo podemos hacer desde la consola con  setwd("Ruta del Directorio"), o también te puedes activar el menú de "Session", y luego dirigirte a "Choose Directory", tal como se aprecia en la figura 2.

Figura 2: Definiendo nuestro directorio de trabajo

Paso 2: Instalar y Activar el Package RSAGA

install.packages("RSAGA")

library ("RSAGA")

Paso 3: Quizás el más importante, porque ahora debemos establecer un ambiente en R. La finalidad es poder usar el objeto, en este caso "env_test" para dirigir las funciones RSAGA. En algunos casos podemos indicar en donde tenemos ubicado nuestro software SAGA.

Paso 4: Vamos a importar nuestro archivo raster que está en formato tif (dem_rimac.tif), para convertirlo en el formato de SAGA (dem_test.sgrd)

rsaga.import.gdal("F:/WORKS/SAMPLES/rtest/dem_rimac.tif", "dem_test.sgrd", env = env_test)

Paso 5: Verificamos la lista de la librería de módulos disponibles con las que podemos trabajar, es decir, todas las herramientas que a través de SAGA podemos ejecutar dentro del ambiente de R. Para nuestro caso vamos a incidir en aquellos relacionados al análisis del territorio. Es importante revisar los nombres, porque deberán ser considerados en nuestra línea de código para cuando sea requerido un módulo en específico.

rsaga.get.libraries(path=env_test$modules)

Figura 3: Lista de la librería de módulos

Paso 6: Vamos a lista la lista de funciones de geoprocesamiento que trabajaremos, teniendo en cuenta lo remarcado en la Figura 3. Por ahora nos dedicaremos a lo relacionado a "ta_morphometry".

rsaga.get.modules(lib = "ta_morphometry", env=env_test)

Figura 4: Funciones relacionados a morfometría

Paso 7: Configurar las herramientas de geoprocesamiento llamando a los módulos dentro de la biblioteca "ta_morphometry". Se pueden emplear cualquiera de los módulos que fueron listados previamente, solo se debe cambiar en "module = _ _ ". Tener en cuenta que también se puede colocar el número de código asignado, por ejemplo para "TPI Based Landform Classification", podemos ingresar "19".

rsaga.get.usage(lib= "ta_morphometry", module = "Slope, Aspect, Curvature", env = env_test)

rsaga.get.usage(lib= "ta_morphometry", module = 19, env = env_test)


Al ejecutar tendremos una vista de la lista de los parámetros (entrada y salida) a considerar, sobre todo, es importante verificar aquí las opciones que tenemos para personalizar nuestros resultados. También nos ayuda a revisar los formatos de ingreso (si es un texto o número), las unidades con los que se trabaja y cuales son las opciones por defecto, esto en caso que no ingresemos un parámetro. Para el primero se tendrá un resultado similar al siguiente.

> rsaga.get.usage(lib= "ta_morphometry", module = "Slope, Aspect, Curvature", env = env_test)
library path: C:\PROGRA~1\SAGA-GIS\modules\
library name: ta_morphometry
library     : Morphometry
Usage: saga_cmd ta_morphometry 0 -ELEVATION [-SLOPE ] [-ASPECT ] [-C_GENE ] [-C_PROF ] [-C_PLAN ] [-C_TANG ] [-C_LONG ] [-C_CROS ] [-C_MINI ] [-C_MAXI ] [-C_TOTA ] [-C_ROTO ] [-METHOD ] [-UNIT_SLOPE ] [-UNIT_ASPECT ]
  -ELEVATION:   Elevation
Grid (input)
  -SLOPE:       Slope
Grid (output)
  -ASPECT:     Aspect
Grid (output)
  -C_GENE:     General Curvature
Grid (optional output)
  -C_PROF:     Profile Curvature
Grid (optional output)
  -C_PLAN:     Plan Curvature
Grid (optional output)
  -METHOD:     Method
Choice
Available Choices:
[0] maximum slope (Travis et al. 1975)
[1] maximum triangle slope (Tarboton 1997)
[2] least squares fitted plane (Horn 1981, Costa-Cabral & Burgess 1996)
[3] 6 parameter 2nd order polynom (Evans 1979)
[4] 6 parameter 2nd order polynom (Heerdegen & Beran 1982)
[5] 6 parameter 2nd order polynom (Bauer, Rohdenburg, Bork 1985)
[6] 9 parameter 2nd order polynom (Zevenbergen & Thorne 1987)
[7] 10 parameter 3rd order polynom (Haralick 1983)
Default: 6
  -UNIT_SLOPE: Slope Units
Choice
Available Choices:
[0] radians
[1] degree
[2] percent
Default: 0
  -UNIT_ASPECT: Aspect Units
Choice
Available Choices:
[0] radians
[1] degree 


Paso 8: Ejecutar el Geoprocesamiento, para lo cual existen dos maneras fundamentalmente, el primero, cuando existe una función específica del RSAGA, tenemos la opción de obtener nuestro resultado. Para nuestro ejemplo, puede darse el caso que deseamos obtener la pendiente, el aspecto y la curvatura en un mismo procedimiento, para ello podemos escribir:

rsaga.slope.asp.curv("dem_test.sgrd", out.slope = "slope_dem1", out.cprof = "cprof_dem1", out.cplan = "cplan_dem1", method = "poly2zevenbergen", env=env_test)

También se pude realizar de manera individual, es decir cuando existe la función para calcula la pendiente:

rsaga.slope("dem_test.sgrd", "slope_dem2", method = "maxslope", env=env_test)

Por último existe la opción, el cual me parece más interesante, porque tenemos la posibilidad de obtener el mismo resultado si ejecutamos la función del Geoprocesador, para lo cual debemos guiarnos de lo mostrado como producto del Paso 7.

rsaga.geoprocessor(lib = "ta_morphometry", module = 0, param = list(ELEVATION = "dem_test.sgrd",SLOPE = "slope_dem3.sgrd", ASPECT = "aspect_dem2.sgrd", C_GENE= "curv_gener1.sgrd", METHOD= "0", UNIT_SLOPE="percent", UNIT_ASPECT="degree"), env = env_test)

Figura 5: Salida del Geoprocesamiento realizado

Paso 9: Ahora vamos a realizar un paso previo para lograr visualizar nuestros resultado. Como vieron en la Figura 5. los archivos se guardaron en formato de SAGA (*.sgrd), el cual lo podemos abrir no solo en ese Software, sino también en QGIS por ejemplo, pero para hacerlo más interesante lo vamos a convertir a un formato tipo ASCII (*.asc). Para ello, volvamos a revisar la Figura 3 e identifiquemos el módulo que requerimos, en este caso hablamos de "io_grid"

rsaga.get.modules(lib="io_grid", env=env_test)

Figura 6: Módulos para gestión de archivos raster

rsaga.get.usage(lib = "io_grid", module =0, env=env_test)

rsaga.geoprocessor(lib="io_grid",module=0,param=list(GRID="dem_fill1.sgrd",FILE="dem_esri.asc"),  env=env_test)

rsaga.geoprocessor(lib="io_grid",module=0,param=list(GRID="slope_dem4.sgrd",FILE="dem_slope.asc"), env=env_test)

rsaga.geoprocessor(lib="io_grid",module=0,param=list(GRID="aspect_dem4.sgrd",FILE="dem_aspect.asc"), env=env_test)


En caso se tengan muchos archivos, existe también una manera más rápida que nos permite realizar este proceso de conversión, a continuación muestro un ejemplo que logré ejecutar, considerando que estuve probando otros módulos.

rsaga.sgrd.to.esri(in.sgrds = c("cprof_dem4", "cplan_dem4", "carea1", "wetindex1", "hillshade1", "direction_dem3","connet_dem3", "order_dem3", "basin_dem3"), out.grids = c("esri_cprof", "esri_cplan", "esri_carea", "esri_wetindex", "esri_hillshade", "esri_direction", "esri_connect", "esri_order", "esri_basin"), out.path = getwd(), format = "ascii", env = env_test)

Paso 10: Vamos a definir los archivos raster para que podamos plotearlos

Es importante tener instalado y activado el Package "raster" y "sp"

plot(elevation)

plot(slope)

plot (aspect)

Visualización en Google Earth

Si bien no lo considero un paso adicional, me encontré con el siguiente enlace, a partir del cual traté de realizar el procedimiento, aunque no me siento muy satisfecho con los resultados, si alguien que lo intente a ver si le sale mejor, lo pueda compartir. En sí, primero debemos remontarnos a nuestra Figura 6, para definir nuestro módulo, en este caso será el "10", porque nuestro objetivo es convertirlo a un archivo PNG (que nos permitirá hacer transparencias).

Es importante tener instalado y activado el Package "rgdal", "maptools"

rsaga.get.modules(lib="io_grid", env=env_test)

rsaga.get.usage(lib= "io_grid", module = "10", env=env_test)


# Generando nuestro archivo tipo PNG

rsaga.geoprocessor(lib="io_grid", module=10, param=list(IMAGE="aspect_dem4.sgrd", FILE="aspect1.png"))

# Definiendo nuestro raster

# Creando nuestro objeto KML

# Creando nuestro archivo KML, con la función kmlOverlay, para mayor detalle revisar aquí.



kmlOverlay(aspect.kml, kmlfile="aspect1.kml", imagefile="aspect1.png", name="Aspecto_rimac")

Figura 7: Salida en Google Earth 

Muy bien, por mi parte seguiré investigando sobre las ventajas que se tiene el poder trabajar dentro de R, sobre todo a partir de nuestros productos obtenidos, aprovechando las ventajas de funciones estadísticas que podemos explotar. Espero que se animen a probarlo.

Teledetección con ORFEO Toolbox y QGIS Parte II

A pesar que ha pasado mucho tiempo, voy a cumplir con lo prometido, es decir mostrar cómo podemos calcular el NDVI (Normalized Difference Vegetation Index) empleando Monteverdi2, el cual viene hacer una de sus principales aplicaciones del ORFEO Toolbox, tal como vimos en la anterior entrada cuando lo describimos.

Dentro del procesamiento digital de imágenes realizaremos lo que vendría hacer una transformación de imagen, debido a que estamos involucrando el uso de múltiples bandas de datos, en general dichas transformaciones generan "nuevas" imágenes de dos o más fuentes que resaltan las características particulares o propiedades de interés, mejor que las imágenes de entrada. 

¿Qué es el NDVI?

Figura. 1: NDVI  - MODIS (Moderate Resolution Imaging Spectroradiometer)

Fuente: http://neo.sci.gsfc.nasa.gov/

Este contenido tiene como objetivo aprender a calcular un índice basado en cocientes, para lo cual vamos a conocer una técnica muy empleada para el monitoreo de las condiciones de vegetación, nos estamos refiriendo al Indice de Vegetación de Diferencia Normalizada (NDVI). Sin ánimos de profundizar el tema, podemos decir que el NDVI indica la abundancia y el estado de la vegetación, basados en el comportamiento reflectivo peculiar de la vegetación. En resumen, la firma espectral característica de la vegetación sana muestra un fuerte contraste entre la baja reflectividad en el rojo y una alta reflectividad en el infrarrojo cercano; esta diferencia es tanto mayor cuanto mayor es la densidad de la vegetación y mejor su estado fitosanitario. 

Los valores del NDVI esta acotado entre -1 a 1 y presenta la siguiente ecuación:

Donde:

IRp = reflectividad en el infrarrojo cercano

R    = reflectividad en el rojo

Calculando el NDVI con el Monteverdi 2

Para realizar este ejercicio emplearemos los mismos datos compartidos en la anterior entrada, es decir lo pueden descargar desde aquí.

Ahora debemos tener claro las bandas del satélite LandSat 8 con las que vamos a trabajar, para lo cual repasemos la siguiente tabla:

De acuerdo a la tabla solo necesitaremos emplear la banda 5 (NIR) y la banda 4 (rojo).

• Como primer paso abrimos el software Monteverdi 2 V0.8, si no lo tienen lo pueden descargar desde aquí.

• Ahora que tenemos los datos y el programa abierto, hacemos clic sobre File --> Import Image y buscamos los dos archivos con los que vamos a trabajar: LC80060662014227_B4 y LC80060662014227_B5.

Figura. 2: Vista del Monteverdi 2 con las bandas 4 y 5

• Nuestro siguiente paso será buscar el algoritmo que nos permita generar el NDVI a partir de ambas bandas, para ello dentro del OTB Applications Browser colocamos la palabra clave "math", porque lo que vamos a realizar será una operación matemática para reproducir la ecuación del NDVI.

Figura 3: Algoritmo del OTB requerido

• Luego con hacer clic sobre Band Math nos aparece en el escritorio del programa una pestaña específica para ingresar los datos requeridos, para ello simplemente arrastramos desde los Datasets, tal como se aprecia en la siguiente figura:

Figura 4: Ingresando las bandas 4 y 5 para ejecutar el algoritmo

• Ahora debemos especificar la ecuación que vamos a considerar y darle también un nombre al archivo raster de salida. Para entender cómo reconoce el Monteverdi las bandas, se debe tener en cuenta que como se trata de bandas independientes cada una de ellas tendría el siguiente esquema "imN1bN2", en donde N1 es el número de imagen considerado y N2 la banda que se considera, es decir que si trabajamos con una imagen que presenta una combinación de bandas, el valor de N2 puede variar desde el 1 hasta el número máximo de bandas que compone la imagen. Para nuestro caso solo se considera una sola banda por imagen, por lo tanto la ecuación sería tal como se aprecia en la figura de abajo.

Figura 5: Generación de la ecuación para obtener el NDVI

• Luego de hacer clic en "Execute", tendremos el siguiente resultado en nuestra vista.

Figura 6: Resultado del  cálculo del NDVI

• Para contar con una mejor visualización, se recomienda seguir los pasos de la entrada anterior cuando aprendimos a realizar composición de bandas, es decir que deberíamos realizar la composición de las bandas B2, B4 y la que acabamos de generar (ndvi_opc1), cuando obtenemos esa composición dirigimos los canales de las bandas adecuadas, es decir que el canal que correspondiente al verde hacer que coincida con el NDVI y con ello tendremos la siguiente imagen.

Figura 7: Ajuste de bandas para visualizar el NDVI

• Como consideración final existe una manera más directa para realizar el cálculo del NDVI, aunque personalmente considero que escribir la ecuación completa es lo más recomendable. En la parte donde se coloca la ecuación debemos incorporar para nuestro caso lo siguiente.

• Recomiendo también que practiquen con imágenes compuestas y generen ecuaciones que permitan no solo el cálculo del NDVI sino otros índices relacionados al monitorio de vegetación.

Calculando el NDVI dentro del QGIS

Para complementar lo ofrecido vamos a ver las opciones que tenemos si nos encontramos dentro del QGIS, en esta oportunidad contamos con la facilidad de poder tomar varios caminos, es decir podemos emplear el GRASS que se ubica dentro de nuestra caja de herramientas de procesado (el complemento Processing debe estar activado), a través del r.mapcalculator , para lo cual tendremos que introducir la formula del NDVI.

Figura 8: Cálculo del NDVI con GRASS

La otra opción que tenemos y que no solo nos proporciona el NDVI sino otros índices de vegetación, lo conseguimos a través del SAGA, presente también dentro de la caja de herramientas de procesado con el algoritmo Vegetation Index (slope based), el cual presenta la característica de calcular hasta siete tipos de índices de vegetación. Para mayor detalle de éstos índices se recomienda revisar a partir de la página 249 de la publicación denominada User Guide for SAGA (version 2.0.5).

Cuando empleamos esta opción, nuestra ventana debe ser llenado tal como se muestra en la Figura 9, teniendo en cuenta que para nuestro caso solo nos interesa que se muestre el NDVI.

Figura 9: Cálculo de Indices de vegetación con SAGA 

Una vez realizado tanto a través del GRASS como por el SAGA nos van a brindar los mismos resultados que los obtenidos con el Monteverdi. Ahora de la misma manera podemos ajustar las bandas para tener una mejor visualización del NDVI, para ello debemos generar una composición de bandas considerando la que acabamos de calcular; por lo generar dentro del QGIS esto se logra fácilmente si nos vamos a herramientas de Ráster ---> Miscelánea --->Combinar...

Figura 10: Herramientas ráster para realizar composición de imágenes

Luego de seleccionar las tres (3) bandas o imágenes, simplemente hay que asegurarse de hacer solamente check en donde indica "pila de capas", los demás deben estar libres. Finalmente dentro de "Estilo" ajustamos las bandas a los respectivos canales, cuidando de que la imagen del NDVI se ubique en el canal de la banda verde.

Figura 11: Ajuste de las bandas a los canales de colores respectivos

Como resultado final tendremos una imagen como la obtenida anteriormente y si activamos el complemento de "Value Tool", podemos desplazarnos con el mouse y ver los valores de píxel de las tres bandas.

Figura 12: Resultado final con la imagen del NDVI calculado

Espero que esta entrada haya despertado el interés del uso de esta herramienta, en una próxima oportunidad se realizará otro tipo de procesamiento digital de imágenes.