PROCESAMIENTO Y ANÁLISIS DE DATOS DE ELEVACIÓN EN R

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 nosotros, estudiantes de Geomática Básica de la Universidad Nacional, a familiarizarse con los DEM y las variables geomorfométricas.

Algunos consejos para escribir tu propio cuaderno:

Configuración

Los paquetes requeridos deben instalarse previamente desde la Consola.

## CONFIGURACIÓN
# INSTALE LOS PAQUETES DESDE LA CONSOLA, NO DESDE ACÁ.
# paquetes = c("rgdal","raster", "rgeos", "terra","elevatr","rasterVis", "rgl", "mapview")
# install.packages(paquetes)

library (raster)

## Loading required package: sp

## 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 (terra)

## terra 1.7.55

library (elevatr)

## 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(rasterVis)

## Loading required package: lattice

library(rgl)
library(mapview)
library(sf)

## Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:terra':
## 
##     intersect, union

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(readr)

Es recomendable empezar a limpiar la memoria:

rm(list=ls())

Ahora, carguemos las bibliotecas:

Al ejecutar un fragmento, preste atención a cualquier mensaje que indique error o advierta 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, necesitamos llamar a dichas funciones usando el modo largo, es decir, nombre_paquete::función.

Después de ejecutar un fragmento por primera vez, es una buena práctica ocultar mensajes y advertencias usando: {r mensaje=FALSE, advertencia=FALSE}.

2. Introducción

Los datos de elevación se utilizan para una amplia gama de aplicaciones, 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 x e y, 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,12,15,etc).

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

Primero, 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 archivo de forma usando una función proporcionada por el paquete sf:

list.files()

##  [1] "ciudadesco.csv"                 "cuaderno_eva"                  
##  [3] "cuaderno4 elevación.rmd"        "CuadernoEVA.nb.html"           
##  [5] "CuadernoEVA.Rmd"                "Cuadernoeva_cundi.html"        
##  [7] "Cuadernoeva_cundi.Rmd"          "CUN_tuberculos_2022.csv"       
##  [9] "Cundimarca_papa.csv"            "Cundinamarca_Papa.csv"         
## [11] "CUNDINARMCA CUADERNO 2.nb.html" "ELEVACIÓN.html"                
## [13] "ELEVACIÓN.nb.html"              "ELEVACIÓN.Rmd"                 
## [15] "EVA_19_22.xlsx"                 "EVA_CUNDI.csv"                 
## [17] "EVA_CUNDI.csv.csv"              "huila_slope.png"               
## [19] "Imagen1.png"                    "imagen2.png"                   
## [21] "imagen3.png"                    "imagen4.png"                   
## [23] "MAPAS-TOLIMA.html"              "MAPAS TOLIMA.Rmd"              
## [25] "MGN_MPIO_POLITICO.cpg"          "MGN_MPIO_POLITICO.dbf"         
## [27] "MGN_MPIO_POLITICO.prj"          "MGN_MPIO_POLITICO.sbn"         
## [29] "MGN_MPIO_POLITICO.sbx"          "MGN_MPIO_POLITICO.shp"         
## [31] "MGN_MPIO_POLITICO.shp.xml"      "MGN_MPIO_POLITICO.shx"         
## [33] "MGN2021_MPIO_POLITICO.rar"      "PRODUCCIÓN-PAPA.html"          
## [35] "PRODUCCIÓN PAPA.Rmd"            "rep_tolima_dem.tif"            
## [37] "rsconnect"                      "Tolima_Algodon.csv"            
## [39] "Tolima_dep.gpkg"                "tolima_mun.gpkg"               
## [41] "tolima_slope.png"

# SpatialPolygonsDataFrame example
munic <-  st_read('MGN_MPIO_POLITICO.shp')

## Reading layer `MGN_MPIO_POLITICO' from data source 
##   `C:\Users\camil\OneDrive\Escritorio\Geomatica\CuadernoEVA\MGN_MPIO_POLITICO.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1121 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -81.73562 ymin: -4.229406 xmax: -66.84722 ymax: 13.39473
## Geodetic CRS:  MAGNA-SIRGAS

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

head(munic)

-Filtremos el departamento de Tolima

munic %>% filter(DPTO_CCDGO=="73")%>% select(MPIO_CNMBR,MPIO_CDPMP, MPIO_NAREA)->
munic1

head(munic1)

Creemos un nuevo objeto con centroides de municipios, pero primero configuramos los datos:

munic1$km2 = munic1$MPIO_NAREA/ 1000000

munic1$km2

##  [1] 1.377251e-03 5.015113e-04 3.443356e-04 2.380448e-04 4.697113e-04
##  [6] 4.406308e-04 1.017631e-03 5.084909e-04 1.908361e-04 1.749306e-04
## [11] 2.102063e-03 3.435981e-04 6.690557e-04 5.075165e-04 6.519504e-04
## [16] 2.164989e-04 1.813222e-04 9.679717e-05 2.196615e-04 5.057438e-04
## [21] 3.227434e-04 2.151819e-04 2.719545e-04 2.823726e-04 2.933432e-04
## [26] 2.022812e-04 4.239322e-04 8.598300e-04 9.538043e-04 6.501113e-05
## [31] 3.544253e-04 1.754130e-03 4.230130e-04 4.022036e-04 2.046219e-03
## [36] 7.758495e-04 7.394762e-04 1.932268e-04 3.834070e-04 4.122262e-04
## [41] 2.701866e-04 1.900782e-04 1.978464e-04 3.346328e-04 2.792155e-04
## [46] 4.304630e-04 3.048869e-04

munic$geometry

## Geometry set for 1121 features 
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -81.73562 ymin: -4.229406 xmax: -66.84722 ymax: 13.39473
## Geodetic CRS:  MAGNA-SIRGAS
## First 5 geometries:

## MULTIPOLYGON (((-75.66974 6.373599, -75.66965 6...

## MULTIPOLYGON (((-75.46938 5.94575, -75.46897 5....

## MULTIPOLYGON (((-76.08351 6.750504, -76.08325 6...

## MULTIPOLYGON (((-75.0332 6.415864, -75.03313 6....

## MULTIPOLYGON (((-75.67587 6.085613, -75.6754 6....

Cambiamos a usar centersxycenters y en lugar de centers$geometry en el código de la visualización con Leaflet por varias razones: La principal es por la compativilidad con ciertos paquetes o funciones que requieren coordenadas numéricas en lugar de objetos geométricos completos a demas de que proporciona una forma más simple y clara de acceder a las coordenadas de los centroides de las regiones. Esto hace que el código sea más legible y fácil de entender, especialmente para quienes puedan no estar familiarizados con la estructura de datos de geometría de sf

(centers <- st_centroid(munic1))

## Warning: st_centroid assumes attributes are constant over geometries

centers$x = st_coordinates(centers)[,1]

centers$x

##  [1] -75.25258 -74.94078 -74.98631 -74.80442 -75.19117 -74.84750 -75.48229
##  [8] -75.49746 -74.74700 -75.18906 -75.58985 -74.91148 -75.14854 -74.68937
## [15] -74.80258 -74.89353 -74.95704 -74.83709 -75.05229 -74.97657 -75.24323
## [22] -74.53964 -74.92337 -75.04763 -74.90668 -74.63783 -75.21994 -75.12202
## [29] -75.27824 -75.01845 -74.92091 -75.81685 -74.87248 -74.87451 -75.85482
## [36] -75.59428 -75.34759 -75.01860 -75.49909 -75.11245 -75.20713 -74.81858
## [43] -75.17148 -74.92089 -75.15583 -74.61897 -74.78303

centers$y = st_coordinates(centers)[,2]

centers$y

##  [1] 4.451919 3.390013 4.582622 4.824907 4.631786 5.006757 3.428434 4.406966
##  [9] 4.123414 5.016691 3.743690 4.354156 3.738624 3.982521 3.597861 4.166477
## [17] 5.079255 4.242813 5.186693 4.076498 5.068553 4.133273 4.866033 4.877932
## [25] 5.235322 4.186067 4.826735 3.540516 3.937790 5.094035 4.485001 3.098974
## [33] 3.729167 3.854282 3.468064 4.097858 4.216730 3.912542 3.955140 4.128605
## [41] 4.729114 4.048927 4.182124 4.709882 4.965859 3.854379 5.179411

munic1$Long <- centers$x 

munic1$Lat <- centers$y 

library(leaflet)

# Crear el mapa
map <- leaflet(munic1) %>%
  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", munic1$km2)(munic1$km2),
    highlightOptions = highlightOptions(
      color = "white", weight = 2,
      bringToFront = TRUE
    )
  ) %>%
  
  # Agregar marcadores (puntos)
  addMarkers(
    lng = munic1$Long, lat = munic1$Lat,
    popup = ~munic1$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'

# Mostrar el mapa
map

library(sf)

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

# Encontrar los centroides de las regiones
centers <- st_centroid(sp.munic1)

## 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.munic1)
centers$munic1 <- munic1$MPIO_CNMBR

centers

library(elevatr)
library(raster)

-Ahora, descarguemos los datos de elevación de nuestro departamento:

# z denota el nivel de zoom de los datos
# cuanto mayor sea el zoom, mayor será la resolución espacial
elevation <- get_elev_raster(munic1, z = 10)

## Mosaicing & Projecting

## Note: Elevation units are in meters.

-¿Qué es el objeto de elevación?

elevation

## class      : RasterLayer 
## dimensions : 4092, 3078, 12595176  (nrow, ncol, ncell)
## resolution : 0.000685414, 0.000685414  (x, y)
## extent     : -76.28906, -74.17936, 2.811272, 5.615986  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
## source     : file563469471e2c.tif 
## names      : file563469471e2c

library(sf)

# Crear un objeto sf a partir de munic
sf_munic1 <- st_as_sf(munic1, coords = c("Long", "Lat"))

# Obtener los centroides de las regiones
centers <- st_centroid(sf_munic1)

## 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_munic1)
centers$munic1 <- munic1$MPIO_CNMBR

# Imprimir el resultado
print(centers)

## Simple feature collection with 47 features and 8 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -75.85482 ymin: 3.098974 xmax: -74.53964 ymax: 5.235322
## Geodetic CRS:  MAGNA-SIRGAS
## First 10 features:
##           MPIO_CNMBR MPIO_CDPMP MPIO_NAREA                   geometry
## 1             IBAGUÉ      73001  1377.2514 POINT (-75.25258 4.451919)
## 2          ALPUJARRA      73024   501.5113 POINT (-74.94078 3.390013)
## 3           ALVARADO      73026   344.3356 POINT (-74.98631 4.582622)
## 4           AMBALEMA      73030   238.0448 POINT (-74.80442 4.824907)
## 5         ANZOÁTEGUI      73043   469.7113 POINT (-75.19117 4.631786)
## 6             ARMERO      73055   440.6308  POINT (-74.8475 5.006757)
## 7              ATACO      73067  1017.6311 POINT (-75.48229 3.428434)
## 8          CAJAMARCA      73124   508.4909 POINT (-75.49746 4.406966)
## 9  CARMEN DE APICALÁ      73148   190.8361   POINT (-74.747 4.123414)
## 10        CASABIANCA      73152   174.9306 POINT (-75.18906 5.016691)
##             km2      Long      Lat region            munic1
## 1  0.0013772514 -75.25258 4.451919      1            IBAGUÉ
## 2  0.0005015113 -74.94078 3.390013      2         ALPUJARRA
## 3  0.0003443356 -74.98631 4.582622      3          ALVARADO
## 4  0.0002380448 -74.80442 4.824907      4          AMBALEMA
## 5  0.0004697113 -75.19117 4.631786      5        ANZOÁTEGUI
## 6  0.0004406308 -74.84750 5.006757      6            ARMERO
## 7  0.0010176311 -75.48229 3.428434      7             ATACO
## 8  0.0005084909 -75.49746 4.406966      8         CAJAMARCA
## 9  0.0001908361 -74.74700 4.123414      9 CARMEN DE APICALÁ
## 10 0.0001749306 -75.18906 5.016691     10        CASABIANCA

4. Visualizar los datos de elevación

-A efectos de visualización, una resolución más baja acelera la tarea:

# Mosaico de los datos de elevación a una resolución más baja (opcional)
elevation2 <- aggregate(elevation, fact = 4, fun = mean)

sf_munic1 <- st_transform(sf_munic1, crs = 4326)

-Una buena paleta es clave para una visualización adecuada. Probemos este:

pal <- colorNumeric(c("forestgreen","yellow","tan","brown"), values(elevation2),
  na.color = "transparent")

# Supongamos que deseas proyectar tus datos a EPSG:4326
sf_munic1 <- st_transform(sf_munic1, crs = 4326)
centers <- st_transform(centers, crs = 4326)

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

leaflet(munic1) %>% 
  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 = ~Long, lat = ~Lat,  # 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 = "Elevation data for Tolima (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'

-Haga clic en el ID de cada región para obtener el nombre y la extensión del municipio en Km^2

5. 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 dimensiones horizontales son grados decimales, PERO la unidad de dimensión vertical son metros.

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

Primero, definimos el CRS objetivo:

#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é resulta?

rep_elev

## class      : RasterLayer 
## dimensions : 3877, 2931, 11363487  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 356719, 591199, 310752.3, 620912.3  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : r_tmp_2023-11-06_171827.329719_22068_65243.grd 
## names      : layer 
## values     : -148.4632, 5370.803  (min, max)

-Ahora reproyectamos el SpatialPolygonsDataFrame que representa los municipios de nuestro departamento.

library(sf)

# Transformar el objeto sf utilizando st_transform
rep_munic1 <- st_transform(sp.munic1, crs = spatialref)

rep_munic1

-Para evitar errores, 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).

library(raster)

# Ruta de destino para guardar el DEM reproyectado
filename <- "./rep_tolima_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)

# Una vez que hayas ejecutado este código, puedes comentarlo nuevamente

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

-Ahora, creemos el histograma de valores de elevación:

hist(rep_elev)

-Estadísticas de 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')

-Creemos un vector de estadísticas unidimensional:

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

-El resultado:

(df_estadisticas <- data.frame(metricas, valores))

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

-Comprobemos el primer resultado:

slope

## class      : RasterLayer 
## dimensions : 3877, 2931, 11363487  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 356719, 591199, 310752.3, 620912.3  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : slope 
## values     : 0, 87.29878  (min, max)

-Ahora, el segundo resultado:

aspect

## class      : RasterLayer 
## dimensions : 3877, 2931, 11363487  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 356719, 591199, 310752.3, 620912.3  (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:

#Este trozo es opcional
#Úselo solo cuando los datos de zoom sean muy altos
slope2 <- aggregate(slope, fact=4, fun=mean)

-Puede resultar conveniente proyectar el conjunto de datos de elevación:

slope3 <- projectRasterForLeaflet(slope2, "ngb")

-Creemos una 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

-hora el codigo para el gráfico:

leaflet(munic1) %>% addTiles() %>%
    addPolygons(color = "white", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.25, fillOpacity = 0.15,
    popup = paste("Region: ", munic1$MPIO_CNMBR, "<br>",
                          "Km2: ", munic1$km2, "<br>")) %>%
  addLabelOnlyMarkers(data = centers,
                   lng = ~Long, lat = ~Lat, 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 Tolima (%)")

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

-Recordemos qué es el objeto pendiente2:

slope2

## class      : RasterLayer 
## dimensions : 970, 733, 711010  (nrow, ncol, ncell)
## resolution : 320, 320  (x, y)
## extent     : 356719, 591279, 310512.3, 620912.3  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : slope 
## values     : 0.0671791, 85.27139  (min, max)

-Tenga en cuenta que existen dos bibliotecas para administrar datos ráster en R: - raster y - terra.

-La biblioteca ráster 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 ráster) en un SpatRaster (es decir, un objeto terra):

(nslope2  <- as(slope2, "SpatRaster"))

## class       : SpatRaster 
## dimensions  : 970, 733, 1  (nrow, ncol, nlyr)
## resolution  : 320, 320  (x, y)
## extent      : 356719, 591279, 310512.3, 620912.3  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        :      slope 
## min value   :  0.0671791 
## max value   : 85.2713875

-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

En tu cuaderno, explica el significado de la matriz anterior.

La matriz que se ha creado se utiliza para reclasificar intervalos de pendientes. Cada fila de esta matriz representa un rango de pendiente y consta de tres columnas. Aquí está el significado de cada columna:

-Columna 1: Este valor representa el límite inferior del rango de pendiente. -Columna 2: Este valor es el límite superior del rango de pendiente. -Columna 3: Esta columna contiene un número que se usa para clasificar el rango de pendiente. Es decir, cuando la pendiente de una región cae dentro de un rango especificado por las columnas 1 y 2, se le asigna el valor correspondiente de la columna 3 como su clase de pendiente. Entonces, por ejemplo, la primera fila de la matriz (0, 3, 1) significa que las pendientes que van desde 0 hasta 3 se clasificarán como “Clase 1”. De manera similar, las siguientes filas definen otros rangos de pendiente y sus respectivas clases.

-Ahora, la tarea de reclasificación:

# para valores >= 0 (en lugar de > 0)
slope_rec <- terra::classify(nslope2, rclmat, right=TRUE)

-Veamos los resultados:

slope_rec

## class       : SpatRaster 
## dimensions  : 970, 733, 1  (nrow, ncol, nlyr)
## resolution  : 320, 320  (x, y)
## extent      : 356719, 591279, 310512.3, 620912.3  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        :    slope 
## min value   :  1.00000 
## max value   : 85.27139

8. Visualización estática

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

Rectificamos la libreria necesaria:

library(tmap)

## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')

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 = "pretty", title = intervals, 
            palette = colors, legend.show = TRUE) +  
  tm_shape(munic1) + 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)

-Veamos el resultado:

slope.map

• Creemos un respaldo:

## Tenga en cuenta que estamos utilizando papel de tamaño carta.
tmap_save(slope.map, "./tolima_slope.png", width = 8.5, height =11, dpi=600)

## Map saved to C:\Users\camil\OneDrive\Escritorio\Geomatica\CuadernoEVA\tolima_slope.png

## Resolution: 5100 by 6600 pixels

## Size: 8.5 by 11 inches (600 dpi)

9. Visualización interactiva

-Ahora queremos crear un mapa de pendientes interactivo. Convirtamos el objeto sf en un objeto SpatVector:

(nmunic1  <- as(rep_munic1, "SpatVector"))

##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 47, 6  (geometries, attributes)
##  extent      : 377130.7, 558289, 317378.1, 587966.7  (xmin, xmax, ymin, ymax)
##  coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
##  names       : MPIO_CNMBR MPIO_CDPMP MPIO_NAREA       km2   Long   Lat
##  type        :      <chr>      <chr>      <num>     <num>  <num> <num>
##  values      :     IBAGUÉ      73001       1377  0.001377 -75.25 4.452
##                 ALPUJARRA      73024      501.5 0.0005015 -74.94  3.39
##                  ALVARADO      73026      344.3 0.0003443 -74.99 4.583

-Convirtamos el objeto ráster en un objeto terra:

(nslope  <- as(slope, "SpatRaster"))

## class       : SpatRaster 
## dimensions  : 3877, 2931, 1  (nrow, ncol, nlyr)
## resolution  : 80, 80  (x, y)
## extent      : 356719, 591199, 310752.3, 620912.3  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        :    slope 
## min value   :  0.00000 
## max value   : 87.29878

-Ahora, calculemos las categorías de pendientes agregadas por municipio:

## Método S4 para la firma 'SpatRaster,SpatVector'
munic1.slope <- zonal(nslope, nmunic1, fun=mean, as.raster=TRUE)

-Ahora, veamos el mapa:

tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(munic1) + tm_borders(col = "darkgray", lwd = 0.5, alpha=0.8) +
  tm_text(text = "MPIO_CNMBR", size = 0.5, alpha=0.7) +
  tm_shape(munic1.slope) +
  tm_raster(style="cont", n=7, alpha=0.6, title="Pendiente media")+
tm_layout(legend.outside = TRUE)

## stars object downsampled to 870 by 1150 cells. See tm_shape manual (argument raster.downsample)

10. Bibliografía

LizsessionInfo()23. Procesamiento y análisis de datos de elevación en R. Disponible en https://rpubs.com/ials2un/dem_analysis Puede encontrar este cuaderno en el siguiente link <http://rpubs.com/alejandrotrian/DEM #Reproducibilidad

sessionInfo()

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