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.

otra forma de instalar los paquetes es usando la función Packages de la ventana izquierda o la función “install.packages” disponible en el Menú Tools de la barra superior de la ventana.

rm(list = ls())
require(xts)
require(lubridate)
require(dplyr)
require(ggplot2)
require(reshape2)
require(openxlsx)
require(raster) #Para que se puedan leer archivos ráster
require(ggspatial) #Para los datos espaciales en las gráficas
require (RSAGA) #Para que R articule con SAGA-GIS
require(rgdal)
require(sf)

Una vez instalados los paquetes requeridos, por favor indique con el comando “setwd” la ubicación de la CARPETA donde tiene los archivos con los que vamos a trabajar.Por favor recuerde que para correr este comando debe usar DOBLE SLASH. Es recomendable que cree una carpeta de trabajo para este curso y que esta carpeta no este guardada en la nube (OneDrive) para evitar problemas de sincronización por falta de internet durante las clases.

# Carpeta de trabajo
setwd("D://Javeriana//Scripts - R") #en mi caso esta es ubicación donde tengo los archivos de trabajo

DEM

A continuación iniciamos con la carga y visualización del DEM. Para este ejercicio puede cargar el raster completo de Colombia o el raster de su zona de influencia descargado de cualquiera de las plataformas trabajadas en el taller de ArcGIS.

dem <- raster("mascara4.tif") #Se carga el DEM, importante que esté dentro de la carpeta de trabajo
plot(dem) #si desea visualizar el archivo que asignó a la varibale, use la función plot.

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=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs"

crs(dem) #Es mucho más detallada

## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1
## +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["MAGNA_Colombia_Bogota",
##     BASEGEOGCRS["GCS_MAGNA",
##         DATUM["MAGNA",
##             ELLIPSOID["GRS_1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1,
##                     ID["EPSG",9001]]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]],
##     CONVERSION["Transverse Mercator",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",4.59620041666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-74.0775079166667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",1000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",1000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["easting",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["northing",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]

dem #Obtener información del DEM, aquí veo las coordenadas para establecer región de trabajo

## class      : RasterLayer 
## dimensions : 3024, 2327, 7036848  (nrow, ncol, ncell)
## resolution : 30.97393, 30.97393  (x, y)
## extent     : 751924.2, 824000.5, 721278.7, 814943.9  (xmin, xmax, ymin, ymax)
## 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 
## source     : mascara4.tif 
## names      : mascara4 
## values     : 0, 5310  (min, max)

Proyecciones Cartográficas y Sistemas de Referencia

La proyección cartográfica es el método que representa la superficie curva de la tierra sobre un plano mediante el uso de modelos matemáticos. Así mismo, los sitemas de referencia son el conjunto de convenciones y conceptos teóricos adecuadamente modelados que permiten definir en cualquier momento, la orientación, ubicación y escala de tres ejes coordenados (X, Y, Z). Si desea mayor detalle por favor consultar: https://www.igac.gov.co/sites/igac.gov.co/files/anexo_8._tipos_de_coordenadas_manejados_en_colombia_igac._ano_2004.pdf, https://docs.qgis.org/3.16/es/docs/gentle_gis_introduction/coordinate_reference_systems.html.

En Colombia, la entidad encargada de la actualización y validación de la información cartográfica es el Instituto Geográfico Agustin Codazzi -IGAC. Hasta el 2020, Colombia usaba una proyección cartográfica Gauß-Krüger (IGAC, 2004), datum MAGNA SIRGAS. Gauß-Krüger es una adaptación de la proyección transversa de Mercator, utiliza como superficie de referencia el área superficial de un cilindro transverso tangente a lo largo del meridiano de origen (meridiano central). Este sistema de proyección contaba con seis (6) puntos de origen y para actualizarnos a un sistema intenacional, era necesario tener un único punto de origen. Por esta razón en 2020 actualizamos el sistema de referencia a Land Administration Domain Model -LADM-, y un marco geográfico único. Sin emabargo, el datum seguirá siendo el datum MAGNA SIRGAS y la proyección se basará en la proyección cartográfica Transverse de Mercator.Para mayor detalle por favor consultar el documento “ABC Nueva Proyección Cartográfica para Colombia - Origen Nacional” disponible en https://origen.igac.gov.co/documentos.html

Si su raster (DEM) no se encuentra en coordenadas planas, puede usar la siguiente línea de código para realizar la reproyección y exportar el nuevo raster. Al usar estas líneas de código, no olvide eliminar el símbolo # que se encuentra al inció de cada renglón.

#dem_repro=projectRaster(dem,crs=("+init=epsg:4326")) #Si se tuviese que reproyectar
#writeRaster(dem_repro, filename="dem_repro.tif", format="GTiff", overwrite=TRUE) #Si ese DEM reproyectado se tuviera que exportar para luego ser cargado en otra herramienta

Adicionalmente, puede generar una gráfica inicial, simple, sólo para confirmar que el DEM quedó cargado con el siguiente código. Si el raster que desea visualizar es el DEM reproyectado por favor no olvide cambiar el nombre “dem” por “del_repro”

#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.

Preliminares en R para habilitar SAGA-GIS

Primero deberá descargar el programa SAGA de la página de internet https://sourceforge.net/projects/saga-gis/files/SAGA%20-%208/SAGA%20-%208.1.1/. y seguir las instrucciones de instalación. Se recomienda que se instale la herramienta en C://Program Files.

Para que R corra SAGA-GIS se debe direccionar el comando a que lea la ubicación del programa. Adicionalmente, es necesario habilitar la librería con el código que se muestra a continuación.Ejecute estas líneas ÚNICAMENTE cuando SAGA este instalado en el computador.

envi = rsaga.env(path="C://Program Files//SAGA")

## 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) #verifica la versión instalada. En este ejercicio recomiendo la 8.1.1

## [1] "8.1.1"

Si desea revisar la información del paquete, antes de continuar con su ejecución puede consultar la página https://cran.r-project.org/web/packages/RSAGA/RSAGA.pdf

Proyectar un DEM dentro de R y generar archivos raster de extensión .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("mascara4.tif")
#Se pregunta en qué proyección se encuentra el DEM
projection(dem)

## [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 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 que se encuentra en 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)

Corrección del DEM - Fill Sink

Para el análisis hidrológico es importante que el DEM no tenga celdas sin información pues la definición de la red de flujo (acumulación y dirección) podría estimarse de manera errónea. Para esto se usan las funciones de preprocesamiento (Fill), las cuales ayudan a interpolar la información de elevación faltante. En este ejercicio, usamos la función sink removal como se muestra a continuación. Si desea consultar otras opciones de mejora del DEM puede consultar el CRAN.

rsaga.sink.removal(in.dem="dem_3116_sgrd.sgrd", out.dem = "demsink_fill.sgrd", method = "fill",env = envi)

Cargar y reproyectar el DEM en extensión .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(".","Est_Paicol") #Cargar archivo en formato vector (punto) en extensión .shp en R

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\maria.arenasb\OneDrive - Pontificia Universidad Javeriana\Javeriana\Scripts - R", layer: "Est_Paicol"
## with 1 features
## It has 19 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 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +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 .x1=x, x2=y, x3=z

##      coords.x1 coords.x2 coords.x3
## [1,]  812725.8  763849.5       785

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).

#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 =812725.8, #copiar las coordenadas del punto aquí
                                TARGET_PT_Y =763849.5,
                                ELEVATION = "demsink_fill.sgrd",
                                AREA = "Cuenca.sgrd",
                                METHOD = 0),env = envi) #Observar que se puso como método el 0

Es así como, se convierte la cuenca en un archivo vector de extensión .shp generada por herramientas de SAGA-GIS desde R.

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

Finalmente, generamos una salida gráfica (mapa) de la cuenca delimitada.

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

#Mapa
ggplot(data = Cuenca_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", 
          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)

Actividad de Clase

Incluir en el mapa la visualización del punto de salida usado para la delimitación
Realizar un plot que contenga el DEM y la cuenca delimitada.

Delimitación Cuencas 1

María Cristina Arenas

2022-06-02

DEM