1. ¿Por qué este cuaderno?

Este cuaderno ilustra varias funcionalidades para obtener, procesar y visualizar modelos digitales de elevación (DEM) en R. Su objetivo es ayudar a los estudiantes de Geomática Básica de la Universidad Nacional a familiarizarse con los DEM y las variables geomorfométricas.

2. Configuraciones

Los paquetes requeridos deben instalarse previamente desde la Consola.

#install.packages("rgdal")
#install.packages("raster")
#install.packages("rgeos")
#install.packages("terra")
#install.packages("elevatr")
#install.packages("rasterVis")
#install.packages("rgl")
#install.packages("mapview")
#install.packages("st")

Ahora se debe limpiar la memoria:

rm(list=ls())

Ahora, carguemos las bibliotecas:

library(rasterVis)
## Warning: package 'rasterVis' was built under R version 4.3.2
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 4.3.2
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, were retired in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
library(raster)
## Loading required package: sp
library(terra)
## terra 1.7.55
library(elevatr)
## Warning: package 'elevatr' was built under R version 4.3.2
## elevatr v0.99.0 NOTE: Version 0.99.0 of 'elevatr' uses 'sf' and 'terra'.  Use 
## of the 'sp', 'raster', and underlying 'rgdal' packages by 'elevatr' is being 
## deprecated; however, get_elev_raster continues to return a RasterLayer.  This 
## will be dropped in future versions, so please plan accordingly.
library(st)
## Warning: package 'st' was built under R version 4.3.2
## Loading required package: sda
## Loading required package: entropy
## Loading required package: corpcor
## Loading required package: fdrtool
## 
## Attaching package: 'sda'
## The following object is masked from 'package:terra':
## 
##     centroids
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2
## ── 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()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(mapview)
## Warning: package 'mapview' was built under R version 4.3.2
library(ggplot2)
library(leaflet)
library(tmap)
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
library(RColorBrewer)
library(sf)
## Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE

3. Introduccion a la elevación

Los datos de elevación se utilizan para una amplia gama de aplicaciones, incluidas, por ejemplo, visualización, hidrología y modelización ecológica. Obtener acceso a estos datos en R no ha tenido una única interfaz, está disponible a través de funciones en muchos paquetes o requiere acceso local a los datos. Esto ya no es necesario porque ahora existe una variedad de API que brindan acceso programático a los datos de elevación. El paquete elevatr se escribió para estandarizar el acceso a los datos de elevación desde las API web.

Para acceder a datos de elevación ráster (por ejemplo, un DEM), el paquete elevatr utiliza Terrain Tiles de Amazon Web Services. Exploraremos esta funcionalidad en este cuaderno.

Existen varias fuentes de modelos de elevación digitales, como Shuttle Radar Topography Mission (SRTM), el conjunto de datos de elevación nacional (NED) del USGS, Global DEM (GDEM) y otros. Cada uno de estos DEM tiene ventajas y desventajas en 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 determinada en un nivel de zoom determinado. Aunque cerrados, estos datos compilados por Mapzen son actualmente accesibles a través de Terrain Tiles en Amazon Web Services (AWS).

La entrada para get_elev_raster() puede ser un data.frame con ubicaciones x (longitud) e y (latitud) como las dos primeras columnas, cualquier objeto sp o cualquier objeto ráster y devuelve un RasterLayer de los mosaicos que se superponen al cuadro delimitador. de la entrada. Si se recuperan varios mosaicos, el resultado resultante es una capa ráster fusionada.

** Usando get_elev_raster() para acceder a Terrain Tiles en AWS**.

Como se mencionó, un marco de datos con columnas xey, un objeto sp o un objeto ráster debe ser la entrada y el src debe configurarse en “mapzen” (este es el valor predeterminado).

No hay diferencia en el uso de los tipos de datos de entrada sp y ráster.

El marco de datos requiere un CRS (es decir, un sistema de referencia de coordenadas). El nivel de zoom (z) está predeterminado en 9 (una compensación entre resolución y tiempo de descarga). Podrías empezar a probar con un nivel de zoom más alto (por ejemplo, 10).

4. Obtenga datos de elevación ráster para su departamento

Ahora, necesitamos cargar un shapefile o un geopaquete que represente a los municipios de nuestro departamento. En este cuaderno utilizaré un shapefile descargado del geoportal del DANE. Leamos el shapefile usando una función proporcionada por el paquete sf:

munic <-  st_read("./datos/santander-mun.shp")
## Reading layer `santander-mun' from data source 
##   `D:\CuadernoAguacate\datos\santander-mun.shp' using driver `ESRI Shapefile'
## Simple feature collection with 87 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.52895 ymin: 5.719626 xmax: -72.5161 ymax: 8.119248
## Geodetic CRS:  MAGNA-SIRGAS

Revisemos la geometría, el cuadro delimitador, el CRS y los atributos del objeto munic

head(munic)
## Simple feature collection with 6 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.72281 ymin: 6.046887 xmax: -73.2271 ymax: 6.384883
## Geodetic CRS:  MAGNA-SIRGAS
##   DPTO_CCDGO MPIO_CCDGO   MPIO_CNMBR MPIO_CDPMP      AREA  LATITUD  LONGITUD
## 1         68        013       AGUADA      68013  75231248 6.182208 -73.53039
## 2         68        673   SAN BENITO      68673  59113022 6.101681 -73.52749
## 3         68        770       SUAITA      68770 285195534 6.114763 -73.36621
## 4         68        397       LA PAZ      68397 284346972 6.217374 -73.62977
## 5         68        245 EL GUACAMAYO      68245  92967630 6.253548 -73.53591
## 6         68        320    GUADALUPE      68320 146330709 6.233278 -73.40977
##   Shape_Leng  Shape_Area                       geometry
## 1  0.4758097 0.006146093 MULTIPOLYGON (((-73.56261 6...
## 2  0.4121582 0.004828582 MULTIPOLYGON (((-73.517 6.0...
## 3  0.9177433 0.023294974 MULTIPOLYGON (((-73.38416 6...
## 4  0.9884018 0.023232222 MULTIPOLYGON (((-73.66446 6...
## 5  0.6013367 0.007596097 MULTIPOLYGON (((-73.5996 6....
## 6  0.5817189 0.011955233 MULTIPOLYGON (((-73.38955 6...

Creemos un nuevo objeto con centroides de municipios, para esto es necesario delimitar el area, confirmar geometria, y usar la funcion centroide de la libreria st_

munic$km2 = munic$AREA / 1000000
munic$km2
##  [1]   75.23125   59.11302  285.19553  284.34697   92.96763  146.33071
##  [7]  166.21697  522.53920  335.20158  251.67855  176.59966   89.59309
## [13]  169.79155  107.96711   39.89918  290.12519  244.22274  484.57231
## [19]   46.66489  447.67673   83.79546   30.15675  137.27581   65.57431
## [25]  146.87727   97.60691  364.41855   19.69444  201.73158 1326.83830
## [31]  431.24871 1402.87001 1116.36061  904.71295  762.34460 1507.98625
## [37]  924.03262  493.02989  550.43324 3174.28854  590.21639 1010.12586
## [43]  398.21775  126.70625 1168.43222  330.72628  235.74857  100.32774
## [49]  152.91569  127.82583   55.53542  169.01757   45.01346  363.00877
## [55]   85.28406   75.30374   71.21135  489.22844   74.81603  258.94615
## [61]  337.11992  102.91751  281.55886  483.97013  180.47170  302.39826
## [67]  421.08144  419.23135   71.42170   75.22128  286.06675   72.98112
## [73]  577.49575  400.43269  608.23307   92.42239  173.21902  373.85726
## [79]  137.17208   67.30036   94.93682   57.18170   58.54112  454.98192
## [85]  221.21709  143.13343   81.11494
munic$geometry
## Geometry set for 87 features 
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.52895 ymin: 5.719626 xmax: -72.5161 ymax: 8.119248
## Geodetic CRS:  MAGNA-SIRGAS
## First 5 geometries:
## MULTIPOLYGON (((-73.56261 6.240321, -73.56222 6...
## MULTIPOLYGON (((-73.517 6.090301, -73.5256 6.08...
## MULTIPOLYGON (((-73.38416 6.194359, -73.38407 6...
## MULTIPOLYGON (((-73.66446 6.380951, -73.66419 6...
## MULTIPOLYGON (((-73.5996 6.318598, -73.5997 6.3...
(centers <- st_centroid(munic))
## Warning: st_centroid assumes attributes are constant over geometries
## Simple feature collection with 87 features and 10 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.16858 ymin: 5.783723 xmax: -72.58294 ymax: 7.539143
## Geodetic CRS:  MAGNA-SIRGAS
## First 10 features:
##    DPTO_CCDGO MPIO_CCDGO      MPIO_CNMBR MPIO_CDPMP      AREA  LATITUD
## 1          68        013          AGUADA      68013  75231248 6.182208
## 2          68        673      SAN BENITO      68673  59113022 6.101681
## 3          68        770          SUAITA      68770 285195534 6.114763
## 4          68        397          LA PAZ      68397 284346972 6.217374
## 5          68        245    EL GUACAMAYO      68245  92967630 6.253548
## 6          68        320       GUADALUPE      68320 146330709 6.233278
## 7          68        020         ALBANIA      68020 166216966 5.778073
## 8          68        773           SUCRE      68773 522539201 5.983204
## 9          68        377      LA BELLEZA      68377 335201582 5.910233
## 10         68        572 PUENTE NACIONAL      68572 251678548 5.831130
##     LONGITUD Shape_Leng  Shape_Area                   geometry       km2
## 1  -73.53039  0.4758097 0.006146093 POINT (-73.53038 6.182207)  75.23125
## 2  -73.52749  0.4121582 0.004828582 POINT (-73.52497 6.110971)  59.11302
## 3  -73.36621  0.9177433 0.023294974 POINT (-73.36621 6.114765) 285.19553
## 4  -73.62977  0.9884018 0.023232222 POINT (-73.62976 6.217368) 284.34697
## 5  -73.53591  0.6013367 0.007596097 POINT (-73.53591 6.253549)  92.96763
## 6  -73.40977  0.5817189 0.011955233 POINT (-73.40978 6.233277) 146.33071
## 7  -73.82846  0.8761299 0.013570466 POINT (-73.82851 5.783723) 166.21697
## 8  -73.95936  1.7366067 0.042677592 POINT (-73.95935 5.983204) 522.53920
## 9  -74.05072  1.1748849 0.027373676 POINT (-74.05072 5.910234) 335.20158
## 10 -73.67842  0.7559087 0.020549103 POINT (-73.68648 5.832654) 251.67855
centers$x = st_coordinates(centers)[,1]
centers$x
##  [1] -73.53038 -73.52497 -73.36621 -73.62976 -73.53591 -73.40978 -73.82851
##  [8] -73.95935 -74.05072 -73.68648 -73.95524 -73.82342 -73.01163 -72.92925
## [15] -73.09961 -73.10785 -73.01873 -73.00517 -73.63474 -73.70244 -73.72174
## [22] -73.58450 -73.21518 -73.25035 -73.11888 -73.16397 -73.27975 -73.28901
## [29] -73.32467 -73.77992 -73.38741 -73.57128 -73.53900 -73.69002 -73.95760
## [36] -73.78818 -73.56475 -73.24870 -73.29137 -74.16858 -73.81869 -74.11162
## [43] -73.92718 -72.99034 -73.33520 -72.95671 -73.05439 -73.06799 -73.11157
## [50] -73.24436 -73.17571 -73.36179 -72.93920 -73.00468 -72.67435 -72.73976
## [57] -72.64322 -72.83067 -72.68028 -72.58294 -72.63002 -72.59819 -72.81614
## [64] -72.95836 -72.81926 -72.83326 -72.65428 -73.15402 -73.23985 -73.12212
## [71] -73.27962 -73.17148 -72.98247 -73.06463 -73.33341 -72.89896 -73.41476
## [78] -73.58207 -73.50211 -73.33115 -73.63239 -72.73841 -73.28273 -73.18258
## [85] -72.91047 -72.85638 -73.10997
centers$x = st_coordinates(centers)[,1]
centers$x
##  [1] -73.53038 -73.52497 -73.36621 -73.62976 -73.53591 -73.40978 -73.82851
##  [8] -73.95935 -74.05072 -73.68648 -73.95524 -73.82342 -73.01163 -72.92925
## [15] -73.09961 -73.10785 -73.01873 -73.00517 -73.63474 -73.70244 -73.72174
## [22] -73.58450 -73.21518 -73.25035 -73.11888 -73.16397 -73.27975 -73.28901
## [29] -73.32467 -73.77992 -73.38741 -73.57128 -73.53900 -73.69002 -73.95760
## [36] -73.78818 -73.56475 -73.24870 -73.29137 -74.16858 -73.81869 -74.11162
## [43] -73.92718 -72.99034 -73.33520 -72.95671 -73.05439 -73.06799 -73.11157
## [50] -73.24436 -73.17571 -73.36179 -72.93920 -73.00468 -72.67435 -72.73976
## [57] -72.64322 -72.83067 -72.68028 -72.58294 -72.63002 -72.59819 -72.81614
## [64] -72.95836 -72.81926 -72.83326 -72.65428 -73.15402 -73.23985 -73.12212
## [71] -73.27962 -73.17148 -72.98247 -73.06463 -73.33341 -72.89896 -73.41476
## [78] -73.58207 -73.50211 -73.33115 -73.63239 -72.73841 -73.28273 -73.18258
## [85] -72.91047 -72.85638 -73.10997
centers$y = st_coordinates(centers)[,2]
centers$y
##  [1] 6.182207 6.110971 6.114765 6.217368 6.253549 6.233277 5.783723 5.983204
##  [9] 5.910234 5.832654 5.801121 5.854289 6.716764 6.752464 6.713780 6.812161
## [17] 6.618346 6.962446 5.954781 6.220505 5.944894 6.033334 6.647059 6.562542
## [25] 6.550777 6.685865 6.834620 6.524753 6.671808 7.055513 7.013201 7.406521
## [33] 6.894287 6.648092 6.685271 7.539143 6.664257 7.036680 7.204050 6.407867
## [41] 6.331489 6.096642 6.096729 7.259038 7.438398 7.168023 7.348772 7.079704
## [49] 7.155830 6.462088 6.513892 6.559878 7.341422 7.458632 6.520198 6.629398
## [57] 6.568121 6.320224 6.658240 6.657662 6.776001 6.505505 6.796478 6.492526
## [65] 6.644299 6.920571 6.891426 6.195528 6.355620 6.355368 6.225593 6.416628
## [73] 6.235449 6.081349 5.942388 7.314570 6.372676 6.408904 6.303117 6.319299
## [81] 6.067944 6.720222 6.391750 7.517433 6.955722 6.456285 6.434578
munic$Long <- centers$x 
munic$Lat <- centers$y

Con la geometria establecida vamos a crear el mapa

library(leaflet)

# Creacion del mapa
map <- leaflet(munic) %>%
  addTiles() %>% 
  
  # Agregar un rectángulo
  addRectangles(
    lng1 = 73.90, lat1 = 4.9,
    lng2 = 72.7, lat2 = 5.9,
    fillColor = "transparent"
  ) %>%
  
  # Agregar polígonos
  addPolygons(
    color = "#444444", weight = 1, smoothFactor = 0.5,
    opacity = 1.0, fillOpacity = 0.4,
    fillColor = colorQuantile("YlOrRd", munic$km2)(munic$km2),
    highlightOptions = highlightOptions(
      color = "white", weight = 2,
      bringToFront = TRUE
    )
  ) %>%
  
  # Agregar marcadores (puntos)
  addMarkers(
    lng = munic$Long, lat = munic$Lat,
    popup = ~munic$NAME_2
  )
## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs).
## Need '+proj=longlat +datum=WGS84'
map
library(sf)

# Convertir las características simples (munic) a objetos sf
sp.munic <- st_as_sf(munic)

# Encontrar los centroides de las regiones
centers <- st_centroid(sp.munic)
## Warning: st_centroid assumes attributes are constant over geometries
# Agregar el nombre de la región y el nombre del municipio a los centroides
centers$region <- row.names(sp.munic)
centers$munic <- munic$MPIO_CNMBR
centers
## Simple feature collection with 87 features and 14 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.16858 ymin: 5.783723 xmax: -72.58294 ymax: 7.539143
## Geodetic CRS:  MAGNA-SIRGAS
## First 10 features:
##    DPTO_CCDGO MPIO_CCDGO      MPIO_CNMBR MPIO_CDPMP      AREA  LATITUD
## 1          68        013          AGUADA      68013  75231248 6.182208
## 2          68        673      SAN BENITO      68673  59113022 6.101681
## 3          68        770          SUAITA      68770 285195534 6.114763
## 4          68        397          LA PAZ      68397 284346972 6.217374
## 5          68        245    EL GUACAMAYO      68245  92967630 6.253548
## 6          68        320       GUADALUPE      68320 146330709 6.233278
## 7          68        020         ALBANIA      68020 166216966 5.778073
## 8          68        773           SUCRE      68773 522539201 5.983204
## 9          68        377      LA BELLEZA      68377 335201582 5.910233
## 10         68        572 PUENTE NACIONAL      68572 251678548 5.831130
##     LONGITUD Shape_Leng  Shape_Area                   geometry       km2
## 1  -73.53039  0.4758097 0.006146093 POINT (-73.53038 6.182207)  75.23125
## 2  -73.52749  0.4121582 0.004828582 POINT (-73.52497 6.110971)  59.11302
## 3  -73.36621  0.9177433 0.023294974 POINT (-73.36621 6.114765) 285.19553
## 4  -73.62977  0.9884018 0.023232222 POINT (-73.62976 6.217368) 284.34697
## 5  -73.53591  0.6013367 0.007596097 POINT (-73.53591 6.253549)  92.96763
## 6  -73.40977  0.5817189 0.011955233 POINT (-73.40978 6.233277) 146.33071
## 7  -73.82846  0.8761299 0.013570466 POINT (-73.82851 5.783723) 166.21697
## 8  -73.95936  1.7366067 0.042677592 POINT (-73.95935 5.983204) 522.53920
## 9  -74.05072  1.1748849 0.027373676 POINT (-74.05072 5.910234) 335.20158
## 10 -73.67842  0.7559087 0.020549103 POINT (-73.68648 5.832654) 251.67855
##         Long      Lat region           munic
## 1  -73.53038 6.182207      1          AGUADA
## 2  -73.52497 6.110971      2      SAN BENITO
## 3  -73.36621 6.114765      3          SUAITA
## 4  -73.62976 6.217368      4          LA PAZ
## 5  -73.53591 6.253549      5    EL GUACAMAYO
## 6  -73.40978 6.233277      6       GUADALUPE
## 7  -73.82851 5.783723      7         ALBANIA
## 8  -73.95935 5.983204      8           SUCRE
## 9  -74.05072 5.910234      9      LA BELLEZA
## 10 -73.68648 5.832654     10 PUENTE NACIONAL
# z denota el nivel de zoom de los datos
# cuanto mayor sea el zoom, mayor será la resolución espacial
elevation <- get_elev_raster(munic, z = 8)
## Mosaicing & Projecting
## Note: Elevation units are in meters.

Ahora miremos la elevacion

elevation
## class      : RasterLayer 
## dimensions : 1020, 1028, 1048560  (nrow, ncol, ncell)
## resolution : 0.002736193, 0.002736193  (x, y)
## extent     : -74.53125, -71.71844, 5.616251, 8.407168  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
## source     : file229c365244a8.tif 
## names      : file229c365244a8
# Crear un objeto sf a partir de munic
sf_munic <- st_as_sf(munic, coords = c("Long", "Lat"))

# Obtener los centroides de las regiones
centers <- st_centroid(sf_munic)
## Warning: st_centroid assumes attributes are constant over geometries
# Agregar el nombre de la región y el nombre del municipio a los centroides
centers$region <- row.names(sf_munic)
centers$munic <- munic$MPIO_CNMBR

# Imprimir el resultado
print(centers)
## Simple feature collection with 87 features and 14 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.16858 ymin: 5.783723 xmax: -72.58294 ymax: 7.539143
## Geodetic CRS:  MAGNA-SIRGAS
## First 10 features:
##    DPTO_CCDGO MPIO_CCDGO      MPIO_CNMBR MPIO_CDPMP      AREA  LATITUD
## 1          68        013          AGUADA      68013  75231248 6.182208
## 2          68        673      SAN BENITO      68673  59113022 6.101681
## 3          68        770          SUAITA      68770 285195534 6.114763
## 4          68        397          LA PAZ      68397 284346972 6.217374
## 5          68        245    EL GUACAMAYO      68245  92967630 6.253548
## 6          68        320       GUADALUPE      68320 146330709 6.233278
## 7          68        020         ALBANIA      68020 166216966 5.778073
## 8          68        773           SUCRE      68773 522539201 5.983204
## 9          68        377      LA BELLEZA      68377 335201582 5.910233
## 10         68        572 PUENTE NACIONAL      68572 251678548 5.831130
##     LONGITUD Shape_Leng  Shape_Area                   geometry       km2
## 1  -73.53039  0.4758097 0.006146093 POINT (-73.53038 6.182207)  75.23125
## 2  -73.52749  0.4121582 0.004828582 POINT (-73.52497 6.110971)  59.11302
## 3  -73.36621  0.9177433 0.023294974 POINT (-73.36621 6.114765) 285.19553
## 4  -73.62977  0.9884018 0.023232222 POINT (-73.62976 6.217368) 284.34697
## 5  -73.53591  0.6013367 0.007596097 POINT (-73.53591 6.253549)  92.96763
## 6  -73.40977  0.5817189 0.011955233 POINT (-73.40978 6.233277) 146.33071
## 7  -73.82846  0.8761299 0.013570466 POINT (-73.82851 5.783723) 166.21697
## 8  -73.95936  1.7366067 0.042677592 POINT (-73.95935 5.983204) 522.53920
## 9  -74.05072  1.1748849 0.027373676 POINT (-74.05072 5.910234) 335.20158
## 10 -73.67842  0.7559087 0.020549103 POINT (-73.68648 5.832654) 251.67855
##         Long      Lat region           munic
## 1  -73.53038 6.182207      1          AGUADA
## 2  -73.52497 6.110971      2      SAN BENITO
## 3  -73.36621 6.114765      3          SUAITA
## 4  -73.62976 6.217368      4          LA PAZ
## 5  -73.53591 6.253549      5    EL GUACAMAYO
## 6  -73.40978 6.233277      6       GUADALUPE
## 7  -73.82851 5.783723      7         ALBANIA
## 8  -73.95935 5.983204      8           SUCRE
## 9  -74.05072 5.910234      9      LA BELLEZA
## 10 -73.68648 5.832654     10 PUENTE NACIONAL

Observe algunas cosas sobre este DEM:

dimensiones: el “tamaño” del archivo en píxeles nrow, ncol: el número de filas y columnas de los datos (imagine una hoja de cálculo o una matriz). ncells: el número total de píxeles o celdas que componen el ráster. resolución: el tamaño de cada píxel (en grados decimales en este caso). Recordemos que 1 grado decimal representa unos 111,11 km en el Ecuador. extensión: la extensión espacial del ráster. Este valor estará en las mismas unidades de coordenadas que el sistema de referencia de coordenadas del ráster. crs: la cadena del sistema de referencia de coordenadas para el ráster. Este ráster está en coordenadas geográficas con un datum de WGS 84. En su cuaderno, escriba el tamaño de celda de su ráster (es decir, la resolución espacial), en metros.

Tenga en cuenta que la elevación es una capa ráster temporal (es decir, sólo existe en la memoria). Sin embargo, es posible escribir el DEM en su directorio local usando la función writeRaster proporcionada por la biblioteca rgdal:

#5 Visualizar los datos de elevacion Para visualizar es mejor usar una resolución más baja para acelera el proceso:

# Mosaico de los datos de elevación a una resolución más baja (opcional)
elevation2 <- aggregate(elevation, fact = 4, fun = mean)
sf_munic <- st_transform(sf_munic, crs = 4326)
pal <- colorNumeric(c("forestgreen","yellow","tan","brown"), values(elevation2),
  na.color = "transparent")
sf_munic <- st_transform(sf_munic, crs = 4326)
centers <- st_transform(centers, crs = 4326)

usaremos la biblioteca de folletos para ver los datos de elevación:

leaflet(munic) %>% 
  addTiles() %>%
  addPolygons(
    color = "white", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.25, fillOpacity = 0.15,
    popup = ~paste("Region: ", MPIO_CNMBR, "<br>Km2: ", km2, "<br>")
  ) %>%
  addLabelOnlyMarkers(
    data = centers,
    lng = ~LONGITUD, lat = ~LATITUD, 
  # Ajusta las variables 'LONGITUD' y 'LATITUD' según tus datos
    label = ~region,
    labelOptions = labelOptions(noHide = TRUE, direction = 'top', textOnly = TRUE)
  ) %>%
  addRasterImage(
    elevation2, colors = pal, opacity = 0.9
  ) %>%
  addLegend(
    pal = pal, values = values(elevation2),
    title = "Datos de elevación Santander(mts)"
  )
## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs).
## Need '+proj=longlat +datum=WGS84'

Se puede hacer clic en el ID de cada municipio para obtener el nombre y la extensión en Km^2

#6 Reproyectar los datos de elevación

Cuando se trabaja con DEM, siempre es una buena idea utilizar coordenadas de mapa en lugar de coordenadas geográficas. Esto se debe al hecho de que, en las coordenadas geográficas, las unidades de dimensión horizontal son grados decimales, PERO la unidad de dimensión vertical son metros.

Para reproducir los datos de elevación, realizamos dos pasos.

Primero, definimos el CRS objetivo:

# WGS 84 / UTM zone 18N
spatialref <- CRS(SRS_string="EPSG:32618")

Ahora reproyectamos el DEM

pr3 <- projectExtent(elevation, crs=spatialref)
res(pr3) <- 80
rep_elev <- projectRaster(elevation, pr3)

Miremos

rep_elev
## class      : RasterLayer 
## dimensions : 3875, 3900, 15112500  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 551603.9, 863603.9, 620826, 930826  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : r_tmp_2023-11-09_214838.452362_8860_22098.grd 
## names      : layer 
## values     : -416.3299, 5233.17  (min, max)
rep_munic <- st_transform(sp.munic, crs = spatialref)
rep_munic
## Simple feature collection with 87 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 552105.4 ymin: 632385.5 xmax: 774559.6 ymax: 897672.5
## Projected CRS: WGS 84 / UTM zone 18N
## First 10 features:
##    DPTO_CCDGO MPIO_CCDGO      MPIO_CNMBR MPIO_CDPMP      AREA  LATITUD
## 1          68        013          AGUADA      68013  75231248 6.182208
## 2          68        673      SAN BENITO      68673  59113022 6.101681
## 3          68        770          SUAITA      68770 285195534 6.114763
## 4          68        397          LA PAZ      68397 284346972 6.217374
## 5          68        245    EL GUACAMAYO      68245  92967630 6.253548
## 6          68        320       GUADALUPE      68320 146330709 6.233278
## 7          68        020         ALBANIA      68020 166216966 5.778073
## 8          68        773           SUCRE      68773 522539201 5.983204
## 9          68        377      LA BELLEZA      68377 335201582 5.910233
## 10         68        572 PUENTE NACIONAL      68572 251678548 5.831130
##     LONGITUD Shape_Leng  Shape_Area                       geometry       km2
## 1  -73.53039  0.4758097 0.006146093 MULTIPOLYGON (((659020.5 68...  75.23125
## 2  -73.52749  0.4121582 0.004828582 MULTIPOLYGON (((664113.8 67...  59.11302
## 3  -73.36621  0.9177433 0.023294974 MULTIPOLYGON (((678782.9 68... 285.19553
## 4  -73.62977  0.9884018 0.023232222 MULTIPOLYGON (((647710.4 70... 284.34697
## 5  -73.53591  0.6013367 0.007596097 MULTIPOLYGON (((654904.7 69...  92.96763
## 6  -73.40977  0.5817189 0.011955233 MULTIPOLYGON (((678146.1 69... 146.33071
## 7  -73.82846  0.8761299 0.013570466 MULTIPOLYGON (((635572.1 63... 166.21697
## 8  -73.95936  1.7366067 0.042677592 MULTIPOLYGON (((599979.8 68... 522.53920
## 9  -74.05072  1.1748849 0.027373676 MULTIPOLYGON (((596996.8 66... 335.20158
## 10 -73.67842  0.7559087 0.020549103 MULTIPOLYGON (((651103.5 65... 251.67855
##         Long      Lat
## 1  -73.53038 6.182207
## 2  -73.52497 6.110971
## 3  -73.36621 6.114765
## 4  -73.62976 6.217368
## 5  -73.53591 6.253549
## 6  -73.40978 6.233277
## 7  -73.82851 5.783723
## 8  -73.95935 5.983204
## 9  -74.05072 5.910234
## 10 -73.68648 5.832654
# Ruta de destino para guardar el DEM reproyectado
filename <- "./rep_Santander_dem.tif"

# Especifica el tipo de dato (datatype) según tus necesidades
datatype <- 'INT4S'  # Ajusta esto según el tipo de datos que desees

# Escribe el DEM reproyectado en el archivo
writeRaster(rep_elev, filename = filename, datatype = datatype, overwrite = TRUE)

#7 Estadísticas básicas de datos de elevación.

Una exploración rápida de las estadísticas de DEM.

Este fragmento “limpia” la elevación inferior a 0,0:

rep_elev[rep_elev < 0.0] <- 0.0

creemos el histograma de valores de elevación:

hist(rep_elev)

Estadisticas DEM

promedio <- cellStats(rep_elev, 'mean')
minimo <- cellStats(rep_elev, 'min')
maximo <- cellStats(rep_elev, 'max')
desviacion  <- cellStats(rep_elev, 'sd')
metricas <- c('mean', 'min', 'max', 'std')
valores <- c(promedio, minimo, maximo, desviacion)

Resultado

(df_estadisticas <- data.frame(metricas, valores))
##   metricas  valores
## 1     mean 1240.658
## 2      min    0.000
## 3      max 5233.170
## 4      std 1102.644

#8 Obtención de variables geomorfométricas.

Antes de continuar, asegúrese de haber comprendido los conceptos básicos sobre: DEM, pendiente y orientación.

Primero, calcule la pendiente, la orientación 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)
slope
## class      : RasterLayer 
## dimensions : 3875, 3900, 15112500  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 551603.9, 863603.9, 620826, 930826  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : r_tmp_2023-11-09_214857.408576_8860_01861.grd 
## names      : slope 
## values     : 0, 70.63123  (min, max)
aspect
## class      : RasterLayer 
## dimensions : 3875, 3900, 15112500  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 551603.9, 863603.9, 620826, 930826  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : aspect 
## values     : 0, 360  (min, max)

Para trazar la pendiente, es una buena idea reducir el tamaño de RasterLayer:

slope2 <- aggregate(slope, fact=4, fun=mean)
slope3 <- projectRasterForLeaflet(slope2, "ngb")

Paleta de colores:

pal2 <- colorNumeric(c("#00FF00", "#FF0000"),   values(slope3),na.color = "transparent")

Esta paleta va desde verde (#00FF00), que puede representar pendientes más suaves, hasta rojo (#FF0000), que puede representar pendientes más empinadas

Ahora el codigo para el grafico:

leaflet(munic) %>% addTiles() %>%
    addPolygons(color = "white", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.25, fillOpacity = 0.15,
    popup = paste("Region: ", munic$MPIO_CNMBR, "<br>",
                          "Km2: ", munic$km2, "<br>")) %>%
  addLabelOnlyMarkers(data = centers,
                    lng = ~LONGITUD, lat = ~LATITUD, label = ~region,
                    labelOptions = labelOptions(noHide = TRUE, direction = 'top', textOnly = TRUE)) %>%
  addRasterImage(slope3, colors = pal2, opacity = 1.0)  %>%
  addLegend(pal = pal2, values = values(slope3),
    title = "Pendiente de Santander (%)")
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs).
## Need '+proj=longlat +datum=WGS84'
slope2
## class      : RasterLayer 
## dimensions : 969, 975, 944775  (nrow, ncol, ncell)
## resolution : 320, 320  (x, y)
## extent     : 551603.9, 863603.9, 620746, 930826  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : slope 
## values     : 0, 66.74435  (min, max)
(nslope2  <- as(slope2, "SpatRaster"))
## class       : SpatRaster 
## dimensions  : 969, 975, 1  (nrow, ncol, nlyr)
## resolution  : 320, 320  (x, y)
## extent      : 551603.9, 863603.9, 620746, 930826  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        :    slope 
## min value   :  0.00000 
## max value   : 66.74435

La clasificación de pendientes es una forma común de entender la variabilidad del relieve. Es necesario incluir aquí una tabla con intervalos de pendientes, nombres de clases y descripción, según lo define el Instituto Geográfico Colombiano (IGAC).

Realicemos una clasificación de pendientes basada en rangos comunes. Primero, creamos 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)
(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
slope_rec <- terra::classify(nslope2, rclmat, right=TRUE)
slope_rec
## class       : SpatRaster 
## dimensions  : 969, 975, 1  (nrow, ncol, nlyr)
## resolution  : 320, 320  (x, y)
## extent      : 551603.9, 863603.9, 620746, 930826  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        : slope 
## min value   :     0 
## max value   :     6

#9 Visualización estática.

Usaremos tmap para crear y guardar un mapa de clasificación de pendientes estático.

intervals <- c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130)

# Definir colores para los intervalos de clase
colors <- c("#f7f7f7", "#d9f0a3", "#addd8e", "#78c679", "#41ab5d",
            "#238443", "#006837", "#005a32", "#00441b", "#00341a",
            "#002d19", "#002616", "#002014")

# Crear un mapa de pendiente con intervalos de clase y colores personalizados
slope.map <- tm_shape(slope_rec) +
  tm_raster(style = "cat", title = "Pendiente", breaks = intervals, 
            palette = colors, legend.show = TRUE) +  
  tm_shape(munic) + tm_borders(col = "darkgray", lwd = 0.5, alpha = 0.8) +
  tm_text(text = "MPIO_CNMBR", size = 0.3, 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.01, legend.position = c("left", "top"), 
            title = 'Clases de pendiente', title.position = c('right', 'top')) +
  tm_credits("Data source: Mapzen & DANE", size = 0.3)

el resultado

slope.map

tmap_save(slope.map, "./Santander_slope.png", width = 8.5, height =11, dpi=600)
## Map saved to D:\CuadernoAguacate\Santander_slope.png
## Resolution: 5100 by 6600 pixels
## Size: 8.5 by 11 inches (600 dpi)

#10 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  : 87, 12  (geometries, attributes)
##  extent      : 552105.4, 774559.6, 632385.5, 897672.5  (xmin, xmax, ymin, ymax)
##  coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
##  names       : DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR MPIO_CDPMP      AREA LATITUD
##  type        :      <chr>      <chr>      <chr>      <chr>     <num>   <num>
##  values      :         68        013     AGUADA      68013 7.523e+07   6.182
##                        68        673 SAN BENITO      68673 5.911e+07   6.102
##                        68        770     SUAITA      68770 2.852e+08   6.115
##  LONGITUD Shape_Leng Shape_Area   km2   Long   Lat
##     <num>      <num>      <num> <num>  <num> <num>
##    -73.53     0.4758   0.006146 75.23 -73.53 6.182
##    -73.53     0.4122   0.004829 59.11 -73.52 6.111
##    -73.37     0.9177    0.02329 285.2 -73.37 6.115
(nslope  <- as(slope, "SpatRaster"))
## class       : SpatRaster 
## dimensions  : 3875, 3900, 1  (nrow, ncol, nlyr)
## resolution  : 80, 80  (x, y)
## extent      : 551603.9, 863603.9, 620826, 930826  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source      : r_tmp_2023-11-09_214857.408576_8860_01861.grd 
## name        :    slope 
## min value   :  0.00000 
## max value   : 70.63123
munic.slope <- zonal(nslope, nmunic, fun=mean, as.raster=TRUE)
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 = "MPIO_CNMBR", size = 0.5, alpha=0.7) +
  tm_shape(munic.slope) +
  tm_raster(style="cont", n=7, alpha=0.6, title="Pendiente media")+
tm_layout(legend.outside = TRUE)
## stars object downsampled to 1003 by 997 cells. See tm_shape manual (argument raster.downsample)

#11 bibliografia

sessionInfo() 
## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## 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    
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] sf_1.0-14          RColorBrewer_1.1-3 tmap_3.3-4         leaflet_2.2.0     
##  [5] mapview_2.11.2     lubridate_1.9.2    forcats_1.0.0      stringr_1.5.0     
##  [9] dplyr_1.1.3        purrr_1.0.2        readr_2.1.4        tidyr_1.3.0       
## [13] tibble_3.2.1       ggplot2_3.4.3      tidyverse_2.0.0    st_1.2.7          
## [17] sda_1.3.8          fdrtool_1.2.17     corpcor_1.6.10     entropy_1.3.1     
## [21] elevatr_0.99.0     terra_1.7-55       raster_3.6-26      sp_2.1-0          
## [25] rasterVis_0.51.6   lattice_0.22-5    
## 
## loaded via a namespace (and not attached):
##  [1] DBI_1.1.3                deldir_1.0-9             tmaptools_3.1-1         
##  [4] s2_1.1.4                 rlang_1.1.1              magrittr_2.0.3          
##  [7] e1071_1.7-13             compiler_4.3.1           png_0.1-8               
## [10] vctrs_0.6.3              pkgconfig_2.0.3          wk_0.8.0                
## [13] crayon_1.5.2             fastmap_1.1.1            ellipsis_0.3.2          
## [16] lwgeom_0.2-13            leafem_0.2.3             utf8_1.2.3              
## [19] rmarkdown_2.24           tzdb_0.4.0               xfun_0.40               
## [22] satellite_1.0.4          cachem_1.0.8             jsonlite_1.8.7          
## [25] progress_1.2.2           jpeg_0.1-10              parallel_4.3.1          
## [28] prettyunits_1.1.1        R6_2.5.1                 bslib_0.5.1             
## [31] stringi_1.7.12           jquerylib_0.1.4          stars_0.6-4             
## [34] Rcpp_1.0.11              knitr_1.44               zoo_1.8-12              
## [37] base64enc_0.1-3          leaflet.providers_1.13.0 timechange_0.2.0        
## [40] tidyselect_1.2.0         rstudioapi_0.15.0        dichromat_2.0-0.1       
## [43] abind_1.4-5              yaml_2.3.7               codetools_0.2-19        
## [46] curl_5.0.2               leafsync_0.1.0           withr_2.5.0             
## [49] evaluate_0.21            units_0.8-3              proxy_0.4-27            
## [52] pillar_1.9.0             KernSmooth_2.23-21       stats4_4.3.1            
## [55] generics_0.1.3           hms_1.1.3                munsell_0.5.0           
## [58] scales_1.2.1             class_7.3-22             glue_1.6.2              
## [61] tools_4.3.1              interp_1.1-4             hexbin_1.28.3           
## [64] XML_3.99-0.14            grid_4.3.1               slippymath_0.3.1        
## [67] crosstalk_1.2.0          latticeExtra_0.6-30      colorspace_2.1-0        
## [70] cli_3.6.1                fansi_1.0.4              viridisLite_0.4.2       
## [73] gtable_0.3.4             sass_0.4.7               digest_0.6.33           
## [76] progressr_0.14.0         classInt_0.4-10          htmlwidgets_1.6.2       
## [79] farver_2.1.1             htmltools_0.5.6          lifecycle_1.0.3         
## [82] httr_1.4.7               mime_0.12