1. Procesos preliminares

Inicialmente, es de tener en cuenta que para poder desarrollar los procedimientos que en la presente actividad se describen, debe tener instalados en su ordenador los siguientes softwares (programas):

  • R y su entorno de desarrollo integrado RStudio.
  • GRASS-GIS, es de considerar que, las únicas versiones compatibles para hacer el acople con R son de la 7.4 a la 7.6 (preferiblemente esta última).
  • SAGA-GIS

Asimismo, es de tener presente que, los datos de entrada (input) serán un modelo de elevación digital (DEM, por sus siglas en inglés) y las coordenadas del punto de desfogue (punto de salida) de la cuenca que, por lo general, se encuentra representado por una estación de medición/muestreo. Para esta actividad, el DEM se descargó/obtuvo del repositorio HydroSHEDS (https://www.hydrosheds.org/) y se determinó como zona de estudio la cuenca hidrográfica del río Coello (SZH 2121), del departamento de Tolima (Colombia).

De acuerdo con lo descrito, se parte antes de comenzar la delimitación con la limpieza del ambiente de trabajo, con el fin de que variables cargadas en procesos anteriores no llegasen a interferir, por lo cual, se usó la siguiente línea de comando:

rm(list=ls(all=TRUE))

Y, por otro lado, se debe direccionar a la carpeta de trabajo, esto se realiza por medio del comando setwd(“direccion_carpeta_de_trabajo”) y a partir de ser aplicado, se verifica que se haya dirigido a través del comando getwd(). No se muestra, ya que la dirección es de un ordenador personal. Además, en ella deben estar todos los archivos a trabajar.

1.1. Cargar e instalar librerías

A su vez, para poder desarrollar los procesos se deben instalar y habilitar las siguientes librerías:

#install.packages("raster")
#install.packages("RSAGA")
#install.packages("sf")
#install.packages("latticeExtra")
#install.packages("rgdal")
#install.packages("ggplot2")
#install.packages("ggspatial")
#install.packages("rgrass7")

library("rgrass7") #Para que R articule con GRASS-GIS
## Error in xmlTreeParse(tr) : empty or no content specified
library("raster") #Para que se puedan leer archivos ráster
library("RSAGA") #Para que R articule con SAGA-GIS
library("sf")
library("latticeExtra")
library("rgdal")
library("ggplot2") #Para las gráficas
library("ggspatial") #Para los datos espaciales en las gráficas

2. Delimitación de una cuenca hidrográfica por GRASS-GIS

2.1. Preliminares en R para habilitar GRASS-GIS y generación de espacio de trabajo

El primer paso consiste en iniciar GRASS-GIS desde R, para ello con la siguiente línea se direcciona a la carpeta donde se encuentra el programa y, también, donde el programa almacena por defecto, que como es conocido es grassdata.

NOTA: tenga en cuenta cada uno de los flag, son importantes, para adquirir información de ellos y, también, a la hora de adquirir datos del comando (su uso y demás), se puede consultar dentro repositorio de GRASS-GIS.

En adición, para el caso, desde aquí se creó una nueva location que, de cierta forma, obliga a que el primer mapset creado sea PERMANENT (sí se corriese desde el programa este al generar la nueva location, crea el mapset por defecto), únicamente con escribir el nombre se logra, como se muestra a continuación. Es de aclarar que, sí la location y el mapset ya existen, simplemente se llaman de la misma manera y, se puede llamar sólo la location (guardará por defecto en PERMANENT), pero no se puede llamar sólo al mapset.

initGRASS(gisBase = "C:/Program Files/GRASS GIS 7.6",
          gisDbase = "C:/Users/usuario/OneDrive - Pontificia Universidad Javeriana/Documents/grassdata",
          location ="WGS84_4326_MPH", mapset="PERMANENT", override=T)
## Warning in system(syscmd, intern = intern, ignore.stderr = ignore.stderr, :
## running command 'g.proj.exe -w' had status 1
## Warning in system(syscmd, intern = intern, ignore.stderr = ignore.stderr, :
## running command 'g.proj.exe -j' had status 1
## gisdbase    C:/Users/usuario/OneDrive - Pontificia Universidad Javeriana/Documents/grassdata 
## location    WGS84_4326_MPH 
## mapset      PERMANENT 
## rows        1 
## columns     1 
## north       1 
## south       0 
## west        0 
## east        1 
## nsres       1 
## ewres       1

Como se pudo apreciar con el anterior paso, se creó la location sin una proyección de coordenadas, por lo tanto, se debe definir la proyección en la que se encuentra este espacio, para ello se usa la siguiente línea, teniendo presente que para este caso se va a utilizar el sistema de coordenadas geográficas (WGS84, EPSG:4326). Sí la carpeta existiese al correr la línea de código anterior hubieran aparecido datos sobre la location y el mapset, en los que se expone la proyección y coordenadas entre las que se encuentra la región de trabajo (por defecto).

execGRASS("g.proj", proj4="+proj=longlat +datum=WGS84 +no_defs", flag="c")

Luego, aunque es opcional, se recomienda crear un nuevo mapset en el que se guarde todo el proceso que se va a realizar:

execGRASS("g.mapset", mapset="Cuenca_Coello", flag=c("c","overwrite"))

Después, se debe verificar que se el sistema se encuentre en la carpeta creada para trabajo (una vez se crea el mapset este lo direcciona por defecto, sin embargo, es mejor confirmar).

execGRASS("g.gisenv", flag="n")

#Sí se tuviera que cambiar de mapset
#execGRASS("g.mapset", mapset="Cuenca_Coello", location ="WGS84_4326_MPH")

De esta manera, se pasa a identificar la región de trabajo en la que se encuentra el ambiente de trabajo.

execGRASS("g.region", flag=c("p")) #Información de la región de trabajo

2.2. Extraer información de archivos ráster y vector desde R

Facilita un poco las cosas interpretar algunos de los metadatos de los archivos tanto ráster como vector desde R.

2.2.1. Para el caso del ráster

Se debe cargar el DEM a las variables del ambiente en R:

dem = raster("hs_cua_dem.tif")
#Se carga el DEM, importante que esté dentro de la carpeta de trabajo

A partir de ello, se obtiene la información del DEM (esta se tendrá en cuenta para establecer la región de trabajo):

#Propiedades del DEM
projection(dem) #Se pregunta en qué proyección viene el DEM
## [1] "+proj=longlat +datum=WGS84 +no_defs"
crs(dem) #Es mucho más detallada
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["WGS 84 (with axis order normalized for visualization)",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]
dem #Obtener información del DEM, aquí veo las coordenadas para establecer región de trabajo
## class      : RasterLayer 
## dimensions : 700, 1072, 750400  (nrow, ncol, ncell)
## resolution : 0.0008333333, 0.0008333333  (x, y)
## extent     : -75.69588, -74.80255, 4.149485, 4.732818  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : hs_cua_dem.tif 
## names      : hs_cua_dem 
## values     : 211, 5212  (min, max)
#dem_repro=projectRaster(dem,crs=("+init=epsg:4326")) #Sí se tuviese que reproyectar
#writeRaster(dem_repro, filename="dem_repro.tif", format="GTiff", overwrite=TRUE) #Sí ese DEM reproyectado se tuviera que sacar para luego ser cargado en GRASSGIS
  • Se genera una gráfica inicial, simple, sólo para confirmar que el DEM quedó cargado:
#Observar cómo quedó cargado (tener en cuenta la proyección)
spplot(dem, col.regions=terrain.colors(255)) #Gráfica con escala de colores

NOTA: el DEM fue procesado previamente, siendo recortado al área de trabajo, no obstante, para realizar el proceso desde aquí, se debe ya haber definido la región de trabajo para que por medio del comando r.mask se haga una máscara que permita realizar el proceso más rápido y sin complicaciones, ya que se reduce el área en la que el programa entraría a correr.

2.2.2. Para el caso del archivo vector

Por otra parte, el archivo vector de geometría punto, que es el punto de salida (para el caso, estación MARANONES [21215090]) se carga al ambiente de la siguiente manera:

Punto_salida = readOGR(".","EM_sal_WGS84") #Cargar archivo .shp en R
class(Punto_salida) #Ver el tipo de .shp y/o datos

Más información del archivo shapefile (de aquí salen las coordenadas que se usarán más adelante para la delimitación de la cuenca):

proj4string(Punto_salida) #Ver la proyección del archivo
## [1] "+proj=longlat +datum=WGS84 +no_defs"
#Punto_salida_wgs84 = spTransform(Punto_salida,CRS(wgs84)) #Sí se tuviera que reproyectar

coordinates(Punto_salida) #Coordenadas del .shp
##      coords.x1 coords.x2
## [1,] -74.93295  4.218239

2.3. Definir región de trabajo

Ya con todos lo datos se entra a definir la región de trabajo para la cual GRASS-GIS recorre para hacer los procesos que se le indiquen. Sin embargo, es de tener en cuenta que esto en primera instancia debería ser hecho a PERMANENT.

#Se retorna al mapset PERMANENT para poder definir la región de trabajo
execGRASS("g.mapset", mapset="PERMANENT", location ="WGS84_4326_MPH")
#Se verifica que se haya cambiado de mapset
execGRASS("g.gisenv", flag="n")
#Y, teniendo como insumo las coordenadas y la resolución del DEM se define la región de trabajo
execGRASS("g.region", n="4.732818", s="4.149485", e="-74.80255", w="-75.69588", res="00:00:03", flag=c("s","p")) #Se pasa la resolución del DEM a grados:minutos:segundos (es opcional)

No obstante, ya que no quedó definido para el mapset que se ha venido desarrollando, se debe realizar el proceso, forzando un poco a que asuma la región de trabajo ya identificada.

#Nuevamente, se regresa al mapset de trabajo
execGRASS("g.mapset", mapset="Cuenca_Coello", location ="WGS84_4326_MPH")
#Se vuelve a verificar que se haya hecho el cambio de mapset
execGRASS("g.gisenv", flag="n")
#Se comprueba que la region de trabajo del mapset haya quedado en la definida
execGRASS("g.region", flag=c("p"))

#Debería definirse la región de trabajo una vez se define dentro de PERMANENT, pero como no se vuelve hacer el mismo proceso
execGRASS("g.region", n="4.732818", s="4.149485", e="-74.80255", w="-75.69588", res="00:00:03", flag=c("d","p"))
execGRASS("g.gisenv", flag="n")
execGRASS("g.region", flag=c("p"))

#Y se verifica la proyección del mapset nuevamente
execGRASS("g.proj", flag=c("p"))

2.4. Cargar DEM dentro del ambiente de GRASS-GIS

Ya habiendo desarrollado todos los procesos anteriores con satisfacción, se es posible cargar el DEM dentro del ambiente de GRASS-GIS para empezar con el proceso de delimitación de la cuenca

#Subir DEM a GRASS-GIS
execGRASS("r.in.gdal", input="C:/Users/usuario/OneDrive - Pontificia Universidad Javeriana/Desktop/RMarkdown_GRASSGIS_y_SAGAGIS/hs_cua_dem.tif",
          output="DEM", flag=c("overwrite"))

#Tener en cuenta sí se deben unir los DEM, para eso usar r.patch (cómo usar el comando en el repositorio)

#Cargar un archivo .shp a GRASSGIS
#execGRASS("v.in.ogr", input="direccion_de_trabajo/archivo.shp",
#          output="nombre_salida")

2.4.1. Corrección del DEM

Se realiza una pequeña corrección al DEM, estableciendo datos nulos o menores a 0 como 0.

execGRASS("r.null", map="DEM", setnull="0")

2.5. Delimitación de la cuenca con GRASS-GIS

Antes de llegar al paso de delimitación, se debe generar un conjunto de ráster que ayudan a dicho proceso y/o que contribuyen en su desarrollo, y auqe juegan un papel importante, estos son: direcciones y acumulaciones de drenajes, aunque también se calculan las corrientes y el mapa de pendientes de inclinación. Asimismo, se define un threshold de 50, dado el tamaño y reolución del ráster, ya que como se entiende será la cantidad de celdas que a partir de las que comenzará a acumular alrededor de una.

execGRASS("r.watershed", elevation="DEM", threshold=50, accumulation="accumulation",
          drainage="drainage_dir", stream="stream", slope_steepness="slope_steepness",
          flag=c("s", "overwrite"))

IMPORTANTE: sí se cuenta con depresiones en el terreno se puede y debe cargar para que sean tenidas en cuenta a la hora de delimitar la cuenca.

También, es de tener presente que, posiblemente la estación no está sobre la celda de acumulación, así que, probablemente el sistema la redistribuye/reubica el punto, pero dentro de lo posible o mejor es prudente conocer sí realmente se encuentra dentro de la corriente.

Es así como se llega ya a la delimitación de la cuenca hidrográfica, donde inicialmente tendrá como entrada el mapa de direcciones de drenaje generado en el anterior proceso y se obtendrá como salida un archivo ráster de la delimitación de la cuenca del río Coello.

## Para el caso, la cuenca del río Coello
execGRASS("r.water.outlet", input="drainage_dir", output="Cuenca_coello", coordinates=c(-74.93295, 4.218239),flag=c("overwrite"))

Ya que, como se mencionó, el archivo de salida es tipo ráster se convierte en un vector de geometría polígono.

# De ráster a vector
execGRASS("r.to.vect", input="Cuenca_coello", output="Cuenca_coello", type="area", flag=c("s", "overwrite"))

Y, para poder realizar su representación por medio de un mapa desde R, el archivo vectorial generado se extrae del GRASS-GIS en una carpeta, donde se le asigna un nombre al archivo.

# Exportar archivos

## Sólo el .shp
execGRASS("v.out.ogr", input="Cuenca_coello", type="area", format="ESRI_Shapefile",
          output="C:/Users/usuario/OneDrive - Pontificia Universidad Javeriana/Desktop/RMarkdown_GRASSGIS_y_SAGAGIS/Cuenca_coello_GRASS.shp",
          flag=c("c", "overwrite"))

Por último, se crea el mapa de la cuenca generada:

#Se lee el archivo .shp y se define como variable
Cuenca_coello_GRASS = st_read("Cuenca_coello_GRASS.shp")

#Mapa
ggplot(data = Cuenca_coello_GRASS) +
  geom_sf(color="black") +
  xlab("Longitud") + ylab("Latitud") +
  ggtitle("Cuenca hidrográfica del río Coello", 
          subtitle = "Departamento de Tolima") +
  annotation_scale() +
  annotation_north_arrow(location='tr')

#Sacar/extraer el mapa de R a archivo .tif
##Relación ancho y alto de la imagen
ancho=1750
alto=ancho/sqrt(2)
## Aquí inicia el proceso de la imagen
tiff(filename="Mapa_cuenca_GRASSGIS.tif",width=ancho,height=alto,
     compression="lzw",pointsize=10,bg="white",res=300)
par(mar=c(4,4,.5,.5))
ggplot(data = Cuenca_coello_GRASS) +
  geom_sf(color="black") +
  xlab("Longitud") + ylab("Latitud") +
  ggtitle("Cuenca hidrográfica del río Coello", 
          subtitle = "Departamento de Tolima") +
  annotation_scale() +
  annotation_north_arrow(location='tr')
dev.off()

3. Delimitación de una cuenca hidrográfica por SAGA-GIS

3.1. Preliminares en R para habilitar SAGA-GIS

Así como en GRASS-GIS, para que R corra SAGA-GIS se debe direccionar el comando a que lea la ubicación del programa.A su vez, es de tener presente que se tiene habilitada la librería.

envi = rsaga.env(path="C:/Users/usuario/Downloads/saga-8.1.1_x64/saga-8.1.1_x64")
## Verify specified path to SAGA command line program...
## Found SAGA command line program. Search for not specified SAGA modules path... 
## Done
rsaga.get.version(env = envi)
## [1] "8.1.1"

3.2. Proyectar un DEM dentro de R y generar archivo dem.sgrd

Dado que para este caso, se pretende hacer la delimitación de la cuenca con una proyección distinta a la WGS84, se carga el DEM y se reproyecta a través de comandos de R, de la siguiente manera:

#Se vuelve a cargar el DEM (es para tenerlo en cuenta, ya que igual para el caso aún se encuentra en el ambiente)
dem = raster("hs_cua_dem.tif")
#Se pregunta en qué proyección se encuentra el DEM
projection(dem)
## [1] "+proj=longlat +datum=WGS84 +no_defs"
#Se reproyecta a la que se desea (para el caso es al sistema cartesiano origen central de Colombia, EPSG: 3116)
dem_3116 = projectRaster(dem,crs=("+init=epsg:3116"))
#Sólo por corroborar se vuelve a preguntar la proyección
projection(dem_3116)
## [1] "+proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs"

Se hace una gráfica simple sólo para corroborar el proceso.

spplot(dem_3116, col.regions=terrain.colors(255))

De esta forma, se continua a cargar el DEM a SAGA-GIS, pero para ello, se debe crear un archivo .sgrd, así que, primero se extrae el DEM reproyectado desde R en un formato .tif a la carpeta en la que se está trabajando y, luego, se procede a generar el archivo. Vale aclarar que, cuendo se crea el archivo .sgrd, también se produce uno de extensión .sdat que es como tal el archivo ráster.

#Exportar DEM en archio .tif
writeRaster(dem_3116, filename="dem_3116.tif", format="GTiff", overwrite=TRUE)

#Crear archivo .sgrd
rsaga.geoprocessor(lib = "io_gdal", 0,
                   param = list(GRIDS ="dem_3116_sgrd",
                                FILES= "dem_3116.tif",
                                TRANSFORM= 1),env = envi)

#Verificar que quedó creado
dem_3116_sgrd = raster("dem_3116_sgrd.sdat")
#Verificar la proyección
projection(dem_3116_sgrd)

3.2.1. Corrección del DEM

Es así como se pasa a una rápida corrección del DEM. No tiene en cuenta el Threshold.

rsaga.geoprocessor(lib = "ta_preprocessor", 2,
                   param = list(DEM ="dem_3116_sgrd.sgrd",
                                DEM_PREPROC= "demsink_fill"),env = envi)

3.3. Cargar y reproyectar el archivo .shp

Así como se hizo con el archivo ráster, se debe cargar el archivo vector de geometría punto que hace las veces de la estación de referencia para el punto de salida de la cuenca y, también debe pasar a ser reproyectado, ya que con esta se sabe las coordenadas que se introducen a la delimitación, Recuerde, dentro de lo posible esta debe estar lo más adentro de la corriente para ser valorada de manera adecuada por el mapa de acumulación.

Punto_salida = readOGR(".","EM_sal_WGS84") #Cargar archivo .shp en R
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\usuario\OneDrive - Pontificia Universidad Javeriana\Desktop\RMarkdown_GRASSGIS_y_SAGAGIS", layer: "EM_sal_WGS84"
## with 1 features
## It has 26 fields
class(Punto_salida) #Ver el tipo de .shp
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
proj4string(Punto_salida) #Proyección en la que se encuentra el archivo vector
## [1] "+proj=longlat +datum=WGS84 +no_defs"
#Reproyección a EPSG: 3116
Punto_salida_3116 = spTransform(Punto_salida,crs("+proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs"))
proj4string(Punto_salida_3116) #Verificar que se haya hecho la reproyección
## [1] "+proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs"
coordinates(Punto_salida_3116) #Preguntar coordenadas del .shp
##      coords.x1 coords.x2
## [1,]  905025.1  958256.9

3.4. Delimitación de la cuenca con SAGA-GIS

Ya con los insumos adquiridos gracias a la implementación de los procedimientos anteriores, se pasa a modelar la divisoria de aguas de la cuenca del río Coello (delimitación de la cuenca hidrográfica). Es de notar que, acá se obvió el threshold.

#rsaga.get.modules("ta_hydrology",env = envi) #muestra las opciones/métodos/parámetros del comando
#rsaga.get.usage('ta_hydrology', 4,env = envi) #muestra las opciones/métodos/parámetros del comando, pero más específico, ya para la opción 4

#Delimitación de la cuenca -> salida = formato ráster
rsaga.geoprocessor(lib = "ta_hydrology", 4,
                   param = list(TARGET_PT_X =905025.1,
                                TARGET_PT_Y =958256.9,
                                ELEVATION = "demsink_fill.sgrd",
                                AREA = "Cuenca_coello.sgrd",
                                METHOD = 0),env = envi) #Observar que se puso como método el 0

Es así como, se convierte a un archivo vector de extensión .shp la cuenca hidrográfica del río Coello generada por herramientas de SAGA-GIS desde R.

rsaga.geoprocessor(lib = "shapes_grid", 6,
                   param = list(GRID = "Cuenca_coello.sgrd",
                                POLYGONS = "Cuenca_coello_SAGA.shp",
                                CLASS_ALL = 1,
                                CLASS_ID = 100,
                                SPLIT = 0),env = envi)

Finalmente, con los procesos ejecutados, se generan los mapas que, a diferencia de los elaborados para el proceso de GRASS-GIS, estos son un poco más desarrollados al haber incluido más parámetros.

#Se lee el archivo .shp y se define como variable
Cuenca_coello_SAGA = st_read("Cuenca_coello_SAGA.shp")

#Mapa
ggplot(data = Cuenca_coello_SAGA) +
  geom_sf(color="black") +
  coord_sf(datum = sf::st_crs(3116))+
  theme (axis.text.x = element_text(size=rel(0.7),vjust=0.5),
         axis.text.y = element_text(size=rel(0.7), angle=90, hjust=0.5))+
  xlab("Este") + ylab("Norte") +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, size=rel(1))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, size=rel(1))) +
  ggtitle("Cuenca hidrográfica del río Coello", 
          subtitle = "Departamento de Tolima") +
  annotation_scale() +
  annotation_north_arrow(location='tr')

# Sacar/extraer el mapa de R a archivo .tif
##Relación ancho y alto de la imagen
ancho=1750
alto=ancho/sqrt(2)
## Aquí inicia el proceso de la imagen
tiff(filename="Mapa_cuenca_SAGAGIS.tif",width=ancho,height=alto,
     compression="lzw",pointsize=10,bg="white",res=300)
par(mar=c(4,4,.5,.5))
ggplot(data = Cuenca_coello_SAGA) +
  geom_sf(color="black") +
  coord_sf(datum = sf::st_crs(3116))+
  theme (axis.text.x = element_text(size=rel(0.7),vjust=0.5),
         axis.text.y = element_text(size=rel(0.7), angle=90, hjust=0.5))+
  xlab("Este") + ylab("Norte") +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, size=rel(1))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, size=rel(1))) +
  ggtitle("Cuenca hidrográfica del río Coello", 
          subtitle = "Departamento de Tolima") +
  annotation_scale() +
  annotation_north_arrow(location='tr')
dev.off()

  1. Ingeniero ambiental, Universidad Santo Tomás de Colombia (Bogotá D. C.). Facultad de Ingenierías, Maestría en Hidrosistemas, Pontificia Universidad Javeriana (Bogotá D. C.).↩︎