Geomorfometría con Whitebox
Daniela Ballari dballari@gmail.com
Universidad del Azuay - (Cuenca-Ecuador)
Introduccion
Este material ha sido adaptado de https://vt-hydroinformatics.github.io/rgeowatersheds.html
En esta actividad, vamos a explorar cómo trabajar con datos raster en R mientras calculamos métricas relevantes desde el punto de vista hidrológico utilizando el paquete R llamado Whitebox Tools. Whitebox Tools es muy potente y tiene un conjunto extenso de herramientas, pero no está disponible en CRAN. Debes instalarlo con la línea de código que está comentada en la parte superior del siguiente fragmento de código.
The following activity is available as a template github repository at the following link: https://github.com/VT-Hydroinformatics/13-Geospatial-Raster-Hydro.git
For more: https://geocompr.robinlovelace.net/spatial-class.html#raster-data
To read in detail about any of the WhiteboxTools used in this activity, check out the user manual: https://jblindsay.github.io/wbt_book/intro.html
Instalar/Cargar paquetes y datos necesarios:
library(tidyverse)
library(raster)
library(sf)
library(whitebox)
library(tmap)
library(stars)
#If installing/using whitebox for the first time
# install_whitebox()
whitebox::wbt_init()
theme_set(theme_classic())Whitebox tools
El funcionamiento de Whitebox Tools puede ser un poco complicado al principio. Puedes desear pasar objetos de R a estas funciones y recibir objetos de R como resultado, pero ese no es el enfoque que siguen.
Básicamente, para la entrada, le proporcionas a la función de Whitebox Tools el nombre del archivo que contiene los datos de entrada. Para la salida, le indicas cómo se debe llamar el archivo de salida.
La función de Whitebox Tools genera los datos calculados en tu directorio de trabajo o en el directorio que le indiques en el argumento de salida.
Esto significa que si deseas hacer algo con la salida de la función de Whitebox Tools, debes leerla por separado.
Cargar el DEM
Primero, configuraremos tmap ya sea como “map” o “view” (interactivo), dependiendo de cómo queremos ver nuestros mapas.
¿Qué significa DEM y qué muestra? ¿Qué significa que el DEM tiene “5, 30 o 90 metros”?
Utilizaremos la función raster() para cargar el DEM. Indicamos que el sistema de coordenadas es WGS84 configurando el argumento crs igual a ‘+init=EPSG:4326’, donde 4326 es el número EPSG para WGS84.
A continuación, un artefacto de la salida del DEM para este análisis es que hay muchas celdas erróneas alrededor del borde que no pertenecen al DEM. Si creamos un mapa con ellas, realmente distorsiona la escala. Por lo tanto, vamos a establecer cualquier valor de elevación por debajo de 1500 pies como NA.
tmap_mode("view")## tmap mode set to interactive viewingdem <- raster("datos/brushDEMsm_5m.tif", crs = '+init=EPSG:4326')
dem[dem < 1500] <- NA
writeRaster(dem, "datos/brushDEMsm_5m_crs.tif", overwrite = TRUE) Visualizar DEM
Ahora vamos a representar el Modelo Digital de Elevación (DEM). Proporciona tm_shape con el raster y luego visualízalo con tm_raster. Le diremos a tmap que la escala en el raster es continua, qué paleta de colores usar, si mostrar o no la leyenda y luego agregaremos una barra de escala.
tm_shape(dem)+
tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
tm_scale_bar()Hillshade
Para generar un hillshade (sombreado de pendientes) a partir del Modelo Digital de Elevación (DEM) usando la función wbt_hillshade de Whitebox Tools en R, puedes seguir estos pasos:
El hillshade (sombreado de pendientes) es una excelente manera de resaltar las características del paisaje en un Modelo Digital de Elevación (DEM) al simular la iluminación del terreno con un sol sintético. Controlar el ángulo del sol y la dirección desde la que brilla te permite destacar diferentes aspectos del paisaje.
En el siguiente fragmento de código, leeremos el DEM y generaremos un hillshade con wbt_hillshade(). Luego, leeremos el hillshade de salida con raster() y crearemos un mapa con tmap. Utilizaremos la paleta “Greys” invertida agregando un signo negativo delante de ella para que el hillshade se vea bien.
wbt_hillshade(dem = "datos/brushDEMsm_5m_crs.tif",
output = "datos/brush_hillshade.tif",
azimuth = 300,
altitude = 45)
hillshade <- raster("datos/brush_hillshade.tif")
tm_shape(hillshade)+
tm_raster(style = "cont", palette = "-Greys", legend.show = FALSE)+
tm_scale_bar()Missing data
Aplicar este paso para cuando se trabaje con los datos de Ecuador
# wbt_fill_missing_data(input = "datos/ASTER_32717_2_crs.tif",
# output = "datos/ASTER_32717_2_crs_f.tif",
# filter = 3)
#
# dem_filter <- raster("datos/ASTER_32717_2_crs_f.tif")
#
# tm_shape(dem_filter)+
# tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
# tm_scale_bar()Parámetros de geomorfometria
Hillshade
wbt_hillshade(dem = "datos/brushDEMsm_5m_crs.tif",
output = "datos/brush_hillshade.tif",
azimuth = 300,
altitude = 45)
hillshade <- raster("datos/brush_hillshade.tif")
tm_shape(hillshade)+
tm_raster(style = "cont", palette = "-Greys", legend.show = FALSE)+
tm_scale_bar()Slope
wbt_slope(dem = "datos/brushDEMsm_5m_crs.tif",
output = "datos/slope.tif",
units = "degrees")
slope <- raster("datos/slope.tif")
tm_shape(hillshade)+
tm_raster(style = "cont", palette = "-Greys", legend.show = FALSE)+
tm_scale_bar()+
tm_shape(slope)+
tm_raster(style = "cont",legend.show = TRUE)+
tm_scale_bar()Al igual que establecimos valores no deseados como NA anteriormente, también podemos utilizar esta técnica como una herramienta de visualización. Crearemos un nuevo ráster en el que las áreas con una pendiente de menos de 30 grados se establecerán como NA. Cuando lo tracemos, mostrará la pendiente solo donde la pendiente es mayor a 30.
En su esencia, esto es una reclasificación. Puedes utilizar la misma estrategia para clasificar las pendientes en categorías. Por ejemplo, puedes asignar un valor de 1 a las pendientes de 0 a 60 por ciento para representar una pendiente baja y un valor de 2 a las pendientes mayores al 60 por ciento para representar una pendiente alta, o puedes crear cualquier cantidad de categorías. Esto puede ser muy útil para encontrar ubicaciones que cumplan con varios criterios. Por ejemplo, puedes utilizar esta técnica para identificar áreas adecuadas para una especie de árbol en particular, aves o sasquatches, si así lo deseas.
slope2 <- slope
slope2[slope2 > 30] <- NA
tm_shape(slope2)+
tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
tm_scale_bar()DEM para hidrología
Preparar un Modelo Digital de Elevación (DEM) para análisis hidrológicos es una etapa importante, ya que muchos de los algoritmos hidrológicos se basan en el flujo de agua siguiendo la pendiente del terreno a partir de las elevaciones en las celdas del DEM. Si a lo largo de una ruta de flujo no hay ninguna celda más baja que una ubicación en particular, el algoritmo se detendrá allí, lo que se conoce como una depresión o sumidero. Es fundamental eliminar estos sumideros antes de realizar análisis hidrológicos.
Correctamente, para tratar con las depresiones en un Modelo Digital de Elevación (DEM), puedes utilizar la técnica de “relleno”. Esto implica elevar las elevaciones de las celdas de las depresiones de manera gradual hasta que la depresión se llene y el agua pueda fluir cuesta abajo. El proceso de relleno garantiza que no haya áreas donde el flujo de agua se detenga debido a depresiones locales en el terreno.
El relleno de depresiones se puede realizar utilizando la función wbt_fill_depression. Esta función identificará las depresiones y ajustará las elevaciones de las celdas para eliminarlas, permitiendo que el flujo de agua sea continuo a lo largo del terreno.
Luego se tratan las depresiones mediante la técnica de “brecha” implica bajar el lado de la característica que está bloqueando el flujo para permitir que el agua fluya cuesta abajo. Este enfoque se utiliza para simular la eliminación de obstáculos artificiales o naturales que impiden el flujo de agua en una cuenca.
Para “brechar” las depresiones en un DEM, puedes utilizar la función wbt_breach_depressions. Esta función identificará las depresiones y realizará brechas en los puntos más bajos de las depresiones para permitir el flujo continuo del agua.
Primero, usaremos wbt_breach_depressions_least_cost() para brechar las depresiones, lo que disminuirá la elevación de las celdas que están bloqueando las depresiones. También definiremos una distancia máxima (dist) para buscar un lugar donde brechar la depresión y decidiremos si llenar cualquier depresión que quede después de realizar la operación.
Luego, para limpiar cualquier depresión restante, utilizaremos wbt_fill_depressions_wang_and_liu(). Es importante asegurarse de que esta función se aplique al resultado de la función de brecha de depresiones, no al DEM original.
wbt_breach_depressions_least_cost(
dem = "datos/brushDEMsm_5m_crs.tif",
output = "datos/bmstationdem_breached.tif",
dist = 5,
fill = TRUE)
wbt_fill_depressions_wang_and_liu(
dem = "datos/bmstationdem_breached.tif",
output = "datos/bmstationdem_filled_breached.tif"
)Visualiza los “filled sinks” y “reached depression”
Para visualizar las áreas donde se realizaron cambios en el DEM debido a las operaciones de llenado y brecha de depresiones, puedes restar el DEM preprocesado con las depresiones llenadas y brechas realizadas del DEM original. Luego, podrás observar áreas donde no ocurrieron cambios (valores iguales a cero), áreas que se llenaron (valores positivos) y áreas donde se disminuyó la elevación para “brechar” una depresión (valores negativos).
Para resaltar las áreas donde ocurrieron cambios, puedes establecer todas las celdas que sean iguales a cero como NA (No Aplicable) y luego representarlas sobre el hillshade para una mejor visualización. Esto te permitirá identificar claramente las ubicaciones donde se realizaron cambios en el DEM.
filled_breached <- raster("datos/bmstationdem_filled_breached.tif")
crs(filled_breached)<- crs(dem)
## What did this do?
difference <- dem - filled_breached
min(filled_breached)## Warning in min(): Nothing to summarize if you provide a single RasterLayer; see
## cellStats## class : RasterLayer
## dimensions : 522, 855, 446310 (nrow, ncol, ncell)
## resolution : 3.925595e-05, 3.925595e-05 (x, y)
## extent : -80.49738, -80.46382, 37.23591, 37.2564 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : bmstationdem_filled_breached.tif
## names : bmstationdem_filled_breacheddifference[difference == 0] <- NA
tm_shape(hillshade)+
tm_raster(style = "cont",palette = "-Greys", legend.show = FALSE)+
tm_scale_bar()+
tm_shape(difference)+
tm_raster(style = "cont",legend.show = TRUE)+
tm_scale_bar()## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.Pasos para delinear una cuenca
El proceso de delimitación de cuencas hidrográficas es una tarea
importante en análisis geoespaciales y en hidrología. Antes de
sumergirnos en los detalles, es útil comprender cómo funciona el proceso
en su conjunto. Aunque se presenta de manera lineal, en realidad
requiere cierta planificación para obtener los resultados deseados.
Utilizaremos la función wbt_watershed() de
Whitebox Tools para delimitar nuestras cuencas hidrográficas, pero hay
algunos pasos previos que debemos realizar.
Generación de la cuadrícula de dirección D8: Esta cuadrícula especifica en qué dirección fluirá el agua desde cada celda. Para que la delimitación de cuencas funcione correctamente, debemos preparar nuestro Modelo Digital de Elevación (DEM) para eliminar fosas o depresiones sin salidas. Esto se hace en dos etapas: primero, llenamos las fosas de una sola celda y luego abrimos las depresiones más grandes. Es importante recordar que debes llenar las depresiones de una sola celda y luego pasar el DEM resultante de ese proceso a la función de apertura de depresiones. Una vez hecho esto, puedes ejecutar la función de dirección D8 para generar tu cuadrícula de dirección.
Creación de puntos de vertido: Los puntos de vertido DEBEN estar ubicados en celdas que forman parte de la red de corrientes, según lo define tu DEM. Si incluso están a una celda de distancia de las celdas de alta acumulación de flujo que denotan la corriente en tu DEM, obtendrás solo una extraña y pequeña porción de una cuenca. Para crear los puntos de vertido, primero debes obtener las coordenadas de tus puntos de vertido. Puedes obtenerlas de Google Earth, utilizar un GPS en el campo o tal vez tienes coordenadas para una ubicación/estructura conocida que define tu “salida”. Una vez que tengas las coordenadas, debes asegurarte de que estén en la red de corrientes. Para hacerlo, usaremos la función de ajuste de puntos de vertido de Jensen en Whitebox Tools. Esta herramienta toma tus puntos y busca dentro de un área definida la red de corrientes, luego mueve los puntos a la red de corrientes. Esto requiere tus puntos de vertido y un archivo de ráster de corrientes. Entonces, necesitas crear el archivo de ráster de corrientes, que se hace utilizando la función de extracción de corrientes en Whitebox Tools. Esta función requiere una cuadrícula de acumulación de flujo D8, por lo que primero debes crearla.
Nota sobre unidades: Asegúrate de conocer las unidades de distancia en tus datos (grados - Latitud / Longitud o metros - UTM). SI está en grados decimales debes especificar cuánto debe buscar la función de puntos de vertido en grados. Si estuvieras utilizando datos UTM, esa cifra sería en metros, y si estuvieras en un sistema de proyección estatal (en VA), esa cifra sería en pies.
Proceso resumido:
Lee el DEM
Llena las fosas de una sola celda y abre las depresiones más grandes
Crea las cuadrículas de acumulación de flujo D8 y de dirección D8
Lee los puntos de vertido
Crea el ráster de corrientes
Ajusta los puntos de vertido al ráster de corrientes
Ejecuta la función de delimitación de cuencas
D8 Flow Accumulation
El análisis hidrológico D8 flow accumulation (acumulación de flujo D8) es una técnica importante para determinar cómo fluye el agua a través del terreno en función de la elevación del terreno en las celdas circundantes. En el algoritmo D8, el flujo de agua se puede dirigir en una de las 8 direcciones posibles, como se muestra en la figura que mencionaste.
El resultado de este análisis es un mapa de flujo de acumulación, que muestra cuántas celdas drenan hacia cada celda en el DEM. Es decir, el valor de una celda en el mapa de flujo de acumulación indica cuántas celdas contribuyen al flujo de agua en esa ubicación específica. Este tipo de mapa es útil para identificar claramente las áreas donde se acumula el flujo de agua y resalta los cursos de agua y los arroyos en el paisaje.
La función wbt_d8_flow_accumulation() toma como entrada el DEM o un archivo de dirección de flujo. Pasaremos nuestro DEM preprocesado con las depresiones llenadas y brechas realizadas. La salida predeterminada es el número de celdas que drenan hacia cada celda, pero también puedes elegir un área de contribución específica.
Visualizaremos nuestra salida trazando el logaritmo de la cuadrícula de acumulación de flujo D8 sobre el hillshade con una opacidad del 0.5 utilizando el parámetro “alpha” en tm_raster. Mapear con el hillshade nos ayuda a ver la acumulación de flujo en el contexto del paisaje. Trazar los valores en escala logarítmica nos permite ver las diferencias en la acumulación de flujo, ya que los valores altos son mucho más altos que los valores bajos en una cuadrícula de acumulación de flujo.
wbt_d8_flow_accumulation(input = "datos/bmstationdem_filled_breached.tif",
output = "datos/D8FA.tif")
d8 <- raster("datos/D8FA.tif")
# d8 <- raster("datos/D8FA.tif")
tm_shape(hillshade)+
tm_raster(style = "cont",palette = "-Greys", legend.show = FALSE)+
tm_shape(log(d8))+
tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE, alpha = .5)+
tm_scale_bar()D infinity flow accumulation
Otro método para calcular la acumulación de flujo es el algoritmo D infinity. Este funciona de manera similar al algoritmo D8, pero con algunas diferencias importantes.
Con D infinity, la dirección del flujo para cada celda puede ser cualquier ángulo, lo que permite una variedad infinita de posibilidades. Si el ángulo no apunta directamente hacia una de las celdas vecinas, el flujo desde la celda focal puede DIVIDIRSE entre las celdas vecinas. Veremos la diferencia resultante cuando observemos la salida de este método.
La función para calcular la acumulación de flujo utilizando el método D infinity es wbt_d_inf_flow_accumulation(), y al igual que la función D8, puede tomar como entrada un archivo de dirección de flujo o un DEM. En este caso, le proporcionaremos nuestro DEM preprocesado con las depresiones llenadas y brechas realizadas.
Al igual que con los datos D8, representaremos los valores de acumulación de flujo en escala logarítmica superpuestos en un hillshade con una opacidad del 50%.
Las diferencias notables entre los resultados de D8 y D infinity pueden incluir la forma en que se manejan los flujos en terrenos con ángulos y cómo se distribuye el flujo en las celdas vecinas. Las diferencias en los resultados pueden deberse a las características específicas del terreno y las características hidrológicas.
La elección de cuál método representa mejor la realidad depende de la naturaleza del terreno y los objetivos del análisis hidrológico. Ambos métodos tienen sus aplicaciones y ventajas en diferentes contextos. D8 tiende a simplificar el flujo en una de las 8 direcciones, lo que puede ser adecuado para terrenos más uniformes, mientras que D infinity permite un flujo más diverso en cualquier dirección, lo que puede ser más preciso en terrenos irregulares. La elección del método dependerá de la precisión requerida y las características del terreno en tu área de estudio.
wbt_d_inf_flow_accumulation("datos/bmstationdem_filled_breached.tif",
"datos/DinfFA.tif")
dinf <- raster("datos/DinfFA.tif")
tm_shape(hillshade)+
tm_raster(style = "cont",palette = "-Greys", legend.show = FALSE)+
tm_shape(log(dinf))+
tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE, alpha = 0.5)+
tm_scale_bar()## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.Map Stream Network
Una cosa interesante y a menudo muy útil que podemos hacer con las cuadrículas de acumulación de flujo que calculamos es mapear la red de arroyos en una cuenca. Si observamos cualquiera de las cuadrículas de acumulación de flujo, veremos que los valores más altos se encuentran en los arroyos. Por lo tanto, si determinamos el valor de acumulación de flujo en el lugar más alto de la red de arroyos con un flujo constante, podemos establecer todas las celdas con una acumulación de flujo inferior a eso como NO y solo tendremos celdas que estén en el arroyo.
A menudo, queremos que nuestra red de arroyos se represente como líneas, por lo que luego debemos convertir esa cuadrícula en un formato vectorial.
Whitebox Tools tiene dos funciones útiles que nos permiten hacer
esto: wbt_extract_streams() crea una
cuadrícula de la red de arroyos utilizando un umbral de acumulación de
flujo que le proporcionas. Toma como entrada una cuadrícula de
acumulación de flujo D8.
Luego, wbt_raster_streams_to_vector()
tomará la salida de wbt_extract_streams()
y un archivo de dirección de flujo D8 y generará un shapefile de tu red
de arroyos.
A continuación, extraemos los arroyos, generamos un archivo de dirección de flujo D8 y luego convertimos los arroyos de la cuadrícula a formato vector. Luego representamos los arroyos sobre el hillshade.
wbt_extract_streams(flow_accum = "datos/D8FA.tif",
output = "datos/raster_streams.tif",
threshold = 500)
wbt_d8_pointer(dem = "datos/bmstationdem_filled_breached.tif",
output = "datos/D8pointer.tif")
wbt_raster_streams_to_vector(streams = "datos/raster_streams.tif",
d8_pntr = "datos/D8pointer.tif",
output = "datos/streams.shp")
streams <- st_read("datos/streams.shp")## Reading layer `streams' from data source
## `G:\Mi unidad\DOCENCIA\MAESTRIA\2023\Recursos_Naturales\Sesion2B\practica_geomorfometria\datos\streams.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 341 features and 2 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -80.49736 ymin: 37.23601 xmax: -80.464 ymax: 37.25619
## CRS: NAtm_shape(hillshade)+
tm_raster(style = "cont",palette = "-Greys", legend.show = FALSE)+
tm_shape(streams)+
tm_lines(col = "blue")+
tm_scale_bar()## Warning: Currect projection of shape streams unknown. Long-lat (WGS84) is
## assumed.Puntos de vertido
El proceso de selección de puntos de vertido es fundamental para delimitar las cuencas hidrográficas de manera precisa. Estos puntos son las ubicaciones para las cuales delinearemos nuestras cuencas. Es crucial que estos puntos estén en la red de corrientes de cada cuenca. Si los puntos están incluso una celda fuera de la red de corrientes, no obtendrás una cuenca válida. En su lugar, obtendrás una pequeña franja que muestra el área que drena hacia ese único punto en el paisaje.
Incluso con ubicaciones GPS altamente precisas, todavía debemos verificar que nuestros puntos de vertido estén en la red de corrientes, porque el Modelo Digital de Elevación (DEM) podría no coincidir perfectamente con las ubicaciones de los puntos.
Afortunadamente, existe una función en Whitebox Tools que se asegura
de que nuestros puntos estén en la red de corrientes. La función
wbt.jenson_snap_pour_points() busca dentro
de una distancia definida desde los puntos que le proporcionas el arroyo
más cercano y luego mueve los puntos a esas ubicaciones. Para usar esta
función, primero debemos crear una red de corrientes.
El proceso para configurar nuestros puntos de vertido se verá de la siguiente manera:
Crear un dataframe con puntos de vertido.
Convertir el dataframe en un shapefile.
Escribir el shapefile en nuestro directorio de datos.
Mover los puntos con la función de ajuste de puntos de vertido.
ppoints <- tribble(
~Lon, ~Lat,
-80.482778, 37.240504,
-80.474464, 37.242990,
-80.471506, 37.244512
)
ppointsSP <- SpatialPoints(ppoints, proj4string = CRS("+proj=longlat +datum=WGS84"))
shapefile(ppointsSP, filename = "datos/pourpoints.shp", overwrite = TRUE)
tmap_mode("plot")## tmap mode set to plottingtm_shape(hillshade)+
tm_raster(style = "cont",palette = "-Greys", legend.show = FALSE)+
tm_shape(ppointsSP)+
tm_dots(col = "red", size = 2)
Finalmente, utilizaremos la función de ajuste de puntos de vertido de Jenson para mover los puntos de vertido a su ubicación correcta.
El parámetro snap_dist le indica a la
función a qué distancia buscar una corriente. Las unidades de los
archivos que estamos utilizando son grados decimales (o metros si
trabajas en UTM), por lo que debemos tener cuidado aquí. Utiliza un
valor de 0.0005 (para grados), que equivale a aproximadamente 50 metros.
Si pusieras 50, estarías buscando en 50 grados de latitud y longitud (lo
hice al hacer esta actividad y hubo muchos bloqueos).
Una vez que obtengas las corrientes y los puntos de vertido ajustados, léelos en tu entorno de R y represéntalos para asegurarte de que los puntos de vertido estén en las corrientes.
wbt_jenson_snap_pour_points(pour_pts = "datos/pourpoints.shp",
streams = "datos/raster_streams.tif",
output = "datos/snappedpp.shp",
snap_dist = 0.0005) #careful with this! Know the units of your data
pp <- shapefile("datos/snappedpp.shp")
pp <- st_read("datos/snappedpp.shp")## Reading layer `snappedpp' from data source
## `G:\Mi unidad\DOCENCIA\MAESTRIA\2023\Recursos_Naturales\Sesion2B\practica_geomorfometria\datos\snappedpp.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -80.48272 ymin: 37.2406 xmax: -80.47149 ymax: 37.24449
## Geodetic CRS: GCS_unknownpp <- st_set_crs(pp, st_crs("+init=EPSG:4326"))
#
streams <- raster("datos/raster_streams.tif")
streams <- st_read("datos/streams.shp")## Reading layer `streams' from data source
## `G:\Mi unidad\DOCENCIA\MAESTRIA\2023\Recursos_Naturales\Sesion2B\practica_geomorfometria\datos\streams.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 341 features and 2 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -80.49736 ymin: 37.23601 xmax: -80.464 ymax: 37.25619
## CRS: NAstreams <- st_set_crs(streams, st_crs("+init=EPSG:4326"))
tm_shape(hillshade)+
tm_raster(style = "cont",palette = "-Greys", legend.show = FALSE)+
tm_shape(streams)+
tm_lines(col = "blue")+
tm_shape(pp)+
tm_dots(col = "red", size=2)
Delinear cuenca hidrográfica
¡Ahora estamos listos para delimitar nuestras cuencas hidrográficas!
Utiliza la función wbt_watershed(), que
toma como entrada un archivo de puntero D8 (d8_pntr) y nuestros puntos
de vertido ajustados (pour_pts). Esta función generará un ráster en el
que cada cuenca estará poblada con un valor único y todas las demás
celdas serán NA.
Luego, lee los resultados de esta función de nuevo y represéntalos sobre el hillshade con alpha establecido en 0.5 para ver qué hizo.
wbt_watershed(d8_pntr = "datos/D8pointer.tif",
pour_pts = "datos/snappedpp.shp",
output = "datos/brush_watersheds.tif")
ws <- raster("datos/brush_watersheds.tif")
tm_shape(hillshade)+
tm_raster(style = "cont",palette = "-Greys", legend.show = FALSE)+
tm_shape(streams)+
tm_lines(col = "blue")+
tm_shape(ws)+
tm_raster(legend.show = TRUE, alpha = 0.5, style = "cat")+
tm_shape(pp)+
tm_dots(col = "red", size = 2)
Convert watersheds to shapefiles
Para realizar mapeos o análisis vectoriales, puede ser muy útil tener
tus cuencas hidrográficas representadas como polígonos. Para lograrlo,
utilizaremos el paquete “stars”.
st_as_stars() convierte nuestro ráster de
cuencas en un objeto que el paquete “stars” puede manejar y luego
st_as_sf() convierte el objeto “stars” en
un objeto vectorial “sf”. También necesitamos establecer
merge en TRUE, lo que indica a
st_as_sf que trate cada conjunto de celdas
con el mismo valor (nuestras cuencas) como su propia entidad.
Ahora podemos representar las versiones vectoriales de nuestras
cuencas y también utilizar filter() para
mostrar solo una a la vez o alguna combinación, en lugar de todas las
cuencas.
wsshape <- st_as_stars(ws) %>% st_as_sf(merge = T)
ws1shp <- wsshape %>% filter(brush_watersheds == "1")
tm_shape(hillshade)+
tm_raster(style = "cont",palette = "-Greys", legend.show = FALSE)+
tm_shape(ws1shp)+
tm_borders(col = "red")
Ahora un DEM de Ecuador
Este es un
Cargar el DEM
Presta atención al sistema de referencia de estos datos y a los datos que se están eliminando. ¿Por qué se está considerando 0?
tmap_mode("view")## tmap mode set to interactive viewingdem <- raster("datos/ASTER_32717_2.tif", crs = '+init=EPSG:32717')
dem[dem < 0] <- NA
writeRaster(dem, "datos/ASTER_32717_2_crs.tif", overwrite = TRUE) Visualizar DEM
tm_shape(dem)+
tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
tm_scale_bar()## stars object downsampled to 998 by 1001 cells. See tm_shape manual (argument raster.downsample)Hillshade
Prueba distintas configuraciones de azimuth y altitud para lograr que los picos y los valles se representen adecuadamente.
wbt_hillshade(dem = "datos/ASTER_32717_2_crs.tif",
output = "datos/aster_hillshade.tif",
azimuth = 115,
altitude = 45)
hillshade <- raster("datos/aster_hillshade.tif")
tm_shape(hillshade)+
tm_raster(style = "cont", palette = "-Greys", legend.show = FALSE)+
tm_scale_bar()## stars object downsampled to 998 by 1001 cells. See tm_shape manual (argument raster.downsample)Missing data
El relleno de datos es fundamental para el DEM que se está utilizando. Compara este resultado con el raster anterior.
wbt_fill_missing_data(input = "datos/ASTER_32717_2_crs.tif",
output = "datos/ASTER_32717_2_crs_f.tif",
filter = 3)
dem_filter <- raster("datos/ASTER_32717_2_crs_f.tif")
tm_shape(dem_filter)+
tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
tm_scale_bar()## stars object downsampled to 998 by 1001 cells. See tm_shape manual (argument raster.downsample)Parámetros de geomorfometria
Hillshade
wbt_hillshade(dem = "datos/ASTER_32717_2_crs_f.tif",
output = "datos/aster_hillshade.tif",
azimuth = 115,
altitude = 45)
hillshade <- raster("datos/aster_hillshade.tif")
tm_shape(hillshade)+
tm_raster(style = "cont", palette = "-Greys", legend.show = FALSE)+
tm_scale_bar()## stars object downsampled to 998 by 1001 cells. See tm_shape manual (argument raster.downsample)Slope
wbt_slope(dem = "datos/ASTER_32717_2_crs_f.tif",
output = "datos/aster_slope.tif",
units = "degrees")
slope <- raster("datos/aster_slope.tif")
tm_shape(hillshade)+
tm_raster(style = "cont", palette = "-Greys", legend.show = FALSE)+
tm_scale_bar()+
tm_shape(slope)+
tm_raster(style = "cont",legend.show = TRUE)+
tm_scale_bar()## stars object downsampled to 998 by 1001 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 998 by 1001 cells. See tm_shape manual (argument raster.downsample)Prueba cambiar el valor límite para visualizar slope (30) y reflexiona sobre el resultado.
slope2 <- slope
slope2[slope2 < 30] <- NA
tm_shape(slope2)+
tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
tm_scale_bar()## stars object downsampled to 998 by 1001 cells. See tm_shape manual (argument raster.downsample)DEM para hidrología
Dilled sinks and reached depressions
wbt_breach_depressions_least_cost(
dem = "datos/ASTER_32717_2_crs_f.tif",
output = "datos/aster_breached.tif",
dist = 5,
fill = TRUE)
wbt_fill_depressions_wang_and_liu(
dem = "datos/aster_breached.tif",
output = "datos/aster_filled_breached.tif"
)filled_breached <- raster("datos/aster_filled_breached.tif")
difference <- dem - filled_breached
difference[difference == 0] <- NA
tm_shape(hillshade)+
tm_raster(style = "cont",palette = "-Greys", legend.show = FALSE)+
tm_scale_bar()+
tm_shape(difference)+
tm_raster(style = "cont",legend.show = TRUE)+
tm_scale_bar()## stars object downsampled to 998 by 1001 cells. See tm_shape manual (argument raster.downsample)
## stars object downsampled to 998 by 1001 cells. See tm_shape manual (argument raster.downsample)## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.