Cuaderno DEM
MARIA ERNANDA SANCHEZ CALDERON
2023-05-18
1. ¿Para qué es este cuaderno?
Este cuaderno ilustra varias funcionalidades para obtener, procesar y visualizar modelos digitales de elevación (MDE) en R. Su objetivo es ayudar a los estudiantes de Geomática Básica de la Universidad Nacional a familiarizarse con los MDE y las variables geomorfométricas.
Algunos consejos para escribir tu propio cuaderno:
Escribe y ejecuta cada fragmento paso a paso. En caso de errores, lee el mensaje, interpreta y trata de solucionarlo. Si el error persiste, consulta al profesor antes de la fecha límite. Una vez que se haya generado tu archivo HTML, es momento de publicarlo en RPubs. Después de publicar el cuaderno, revísalo, corrige la redacción deficiente y vuelve a publicarlo utilizando la opción de actualización. En caso de dificultades, por favor, consulta la documentación de R.
Setup
Los paquetes requeridos deben instalarse previamente desde la Consola.
## SETUP
# INSTALL THESE PACKAGES FROM THE CONSOLE, NOT FROM THIS CHUNK
#paquetes = c("rgdal","raster", "rgeos", "terra","elevatr","rasterVis", "rgl", "mapview")
#install.packages(paquetes)Es recomendable empezar a limpiar la memoria:
rm(list=ls())Ahora, cargaremos las librerias
library(rasterVis)## Warning: package 'rasterVis' was built under R version 4.2.3## Loading required package: latticelibrary(raster)## Warning: package 'raster' was built under R version 4.2.3## Loading required package: sp## Warning: package 'sp' was built under R version 4.2.3library(terra)## Warning: package 'terra' was built under R version 4.2.3## terra 1.7.29library(rgeos)## Warning: package 'rgeos' was built under R version 4.2.3## rgeos version: 0.6-3, (SVN revision 696)
## GEOS runtime version: 3.9.3-CAPI-1.14.3
## Please note that rgeos will be retired during October 2023,
## plan transition to sf or terra functions using GEOS at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html for details.
## GEOS using OverlayNG
## Linking to sp version: 1.6-0
## Polygon checking: TRUElibrary(rgdal)## Warning: package 'rgdal' was built under R version 4.2.3## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-6, (SVN revision 1201)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/ferxx/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/ferxx/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.##
## Attaching package: 'rgdal'## The following object is masked from 'package:terra':
##
## projectlibrary(elevatr)## Warning: package 'elevatr' was built under R version 4.2.3library(sf)## Warning: package 'sf' was built under R version 4.2.3## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUElibrary(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.0 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.1 ✔ tibble 3.1.8
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks terra::extract(), raster::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::select() masks raster::select()
## ✖ dplyr::symdiff() masks rgeos::symdiff()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorslibrary(mapview)## Warning: package 'mapview' was built under R version 4.2.3library(ggplot2)
library(leaflet)## Warning: package 'leaflet' was built under R version 4.2.3library(tmap)## Warning: package 'tmap' was built under R version 4.2.3library(RColorBrewer)Al ejecutar un fragmento, presta atención a cualquier mensaje que indique un error o advertencia sobre posibles conflictos.
Los conflictos pueden ser causados por funciones de diferentes paquetes que tienen los mismos nombres (por ejemplo, raster::extract y tidyr::extract). Para evitar problemas, debemos llamar a dichas funciones utilizando la forma larga, es decir, nombre_paquete::funcion.
Después de ejecutar un fragmento por primera vez, es una buena práctica ocultar los mensajes y las advertencias usando: {r message=FALSE, warning=FALSE}.
2. Introduccion
Los datos de elevación se utilizan en una amplia variedad de aplicaciones, incluyendo la visualización, la hidrología y la modelización ecológica, por ejemplo. Acceder a estos datos en R no ha tenido una única interfaz, se ha puesto a disposición a través de funciones en muchos paquetes o requiere acceso local a los datos. Esto ya no es necesario, ya que ahora existen varias API que proporcionan acceso programático a datos de elevación. El paquete elevatr se creó para estandarizar el acceso a datos de elevación desde API web.
Para acceder a datos de elevación en formato raster (por ejemplo, un MDE), el paquete elevatr utiliza los Azulejos de Terreno de Amazon Web Services. Exploraremos esta funcionalidad en este cuaderno.
Existen varias fuentes de modelos digitales de elevación, como la Misión Topográfica del Radar del Transbordador Espacial (SRTM), el Conjunto de Datos de Elevación Nacional (NED) del USGS, el MDE Global (GDEM) y otros. Cada uno de estos MDE tiene ventajas y desventajas para su uso. Antes de su cierre en enero de 2018, Mapzen combinó varias de estas fuentes para crear un producto de elevación de síntesis que utiliza los mejores datos de elevación disponibles para una región y un nivel de zoom determinados. Aunque Mapzen cerró, estos datos compilados por Mapzen siguen siendo accesibles a través de los Azulejos de Terreno en Amazon Web Services (AWS).
La entrada para get_elev_raster() puede ser un data.frame con columnas x (longitud) e y (latitud) como las dos primeras columnas, cualquier objeto sp o cualquier objeto raster, y devuelve una RasterLayer de los azulejos que se superponen a la caja delimitadora de la entrada. Si se recuperan varios azulejos, la salida resultante es una Raster Layer fusionada.
Uso de get_elev_raster() para acceder a los Azulejos de Terreno en AWS.
Como se mencionó, se debe proporcionar un data frame con columnas x e y, un objeto sp o un objeto raster como entrada, y el argumento src debe establecerse en “mapzen” (esto es lo predeterminado).
No hay diferencia en el uso de los tipos de datos de entrada sp y raster.
El data frame requiere un CRS (es decir, un sistema de referencia de coordenadas). El nivel de zoom (z) se establece en 9 de forma predeterminada (un compromiso entre la resolución y el tiempo de descarga). Puedes intentar usar un nivel de zoom más alto (por ejemplo, 10).
3. Obtener datos de elevacion para nuestro departamento
Primero, necesitamos cargar un shapefile o un geopackage que represente los municipios de nuestro departamento. En este cuaderno, usaré un shapefile descargado del geoportal del DANE.
Vamos a leer el shapefile utilizando una función proporcionada por el paquete sf:
munic <- st_read("C:/EJERCICIOS/EJERCICIO2/DATOS/MUNICIPIOS.gpkg")## Reading layer `MUNICIPIOS' from data source
## `C:\EJERCICIOS\EJERCICIO2\DATOS\MUNICIPIOS.gpkg' using driver `GPKG'
## Simple feature collection with 14 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.3595 ymin: 4.670201 xmax: -75.411 ymax: 5.475137
## Geodetic CRS: WGS 84Let’s check geometry, bounding box, CRS, and attributes of the munic object:
- si es necesario puede cambiar los nombres de las columnas o eliminar las que no se necesiten
# Eliminar columnas de un dataframe
borrar <- c("NL_NAME_2","VARNAME_2")
munic <- munic[ , !(names(munic) %in% borrar)]head(munic)## Simple feature collection with 6 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.1697 ymin: 4.751701 xmax: -75.6503 ymax: 5.370399
## Geodetic CRS: WGS 84
## ID_0 ISO NAME_0 ID_1 NAME_1 ID_2 NAME_2 TYPE_2 ENGTYPE_2
## 1 53 COL Colombia 25 Risaralda 838 Apía Municipio Municipality
## 2 53 COL Colombia 25 Risaralda 839 Balboa Municipio Municipality
## 3 53 COL Colombia 25 Risaralda 840 Belén de Umbría Municipio Municipality
## 4 53 COL Colombia 25 Risaralda 841 Dosquebradas Municipio Municipality
## 5 53 COL Colombia 25 Risaralda 842 Guática Municipio Municipality
## 6 53 COL Colombia 25 Risaralda 843 La Celia Municipio Municipality
## AREA AREA1 geom
## 1 255 255.15373 MULTIPOLYGON (((-76.1697 4....
## 2 149 148.75731 MULTIPOLYGON (((-76.0186 4....
## 3 175 175.30133 MULTIPOLYGON (((-75.986 5.2...
## 4 71 70.88232 MULTIPOLYGON (((-75.6566 4....
## 5 101 101.06658 MULTIPOLYGON (((-75.8978 5....
## 6 73 73.21514 MULTIPOLYGON (((-76.05933 4...Creemos un nuevo objeto con centroides de municipios:
# The task is to find a center point for each region
# We convert the simple features to spatial polygons
sp.munic = as_Spatial(munic)
# Now we use the *gCentroid* function from the rgeos package
centers <- data.frame(gCentroid(sp.munic, byid = TRUE))
centers$region <- row.names(sp.munic)
centers$munic <- sp.munic$NAME_2centers## x y region munic
## 1 -76.03789 5.121155 1 Apía
## 2 -75.98727 4.917795 2 Balboa
## 3 -75.91388 5.194985 3 Belén de Umbría
## 4 -75.68981 4.825198 4 Dosquebradas
## 5 -75.85704 5.312344 5 Guática
## 6 -76.08330 4.991792 6 La Celia
## 7 -75.88906 4.895382 7 La Virginia
## 8 -75.78663 4.931795 8 Marsella
## 9 -76.00821 5.368507 9 Mistrato
## 10 -75.73321 4.770035 10 Pereira
## 11 -76.19165 5.205294 11 Pueblo Rico
## 12 -75.75777 5.319783 12 Quinchía
## 13 -75.57197 4.832214 13 Santa Rosa de Cabal
## 14 -76.02013 5.027967 14 SantuarioAhora, descarguemos los datos de elevación para nuestro departamento:
# z denotes the zoom level of the data
# the higher the zoom the higher the spatial resolution
#elevation <- get_elev_raster(munic, z = 12)
elevation <- get_elev_raster(munic, z = 10)## Mosaicing & Projecting## Note: Elevation units are in meters.¿Qué es este elemento?
elevation## class : RasterLayer
## dimensions : 1532, 2051, 3142132 (nrow, ncol, ncell)
## resolution : 0.0006856682, 0.0006856682 (x, y)
## extent : -76.64062, -75.23432, 4.565542, 5.615986 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : file4508777a267f.tif
## names : file4508777a267fObserva algunas cosas sobre este MDE:
dimensiones: el “tamaño” del archivo en píxeles. nrow, ncol: el número de filas y columnas en los datos (imagina una hoja de cálculo o una matriz). ncells: el número total de píxeles o celdas que componen el raster. resolución: el tamaño de cada píxel (en grados decimales en este caso). Recordemos que 1 grado decimal representa aproximadamente 111.11 km en el ecuador. extensión: la extensión espacial del raster. Este valor estará en las mismas unidades de coordenadas que el sistema de referencia de coordenadas del raster. crs: la cadena del sistema de referencia de coordenadas para el raster. Este raster está en coordenadas geográficas con un datum de WGS 84.
En tu cuaderno, anota el tamaño de celda de tu raster (es decir, la resolución espacial), en metros. En este caso es de 0.0006856682, 0.0006856682 (x, y)
Ten en cuenta que elevation es una capa de raster temporal (es decir, solo existe en la memoria). Sin embargo, es posible escribir el MDE en tu directorio local utilizando la función writeRaster proporcionada por la biblioteca rgdal:
# uncomment if needed
# write to geotiff file (depends on rgdal)
# NOTE THAT z10 means Zoom Level 10
# CHANGE IT IF NEEDED
if (require(rgdal)) {
rf <- writeRaster(elevation, filename=file.path("C:/EJERCICIOS/EJERCICIO2/DATOS/RIS_DEM_Z10.tiff"), format="GTiff", overwrite=TRUE)
}En mi caso, los datos de elevación en el zoom 12 fueron de unos 26,4 MB. Cuando el zoom era 10, el tamaño del archivo se redujo a unos 7,44 MB.
Anote el tamaño de su DEM.
4 VISUALIZACION
Para fines de visualización, una resolución más baja acelera la tarea:
#This chunk is optional
#Use it only when zoom data was very high
elevation2 <- aggregate(elevation, fact=4, fun=mean)crs(elevation2) <- "+proj=longlat +datum=WGS84"Una buena paleta es clave para una correcta visualización. Probemos con este:
pal <- colorNumeric(c("forestgreen","yellow","tan","brown"), values(elevation2),
na.color = "transparent")Ahora, usaremos la biblioteca de folletos para ver los datos de elevación:
leaflet(munic) %>% addTiles() %>%
addPolygons(color = "white", weight = 2.0, smoothFactor = 0.5,
opacity = 0.25, fillOpacity = 0.15,
popup = paste("Region: ", munic$NAME_2, "<br>",
"Km2: ", munic$AREA1, "<br>")) %>%
addLabelOnlyMarkers(data = centers,
lng = ~x, lat = ~y, label = ~region,
labelOptions = labelOptions(noHide = TRUE, direction = 'top', textOnly = TRUE)) %>%
addRasterImage(elevation2, colors = pal, opacity = 0.9) %>%
addLegend(pal = pal, values = values(elevation2),
title = "Elevation data for Risaralda (mts)")Haga clic en cada ID de región para obtener el nombre y la extensión del municipio. Cambie las opciones de trazado para que su mapa sea más legible.
5. Reproyectar los datos de elevación
Cuando se trabaja con MDE, siempre es buena idea utilizar coordenadas de mapa en lugar de coordenadas geográficas. Esto se debe a que, en coordenadas geográficas, las unidades de las dimensiones horizontales son grados decimales, PERO la unidad de la dimensión vertical es metros.
Para reproyectar los datos de elevación, llevamos a cabo dos pasos.
Primero, definimos el sistema de referencia de coordenadas (CRS) de destino:
#library(sp)
# WGS 84 / UTM zone 18N
spatialref <- CRS(SRS_string="EPSG:32618")Luego, reproyectamos el DEM usando la función projectRaster del paquete raster:
# Project Raster
# First, we create one raster object
# using projectExtent (no values are transferred)
pr3 <- projectExtent(elevation, crs=spatialref)
# Adjust the cell size to your cell size
res(pr3) <- 80
# now project
rep_elev <- projectRaster(elevation, pr3)¿Qué obtenemos?
rep_elev## class : RasterLayer
## dimensions : 1455, 1951, 2838705 (nrow, ncol, ncell)
## resolution : 80, 80 (x, y)
## extent : 317990.2, 474070.2, 504609.7, 621009.7 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : file4508777a267f
## values : -116.0929, 5278.495 (min, max)Ahora reproyectemos el SpatialPolygonsDataFrame que representa a los municipios de nuestro departamento.
(rep_munic = spTransform(sp.munic,spatialref))## class : SpatialPolygonsDataFrame
## features : 14
## extent : 349337, 454422.6, 516228.8, 605258.1 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## variables : 11
## names : ID_0, ISO, NAME_0, ID_1, NAME_1, ID_2, NAME_2, TYPE_2, ENGTYPE_2, AREA, AREA1
## min values : 53, COL, Colombia, 25, Risaralda, 838, Apía, Municipio, Municipality, 31, 30.5922403667703
## max values : 53, COL, Colombia, 25, Risaralda, 851, Santuario, Municipio, Municipality, 930, 929.547212212389Para evitar un dolor de cabeza, guardemos nuestro DEM. Asegúrese de indicar un nombre de ruta y un nombre de archivo que se adapte a sus necesidades (es decir, apropiado para su plataforma y su departamento).
# Write the reprojected DEM to disk
#uncomment and run the following line
#once you run the line, comment it again
writeRaster(rep_elev, filename="C:/EJERCICIOS/EJERCICIO2/DATOS/REP_RISA_DEM_Z10.tiff", datatype='INT4S', overwrite=TRUE)6. Estadísticas básicas de datos de elevación
Una exploración rápida de las estadísticas DEM.
Este fragmento “limpia” la elevación por debajo de 0.0:
rep_elev[rep_elev < 0.0] <- 0.0Ahora, vamos a crear un histograma con los valores de elevacion
# histograma
hist(rep_elev)
Podemos calcular un resumen de las estadísticas DEM:
# works for large files
promedio <- cellStats(rep_elev, 'mean')
minimo <- cellStats(rep_elev, 'min')
maximo <- cellStats(rep_elev, 'max')
desviacion <- cellStats(rep_elev, 'sd')Vamos a crear un vector de estadísticas unidimensional:
metricas <- c('mean', 'min', 'max', 'std')
valores <- c(promedio, minimo, maximo, desviacion)Ahora, el resultado
(df_estadisticas <- data.frame(metricas, valores))## metricas valores
## 1 mean 1552.790
## 2 min 0.000
## 3 max 5278.495
## 4 std 1052.7037. Obtención de variables geomorfométricas
Antes de continuar, asegúrese de haber entendido los conceptos básicos sobre: DEM, pendiente y orientación.
Primero, calcule la pendiente, el aspecto y el sombreado:
slope = terrain(rep_elev,opt='slope', unit='degrees')
aspect = terrain(rep_elev,opt='aspect',unit='degrees')
hill = hillShade(slope,aspect,40,315)Veamos a primera salida:
slope## class : RasterLayer
## dimensions : 1455, 1951, 2838705 (nrow, ncol, ncell)
## resolution : 80, 80 (x, y)
## extent : 317990.2, 474070.2, 504609.7, 621009.7 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : slope
## values : 0, 85.54764 (min, max)Ahora la segunda:
aspect## class : RasterLayer
## dimensions : 1455, 1951, 2838705 (nrow, ncol, ncell)
## resolution : 80, 80 (x, y)
## extent : 317990.2, 474070.2, 504609.7, 621009.7 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : aspect
## values : 0, 359.9998 (min, max)Para trazar la pendiente, es una buena idea reducir el tamaño de RasterLayer:
#This chunk is optional
#Use it only when zoom data was very high
slope2 <- aggregate(slope, fact=4, fun=mean)Puede ser conveniente proyectar el conjunto de datos de elevación:
slope3 <- projectRasterForLeaflet(slope2, "ngb")Vamos a crear una paleta de colores:
pal2 <- colorNumeric(c("#537188", "#CBB279", "#E1D4BB", "#EEEEEE"), values(slope3),na.color = "transparent")Ahora, el mapa:
leaflet(munic) %>% addTiles() %>%
addPolygons(color = "white", weight = 1.0, smoothFactor = 0.5,
opacity = 0.25, fillOpacity = 0.15,
popup = paste("Region: ", munic$NAME_2, "<br>",
"Km2: ", munic$AREA1, "<br>")) %>%
addLabelOnlyMarkers(data = centers,
lng = ~x, lat = ~y, label = ~region,
labelOptions = labelOptions(noHide = TRUE, direction = 'top', textOnly = TRUE)) %>%
addRasterImage(slope3, colors = pal2, opacity = 1.0) %>%
addLegend(pal = pal2, values = values(slope3),
title = "Slope for Risaralda (%)")Recordemos qué es el objeto slope2:
slope2## class : RasterLayer
## dimensions : 364, 488, 177632 (nrow, ncol, ncell)
## resolution : 320, 320 (x, y)
## extent : 317990.2, 474150.2, 504529.7, 621009.7 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : slope
## values : 0.1369883, 83.65871 (min, max)Tenga en cuenta que hay dos bibliotecas para administrar datos ráster en R: - raster y - terra
La biblioteca de Raster es bien conocida (pero, en algunas tareas, bastante lenta). La biblioteca terra es más reciente (pero más rápida).
Convertiremos un RasterLayer (es decir, un objeto raster) en un SpatRaster (es decir, un objeto terra):
(nslope2 <- as(slope2, "SpatRaster"))## class : SpatRaster
## dimensions : 364, 488, 1 (nrow, ncol, nlyr)
## resolution : 320, 320 (x, y)
## extent : 317990.2, 474150.2, 504529.7, 621009.7 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
## source(s) : memory
## name : slope
## min value : 0.1369883
## max value : 83.6587089La clasificación de pendientes es una forma común de entender la variabilidad del relieve. Aquí debes incluir una tabla con intervalos de pendiente, nombres de clase y descripciones, según lo definido por el Instituto Geográfico Agustín Codazzi (IGAC) de Colombia.
Realicemos una clasificación de pendientes basada en rangos comunes.
La tabla de clasificación de pendientes según el IGAC es la siguiente:
| Intervalo de Pendiente | Nombre de Clase | Descripción |
|---|---|---|
| 0% - 5% | Plana | Terreno prácticamente plano |
| 5% - 15% | Suave | Pendiente suave, ligeramente inclinada |
| 15% - 30% | Moderada | Pendiente moderada, terreno ondulado |
| 30% - 45% | Fuerte | Pendiente fuerte, terreno empinado |
| > 45% | Muy Fuerte | Pendiente muy fuerte, terreno escarpado |
Utilizando esta tabla, podemos realizar una clasificación de pendientes basada en estos rangos. Primero, crearemos una matriz de reclasificación.
m <- c(0, 3, 1, 3, 7, 2, 7,12,3, 12,
25, 4, 25, 50, 5, 50, 75, 6, 75, 90, 7)
(rclmat <- matrix(m, ncol=3, byrow = TRUE))## [,1] [,2] [,3]
## [1,] 0 3 1
## [2,] 3 7 2
## [3,] 7 12 3
## [4,] 12 25 4
## [5,] 25 50 5
## [6,] 50 75 6
## [7,] 75 90 7Ahora, la tarea de reclasificación:
# for values >= 0 (instead of > 0), do
slope_rec <- terra::classify(nslope2, rclmat, right=TRUE)Veamos el resultado
slope_rec## class : SpatRaster
## dimensions : 364, 488, 1 (nrow, ncol, nlyr)
## resolution : 320, 320 (x, y)
## extent : 317990.2, 474150.2, 504529.7, 621009.7 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
## source(s) : memory
## name : slope
## min value : 1
## max value : 78.visualización estática
Usaremos tmap para crear y guardar un mapa estático de clasificación de pendientes.
slope.map <- tm_shape(slope_rec) +
tm_raster(alpha = 0.6, style= "cat",
title="Slope") +
# Borders only
tm_shape(munic) + tm_borders(col = "darkgray", lwd = 0.8, alpha=0.8) +
tm_text(text = "NAME_2", size = 0.8, alpha=0.7) +
tm_scale_bar(text.size = 0.5) +
tm_compass(position = c("right", "top"), size = 2) +
tm_graticules(alpha = 0.2) +
tm_layout(outer.margins = 0.05, legend.position = c("left", "top"), title= 'Slope classes', title.position = c('right', 'top')) +
tm_credits("Data source: Mapzen & DANE", size=0.4) Veamos la salida
slope.map
Se pueden encontrar más consejos para trazar objetos ráster y vectoriales en [este enlace] (https://r-tmap.github.io/tmap-book/). Asegúrese de probar opciones adicionales a las incluidas en este cuaderno.
Ahora, tiempo de copia de seguridad:
## Note we are using letter paper size
tmap_save(slope.map, "C:/EJERCICIOS/EJERCICIO2/DATOS/risa_slope.png", width = 8.5, height =11, dpi=600)## Map saved to C:\EJERCICIOS\EJERCICIO2\DATOS\risa_slope.png## Resolution: 5100 by 6600 pixels## Size: 8.5 by 11 inches (600 dpi)Asegúrese de verificar la salida en su directorio de trabajo
8.Visualización interactiva
Ahora, queremos crear un mapa de pendientes interactivo.
Convirtamos el objeto sf en un objeto SpatVector:
(nmunic <- as(rep_munic, "SpatVector"))## class : SpatVector
## geometry : polygons
## dimensions : 14, 11 (geometries, attributes)
## extent : 349337, 454422.6, 516228.8, 605258.1 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## names : ID_0 ISO NAME_0 ID_1 NAME_1 ID_2 NAME_2
## type : <num> <chr> <chr> <num> <chr> <num> <chr>
## values : 53 COL Colombia 25 Risaralda 838 Apía
## 53 COL Colombia 25 Risaralda 839 Balboa
## 53 COL Colombia 25 Risaralda 840 Belén de Umbría
## TYPE_2 ENGTYPE_2 AREA AREA1
## <chr> <chr> <int> <num>
## Municipio Municipality 255 255.2
## Municipio Municipality 149 148.8
## Municipio Municipality 175 175.3Convirtamos el objeto raster en un objeto terra:
(nslope <- as(slope, "SpatRaster"))## class : SpatRaster
## dimensions : 1455, 1951, 1 (nrow, ncol, nlyr)
## resolution : 80, 80 (x, y)
## extent : 317990.2, 474070.2, 504609.7, 621009.7 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
## source(s) : memory
## name : slope
## min value : 0.00000
## max value : 85.54764Ahora, calculemos las categorías de pendientes agregadas por municipio:
## S4 method for signature 'SpatRaster,SpatVector'
munic.slope <- zonal(nslope, nmunic, fun=mean, as.raster=TRUE)Tenga en cuenta que el resultado es un objeto terra:
munic.slope## class : SpatRaster
## dimensions : 1455, 1951, 1 (nrow, ncol, nlyr)
## resolution : 80, 80 (x, y)
## extent : 317990.2, 474070.2, 504609.7, 621009.7 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
## source(s) : memory
## name : slope
## min value : 2.207806
## max value : 28.232075Ahora, veamos el mapa
tmap_mode("view")## tmap mode set to interactive viewing tm_shape(munic) + tm_borders(col = "darkgray", lwd = 0.5, alpha=0.8) +
tm_text(text = "NAME_2", size = 0.5, alpha=0.7) +
tm_shape(munic.slope) +
tm_raster(style="cont", n=7, alpha=0.6, title="Mean slope")+
tm_layout(legend.outside = TRUE)## stars object downsampled to 1158 by 863 cells. See tm_shape manual (argument raster.downsample)10. Bibliografía
Si reutiliza el código de este cuaderno, cite este trabajo como:
Lizarazo, I., 2023. Procesamiento y análisis de datos de elevación en R. Disponible en https://rpubs.com/ials2un/dem_analysis
sessionInfo()## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8 LC_CTYPE=Spanish_Colombia.utf8
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Colombia.utf8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] RColorBrewer_1.1-3 tmap_3.3-3 leaflet_2.1.2 mapview_2.11.0
## [5] lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0 dplyr_1.1.0
## [9] purrr_1.0.1 readr_2.1.4 tidyr_1.3.0 tibble_3.1.8
## [13] ggplot2_3.4.1 tidyverse_2.0.0 sf_1.0-12 elevatr_0.4.2
## [17] rgdal_1.6-6 rgeos_0.6-3 terra_1.7-29 raster_3.6-20
## [21] sp_1.6-0 rasterVis_0.51.5 lattice_0.20-45
##
## loaded via a namespace (and not attached):
## [1] satellite_1.0.4 httr_1.4.5 progress_1.2.2
## [4] webshot_0.5.4 tools_4.2.2 bslib_0.4.2
## [7] utf8_1.2.3 R6_2.5.1 KernSmooth_2.23-20
## [10] DBI_1.1.3 colorspace_2.1-0 withr_2.5.0
## [13] prettyunits_1.1.1 tidyselect_1.2.0 curl_5.0.0
## [16] compiler_4.2.2 progressr_0.13.0 leafem_0.2.0
## [19] cli_3.6.0 sass_0.4.5 scales_1.2.1
## [22] classInt_0.4-9 hexbin_1.28.3 proxy_0.4-27
## [25] digest_0.6.31 rmarkdown_2.20 dichromat_2.0-0.1
## [28] base64enc_0.1-3 jpeg_0.1-10 pkgconfig_2.0.3
## [31] htmltools_0.5.4 highr_0.10 fastmap_1.1.0
## [34] htmlwidgets_1.6.2 rlang_1.0.6 rstudioapi_0.14
## [37] farver_2.1.1 jquerylib_0.1.4 generics_0.1.3
## [40] zoo_1.8-12 jsonlite_1.8.4 crosstalk_1.2.0
## [43] magrittr_2.0.3 s2_1.1.3 interp_1.1-4
## [46] Rcpp_1.0.10 munsell_0.5.0 fansi_1.0.4
## [49] abind_1.4-5 lifecycle_1.0.3 stringi_1.7.12
## [52] leafsync_0.1.0 yaml_2.3.7 tmaptools_3.1-1
## [55] grid_4.2.2 parallel_4.2.2 crayon_1.5.2
## [58] deldir_1.0-6 stars_0.6-1 hms_1.1.2
## [61] knitr_1.42 pillar_1.8.1 markdown_1.5
## [64] codetools_0.2-18 stats4_4.2.2 wk_0.7.2
## [67] XML_3.99-0.14 glue_1.6.2 evaluate_0.20
## [70] leaflet.providers_1.9.0 latticeExtra_0.6-30 png_0.1-8
## [73] vctrs_0.5.2 tzdb_0.3.0 gtable_0.3.1
## [76] slippymath_0.3.1 cachem_1.0.6 xfun_0.37
## [79] mime_0.12 lwgeom_0.2-11 e1071_1.7-13
## [82] class_7.3-20 viridisLite_0.4.1 units_0.8-2
## [85] timechange_0.2.0 ellipsis_0.3.2