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

Los procesos que se desarrollarán a continuación están realizados con información que podrá descargar AQUÍ (haga click sobre la palabra, es un hipervínculo). Note que únicamente requieren como entradas un DEM y la referencia del punto de cierre (para el caso es una estación de monitoreo limnimétrica).

1. Pre-procesamientos/Pre-liminares

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 para esta actividad.

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

Ahora, se recomienda guardar el Script en el que se está trabajando (sí hasta ahora no se ha hecho) y, dado que ya se definió la carpeta de trabajo, por defecto al salvar el progreso será direccionado por defecto a esta, lo único que se debe hacer es colocar un nombre.

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 ello, cuenta con varias opciones para hacerlo, a continuación se muestra la más cómoda. Para realizarlo, simplemente tendrá que 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 logró 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

Como ya se mencionó, este proceso se realiza con los datos que puede descargar AQUÍ (haga click sobre la palabra, es un hipervínculo). Recuerde, en lo posible, debe organizarlos dentro de la carpeta de trabajo. Asimismo, tenga en cuenta que el DEM ha tenido un pre-procesamiento, en este lo que se realizó fue recortar el ráster al tamaño del área de estudio o zona de influencia.

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

# Note que antes del archivo se encuentran dos (2) carpetas, a esto hace referencia
# lo descrito anteriormente en el direccionamiento de la carpeta de trabajo, los
# archivos se están llamando desde carpetas internas dentro de la carpeta de trabajo

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)

A continuación, se asegura y/o 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á todo 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. Sin embargo, 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)

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

# Exportar vector
writeOGR(Estacion_cierre, dsn = "C:/Users/ASUS/Downloads/prueba3", 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)

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/ e instalar en archivos de programa.

A continuación, se detalla uno a uno los pasos para llevar a cabo procesos dentro de GRASS-GIS con el objeto de delimitar una cuenca hidrográfica que, para el caso, será de ejemplo la Subzona hidrográfica (SZH) del Río Güejar ubicada en el departamento del Meta, de Colombia.

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
  • Ahora, se debe iniciar GRASS-GIS desde R, para eso tiene que seguir:
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
          gisDbase = "C:/Users/ASUS/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/ASUS/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 de la Location. 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 modelada, 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 que se vuelve algo manual se proponen dos (2) opciones, estas son: la primera es que exporte el ráster de acumulaciones generado y lo cargue junto con la estación de cierre en un SIG de su manejo (ejemplo, ArcMap de ArcGIS o QGIS) desplazando el punto hasta que se encuentre con la red y tomando la coordenada.
# 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 es sólo para ver la diferencia
Coord_Estacion_cierre=Estacion_cierre@coords

# Aquí se pega la coordenada tomada. 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(4924940.005874055, 1886190.791635942)
Coord_Estacion_cierre=c(4692843.59, 1830194.005)
  • Ya con la coordenada, se procede a modelar 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"))
  • Seguido, se convierte en la salida ráster de la cuenca modelada 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 generada.
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.
# 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.
# 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\ASUS\Downloads\prueba3\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\ASUS\Downloads\prueba3\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