Este cuaderno es una guía para facilitar el desarrollo de la actividades en clase de los cursos: Hidroclimatológia (pregrado y posgrado) y Modelación de procesos hidrológicos de la FEAR-PUJ, y fue consolidado con la ayuda de Jhon Pulgarin (estudiante de Maestría en Hidrosistemas en la PUJ).

El presente cuaderno inicia con el preprocesamiento del Modelos de Elevación Digital (DEM) y termina con la delimitación de la cuenca.

Antes de iniciar, es recomendable iniciar limpiando el espacio de trabajo y llamando de una vez todos las librerías que se vayan a necesitar. Si no se tiene instalada alguna de las librerías es necesario hacerlo con el comando “install.packages”. Así mismo, se selecciona el directorio de trabajo donde se tiene el archivo con la información a procesar.

Software instalados

Para poder desarrollar las actividades que aquí se proponen es necesario que en su equipo de cómputo tenga descargados e instalados de manera correcta los siguientes programas:

Datos para desarrollar ejercicio propuesto

Para este ejercicio práctico únicamente requieren como entradas un DEM (recortado) y las coordenadas del punto de cierre (para el caso es una estación de monitoreo limnimétrica ó limnigráfica).

1. Pre-procesamientos

Para iniciar, es importante que ninguna variable de procesos anteriores o alojada en el ambiente interfiera en los procesos a desarrollar, así que, para empezar se realiza una limpieza del ambiente de trabajo con el uso de la siguiente línea de comando:

rm(list=ls(all=TRUE))

Luego, se debe definir la carpeta de trabajo, para el caso, esta será la que almacenará tanto las entradas (DEM y punto de cierre) para llevar a cabo el proceso, como la que recibirá las salidas. Es por ello que, se aconseja que cree una carpeta nueva para esta actividad.

Conforme con lo mencionado, en las siguientes líneas se expresa cómo se direccionar la carpeta al entorno de trabajo:

setwd("C://Users//maria.arenasb//OneDrive - Pontificia Universidad Javeriana//Javeriana//Scripts - R")
getwd() #este comando le permite verificar que esta en la carpeta que indicó en el comando anterior
## [1] "C:/Users/maria.arenasb/OneDrive - Pontificia Universidad Javeriana/Javeriana/Scripts - R"

1.1. Instalar y habilitar librerías/paquetes de soporte

Recuerde, antes de iniciar cualquier proceso en un medio de programación (para el caso es R), siempre tendrá que cargar las respectivas librerías, paquetes o subseries que contienen los códigos de soporte de cada uno de los comandos.

Para ello, inicialmente deberá instalarlas (en caso de no tenerlas aún, se recomienda instalarlas una única vez, a menos de que R lo solicite, por ejemplo, para ser actualizada).

Para iniciar el proceso de isntalación de paquetes se recomienda quitar el comentario (eliminar el signo #) y correr cada una de las líneas. Al terminar, se aconseja volver a colocar en comentario.

# install.packages("rgrass")
# install.packages("RSAGA")
# install.packages("raster")
# install.packages("rgdal")
# install.packages("sf")
# install.packages("sp")
# install.packages("prettymapr")
# install.packages("ggplot2")
# install.packages("ggspatial")
# install.packages("shapefiles")

Una vez se logra instalar las librerías con éxito, se procede a habilitarlas de la siguiente manera:

library("rgrass")
library("RSAGA")
library("raster")
library("rgdal")
library("sf")
library("sp")
library("prettymapr")
library("ggplot2") # Para las gráficas
library("ggspatial") #Para los datos espaciales en las gráficas
library("shapefiles")

1.2. Cargar datos

Recuerde, en lo posible, que debe organizar todos los archivos dentro de la misma carpeta de trabajo. Asimismo, tenga en cuenta que el DEM ha tenido un pre-procesamiento, en este lo que se realizó fue un recorte del ráster al tamaño del área de estudio o zona de influencia para optimizar lso tiempos de ejecución.

Estacion_cierre=shapefile("Est_Paicol.shp")
## Warning: [vect] Z coordinates ignored
DEM=raster("mascara4.tif")

1.2.1. Reproyección de los archivos espaciales

Inicialmente, explore el sistema de coordenadas en el que se encuentran los archivos espaciales. Por medio de las siguientes líneas podrá ver la estructura PROJ4 del sistema de coordenadas en el que se encuentra.

Estacion_cierre@proj4string@projargs
# projection(Estacion_cierre) #Es similar al anterior
DEM@crs@projargs
# projection(DEM) #Es similar al anterior

A continuación, se llevan los archivos a una misma proyección cartográfica para realizar los procesos sin algún tipo de interferencia/error por esto.

Para el caso, se trabajá con el sistema de coordenadas plano MAGNA-SIRGAS / Origen-Nacional (EPSG: 9377), oficial para Colombia. Para mayor información dirijase a https://origen.igac.gov.co/docs/ABC_Nueva_Proyeccion_Cartografica_Colombia.pdf.

Tenga en cuenta que para este ejercicio, la reproyección no es tan representativa, ya que los archivos ya se encuentran proyectados a dicho sistema de coordenadas.

Estacion_cierre=spTransform(Estacion_cierre,crs("+init=epsg:9377"))
# projection(Estacion_cierre) = crs("+init=epsg:9377") # Otra forma de reproyectar shape
DEM=projectRaster(DEM,crs=crs("+init=epsg:9377"))
# projection(DEM) = crs("+init=epsg:9377") # Otra forma de reproyectar raster

Sí desea verificar lo hecho hasta ahora construya una gráfica rápida de la siguiente manera:

plot(DEM)+
plot(Estacion_cierre, add = TRUE)

## integer(0)

Finalmente, los archivos reproyectados se exportan a la carpeta de trabajo, esto para ser usados como insumo dentro de los procesos.

# Exportar vector
writeOGR(Estacion_cierre, dsn = "C://Users//maria.arenasb//OneDrive - Pontificia Universidad Javeriana//Javeriana//Scripts - R", 
         layer = "Estacion_cierre_paicol1",
         driver = "ESRI Shapefile",overwrite_layer=T)
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
# Exportar Ráster
writeRaster(DEM,"DEM_EPSG9377.tif",overwrite = T) #puede cambiar el nombre que se encuentra en verde

2. Uso de GRASS-GIS desde R

Para realizar el modelado de cuencas hidrográficas se usa el paquete rgrass (consultar en el CRAN: https://cran.r-project.org/web/packages/rgrass/rgrass.pdf).

Este se encuentra actualizado y aparentemente no genera conflicto con las versiones recientes o anteriores, como sí lo hacía rgrass7.

Grass se debe descargar desde este link https://grass.osgeo.org/news/2022_01_30_grass_gis_8_0_0_released/. Dirijase a la parte inferior de la página web y seleccione la versión Standalone para windows. Una vez descargado el archivo, haga doble clic e instalelo en la carpeta archivos de programa del disco C://.

2.1. Procesos iniciales en GRASS-GIS desde R

  • Comience nombrando las siguientes variables:
Location ="Loc_EPSG9377"
Mapset="Delimitacion_cuenca"
  • Continúe definiendo el Threshold (umbral), para establecerlo se puede determinar una relación, por la cual, se espera que a un (1) \(Km^2\) se empiecen a acumular píxeles para establecer la red del drenaje, entonces, esto se traduciría en 1 \(Km^2\) sobre la resolución del DEM también en \(Km^2\), que para el caso este último valor es de 90 m, ya que se está usando un DEM del proyecto HydroSHED de resolución 3 arc-second. Como se puede evidenciar a continuación en la ecuación (1).

\[\begin{equation} Threshold = \frac{1 Km^2}{(0.09^2) Km^2} = 123.46 \end{equation}\]

Aproximando por encima, el Threshold a usar será de 124 unidades.

Threshold=124 #este lo usamos por defecto para el raster de 90 m
  • Ahora, se debe iniciar GRASS-GIS desde R, para eso tiene que seguir los pasos descritos a continuación:
initGRASS(# Direccionar al ejecutable de GRASS-GIS, para el caso se trabajará con la
          # versión 8.0
          gisBase = "C://Program Files//GRASS GIS 8.0", # Puede ponerse como NULL
          
          # Direccionar la carpeta GRASSDATA, allí se guarda todo el trabajo realizado. Esta carpeta debe crearla antes de iniciar a correr el código.
          gisDbase = "C://Users//maria.arenasb//OneDrive - Pontificia Universidad Javeriana//Javeriana//Hidroclimatologia - Ecologia//grassdata",
          
          # Acá se crea una nueva Location, por ese se definió el nombre anterior
          location =Location,
          
          # El mapset creado por defecto es PERMANENT, pero por sí acaso se establece
          mapset="PERMANENT", # Se recomienda no modificar
          
          # Para que sobreescriba y no genere conflicto
          override=T)
## Warning in system(syscmd, intern = intern, ignore.stderr = ignore.stderr, :
## comando ejecutado 'g.proj.exe -w' tiene estatus 884
## gisdbase    C://Users//maria.arenasb//OneDrive - Pontificia Universidad Javeriana//Javeriana//Hidroclimatologia - Ecologia//grassdata 
## location    Loc_EPSG9377 
## mapset      PERMANENT 
## rows        1 
## columns     1 
## north       1 
## south       0 
## west        0 
## east        1 
## nsres       1 
## ewres       1 
## projection  +proj=tmerc +lat_0=4 +lon_0=-73 +k=0.9992 +x_0=5000000 +y_0=2000000
## +ellps=GRS80 +units=m +no_defs +type=crs
  • Definir proyección. Se continúa trabajando con MAGNA-SIRGAS / Origen-Nacional (EPSG: 9377), para ello, se llevaron a cabo los procesos anteriores.
execGRASS("g.proj", epsg=9377, flag=c("c"))

# Imprimir las propiedades de la location para comprobar que se proyectó al sistema de
# coordenadas
execGRASS("g.proj", flag=c("p"))
  • Continuando, aunque es opcional, es preferible que genere su propio mapset de trabajo, por eso definió el nombre atrás.
execGRASS("g.mapset", mapset=Mapset, flag=c("c","overwrite"))

# Verifique dónde se encuentra parado
execGRASS("g.gisenv", flag="n")

# Sí se tuviera que cambiar de mapset
# execGRASS("g.mapset", mapset=Mapset, location =Location)
  • Antes de finalizar, defina la región de trabajo, para poder hacer procesos en GRASS-GIS es muy importante y, para el caso, siendo un poco más fácil, esta se definirá a partir del DEM.
# Información de la region de trabajo hasta ahora (note los valores)
execGRASS("g.region", flag=c("p"))

# Cargar DEM al Mapset
execGRASS("r.in.gdal", input="DEM_EPSG9377.tif",output="DEM",
          flag=c("overwrite"))

# Definir region de trabajo a partir del DEM
execGRASS("g.region", raster="DEM", flag=c("p"))
  • Por último, con el fin de comprobar que los procesos hasta ahora se hayan hecho de la manera deseada corra las siguientes líneas.
execGRASS("g.gisenv", flag="n")
execGRASS("g.region", flag=c("p"))
execGRASS("g.proj", flag=c("p"))

2.2. Delimitación de cuencas hidrográficas con GRASS-GIS desde R

Logrado lo anterior, se procede a los siguientes pasos:

execGRASS("r.fillnulls", input="DEM",output="DEM_fill", method="bilinear",
          flag=c("overwrite"))
  • Generar las direcciones y acumulaciones del flujo.
execGRASS("r.watershed", elevation="DEM_fill", threshold=Threshold, accumulation="accum",
          drainage="drainage_dir", stream="stream", slope_steepness="slope_steepness",
          flag=c("s", "overwrite"))
  • Importar punto de desfogue/cierre de referencia a GRASS-GIS.
execGRASS("v.in.ogr", input="Estacion_cierre_paicol1.shp",
          output="Estacion_cierre", flag=c("overwrite"))
  • Después, dada la resolución del DEM, podrá notar que la estación de cierre no está quedando sobre la red principal del río, por ello, tendrá que tomar como referencia el punto de cierre y lo desplazará hasta encontrarse con la red principal tomando lo más central posible una coordenada.

De acuerdo con lo expresado, para realizar esta corrección se propone: (i) exporte el ráster de acumulaciones generado (F_accum.tif) y carguelo junto con la estación de cierre en un SIG (ArcGIS o QGIS) desplazando el punto hasta que se encuentre con la red (como lo vimos en clase) y tomando esa coordenada como referencia.

# Exportar los flujos de las acumulaciones
execGRASS("r.out.gdal", input="accum", format="GTiff",
          output="f_accum.tif",
          flag=c("overwrite"))
# Esta sería la coordenada de la estación sin colocarla en la línea de flujo. Es sólo para ver la diferencia
Coord_Estacion_cierre=Estacion_cierre@coords

# Aquí se pega la coordenada tomada cuando ya esta sobre la línea de acumulación de flujo. En este punto verifique que la estación esta sobre la red de drenaje usando el F_accum generado en unos pasos anteriores y use esa coordenada.
Coord_Estacion_cierre=c(4692843.59, 1830194.005)
  • Ya con la coordenada, se procede a delimitar la cuenca a partir de un punto de cierre.
execGRASS("r.water.outlet", input="drainage_dir", output="Cuenca_9377",
          coordinates=Coord_Estacion_cierre,flag=c("overwrite"))
  • Luego, se convierte el ráster de la cuenca delimitada en un formato vectorial.
execGRASS("r.to.vect", input="Cuenca_9377", output="Cuenca_9377", type="area",
          flag=c("s", "overwrite"))
  • Finalmente, se exporta la cuenca en formato vector.
execGRASS("v.out.ogr", input="Cuenca_9377", type="area", format="ESRI_Shapefile",
          output="Cuenca_9377_GRASS.shp",
          flag=c("c", "overwrite"))

2.2.1. Generar salidas gráficas (mapas)

Aquí se exponen dos maneras de generar los mapas (a cada una de las opciones le podrá agregar más cosas):

  • La primera a partir del comando plot.Esta es la opción base
# Cargar el archivo vector generado de la cuenca como shp
Cuenca_GRASS=shapefile("Cuenca_9377_GRASS.shp")

plot(Cuenca_GRASS,axes = TRUE,main = "Cuenca hidrográfica modelada",
     xlab = "Este (m.)",
     ylab = "Norte (m.)",
     las = 2, cex.axis=0.6, col="blue")
plot(Estacion_cierre, add = TRUE, col="red")
grid(col = "gray")
addnortharrow("topleft",scale = 0.4)
addscalebar(plotunit = "km",pos = "bottomright", unitcategory = "metric")

  • Y, la segunda, por medio de ggplot. Esta es la opción VIP
# Cargar el archivo vector generado de la cuenca como sf
Cuenca_GRASS=st_read("Cuenca_9377_GRASS.shp")
## Reading layer `Cuenca_9377_GRASS' from data source 
##   `C:\Users\maria.arenasb\OneDrive - Pontificia Universidad Javeriana\Javeriana\Scripts - R\Cuenca_9377_GRASS.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 4632052 ymin: 1789379 xmax: 4701492 ymax: 1881418
## Projected CRS: MAGNA-SIRGAS / Origen-Nacional
ggplot(data = Cuenca_GRASS) +
  geom_sf(fill = "greenyellow",color="black") +
  # coord_sf(crs = sf::st_crs(9377)) + xlab("Este") + ylab("Norte") +
  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 modelada", 
          subtitle = "Delimitación de cuencas hidrográficas") +
  annotation_scale() +
  annotation_north_arrow(location='tr')

Y con el fin de exportar cada uno de los mapas en diferentes formatos:

  • Para el caso de la salida generada con el comando plot en formato .pdf.
pdf("Mapa_plot.pdf", onefile=TRUE , height=12 , width=8 , paper="a4")
Cuenca_GRASS=shapefile("Cuenca_9377_GRASS.shp")

# Margen del mapa abajo, izquierda, arriba y derecha
par(mar = c(4, 4, 2.5, 2), xpd = FALSE)
plot(Cuenca_GRASS,axes = TRUE,main = "Cuenca hidrográfica modelada",
     xlab = "Este (m.)",
     ylab = "Norte (m.)",
     las = 2, cex.axis=0.6, col="blue")
plot(Estacion_cierre, add = TRUE, col="red")
grid(col = "gray")
addnortharrow("topleft",scale = 0.4)
addscalebar(plotunit = "km",pos = "bottomright", unitcategory = "metric")
dev.off()
## png 
##   2
  • Mientras que, para el caso de la salida generada con el comando ggplot en formato .jpeg.
jpeg("Mapa_ggplot.jpeg", quality = 500)
Cuenca_GRASS=st_read("Cuenca_9377_GRASS.shp")
## Reading layer `Cuenca_9377_GRASS' from data source 
##   `C:\Users\maria.arenasb\OneDrive - Pontificia Universidad Javeriana\Javeriana\Scripts - R\Cuenca_9377_GRASS.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 4632052 ymin: 1789379 xmax: 4701492 ymax: 1881418
## Projected CRS: MAGNA-SIRGAS / Origen-Nacional
par(mar=c(4,4,.5,.5))
ggplot(data = Cuenca_GRASS) +
  geom_sf(fill = "greenyellow",color="black") +
  # coord_sf(crs = sf::st_crs(9377)) + xlab("Este") + ylab("Norte") +
  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 modelada", 
          subtitle = "Delimitación de cuencas hidrográficas") +
  annotation_scale() +
  annotation_north_arrow(location='tr')
dev.off()
## png 
##   2