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.