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 nosotors, 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:

• Escriba y ejecute cada fragmento paso a paso
• En caso de errores, lee el mensaje, interpretalo e intenta solucionarlo.
• En caso de que el error persista, consultar al profesor antes de la fecha límite.
• Una vez que se produce su archivo html, es hora de publicarlo en RPubs
  
• Después de publicar el cuaderno, revíselo, corrija la redacción de mala calidad y vuelva a publicarlo usando la opción de actualización.
• En caso de problemas, revise la documentación de R o busque que las librerias sean compatibles con su version de Rstudio.

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)
library(dplyr)
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] "huila_slope.png"                     "Modelo de Elevacion Digital.nb.html"
##  [3] "Modelo de Elevacion Digital.Rmd"     "Modelo_de_Elevacion_Digital.html"   
##  [5] "Modelo_de_Elevacion_Digital.nb.html" "Modelo_de_Elevacion_Digital.Rmd"    
##  [7] "NUEVOHUILA.qgz"                      "NUEVOS_HUILA.cpg"                   
##  [9] "NUEVOS_HUILA.dbf"                    "NUEVOS_HUILA.prj"                   
## [11] "NUEVOS_HUILA.qmd"                    "NUEVOS_HUILA.shp"                   
## [13] "NUEVOS_HUILA.shx"                    "rep_huila_dem.tif"                  
## [15] "rsconnect"

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

## Reading layer `NUEVOS_HUILA' from data source 
##   `D:\Sexta Matricula\Geomatica\mapa DEM\NUEVOS_HUILA.shp' using driver `ESRI Shapefile'
## Simple feature collection with 37 features and 7 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.62466 ymin: 1.552125 xmax: -74.41303 ymax: 3.843208
## Geodetic CRS:  MAGNA-SIRGAS

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

head(munic)

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

munic$km2 = munic$AREA / 1000000

munic$km2

##  [1] 1269.82717  503.70310  249.60734  194.23340  356.51548  370.65657
##  [7]  356.42110  814.65092  131.94003  155.81261  278.31076  884.26636
## [13]  178.21746  193.66005  630.66594  251.23017  466.64491 1390.43219
## [19]  339.37959  429.66717  362.81184  367.02539  530.39939  471.73435
## [25]  185.72709  543.92692  332.91538  606.52117   80.44174 1584.96992
## [31]  462.79934  786.00202  180.97556  589.68899  795.63062  273.73165
## [37]  540.51687

munic$geometry

## Geometry set for 37 features 
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.62466 ymin: 1.552125 xmax: -74.41303 ymax: 3.843208
## Geodetic CRS:  MAGNA-SIRGAS
## First 5 geometries:

## POLYGON ((-75.35274 3.282034, -75.35222 3.28198...

## POLYGON ((-75.52383 2.512427, -75.52346 2.51242...

## POLYGON ((-75.70495 2.08465, -75.70472 2.084607...

## POLYGON ((-75.44193 2.677073, -75.44132 2.67665...

## POLYGON ((-75.67548 2.797141, -75.6749 2.796984...

Cambiamos a usar centers\(x y centers\)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(munic))

## Warning: st_centroid assumes attributes are constant over geometries

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

centers$x

##  [1] -75.27233 -75.51700 -75.69694 -75.44781 -75.67383 -76.25636 -76.11782
##  [8] -76.00118 -75.78694 -76.04110 -75.72803 -75.44066 -76.14786 -75.83195
## [15] -76.11920 -75.23554 -76.22897 -76.41656 -75.63360 -75.80798 -75.87535
## [22] -75.68260 -75.08536 -75.71434 -75.91758 -75.13839 -75.51593 -75.57131
## [29] -75.95092 -74.67833 -75.32944 -74.95797 -75.76858 -75.29136 -75.35263
## [36] -75.71632 -75.99665

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

centers$y

##  [1] 2.993357 2.371780 1.984393 2.541653 2.691515 2.006780 2.181456 2.328075
##  [9] 2.576064 2.061191 2.406899 2.914488 1.684580 2.259601 1.827902 2.782224
## [17] 2.100383 1.924161 2.926266 1.870991 2.130281 2.526854 3.030613 2.824119
## [25] 1.950715 3.287054 2.646561 2.172440 2.006248 3.451956 2.658666 3.135422
## [33] 2.079228 2.497847 3.262372 2.269749 1.709113

munic$Long <- centers$x 

munic$Lat <- centers$y 

library(leaflet)

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

# Mostrar el mapa
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

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(munic, z = 10)

## Mosaicing & Projecting

## Note: Elevation units are in meters.

¿Qué es el objeto de elevación?

elevation

## class      : RasterLayer 
## dimensions : 3582, 3586, 12845052  (nrow, ncol, ncell)
## resolution : 0.0006862561, 0.0006862561  (x, y)
## extent     : -76.64062, -74.17971, 1.406085, 3.864255  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
## source     : file278069ef559.tif 
## names      : file278069ef559

library(sf)

# 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 37 features and 12 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -76.41656 ymin: 1.68458 xmax: -74.67833 ymax: 3.451956
## Geodetic CRS:  MAGNA-SIRGAS
## First 10 features:
##    DPTO_CCDGO MPIO_CCDGO   MPIO_CNMBR MPIO_CDPMP       AREA  LATITUD  LONGITUD
## 1          41        001        NEIVA      41001 1269827167 2.993360 -75.27236
## 2          41        306      GIGANTE      41306  503703102 2.371782 -75.51700
## 3          41        319    GUADALUPE      41319  249607339 1.984394 -75.69694
## 4          41        349         HOBO      41349  194233398 2.541654 -75.44781
## 5          41        357       ÍQUIRA      41357  356515480 2.691518 -75.67383
## 6          41        359        ISNOS      41359  370656566 2.006784 -76.25636
## 7          41        378 LA ARGENTINA      41378  356421101 2.181453 -76.11782
## 8          41        396     LA PLATA      41396  814650921 2.328077 -76.00117
## 9          41        483       NÁTAGA      41483  131940031 2.576066 -75.78694
## 10         41        503      OPORAPA      41503  155812607 2.061191 -76.04110
##                      geometry       km2      Long      Lat region        munic
## 1  POINT (-75.27233 2.993357) 1269.8272 -75.27233 2.993357      1        NEIVA
## 2     POINT (-75.517 2.37178)  503.7031 -75.51700 2.371780      2      GIGANTE
## 3  POINT (-75.69694 1.984393)  249.6073 -75.69694 1.984393      3    GUADALUPE
## 4  POINT (-75.44781 2.541653)  194.2334 -75.44781 2.541653      4         HOBO
## 5  POINT (-75.67383 2.691515)  356.5155 -75.67383 2.691515      5       ÍQUIRA
## 6   POINT (-76.25636 2.00678)  370.6566 -76.25636 2.006780      6        ISNOS
## 7  POINT (-76.11782 2.181456)  356.4211 -76.11782 2.181456      7 LA ARGENTINA
## 8  POINT (-76.00118 2.328075)  814.6509 -76.00118 2.328075      8     LA PLATA
## 9  POINT (-75.78694 2.576064)  131.9400 -75.78694 2.576064      9       NÁTAGA
## 10  POINT (-76.0411 2.061191)  155.8126 -76.04110 2.061191     10      OPORAPA

4. Visualizar los datos de elevacion

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_munic <- st_transform(sf_munic, 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_munic <- st_transform(sf_munic, crs = 4326)
centers <- st_transform(centers, crs = 4326)

Ahora 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 = "Elevation data for Huila (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)

Que resulta ?

rep_elev

## class      : RasterLayer 
## dimensions : 3399, 3422, 11631378  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 317469, 591229, 155378.5, 427298.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : file278069ef559 
## values     : -666.9117, 5367.665  (min, max)

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

library(sf)

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

rep_munic

To avoid a headache, let’s save our DEM. Make sure to indicate a pathname and a filename which suits your needs (i.e. appropriate for your platform, and your department).

library(raster)

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

Podemos calcular un resumen de las 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 : 3399, 3422, 11631378  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 317469, 591229, 155378.5, 427298.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : slope 
## values     : 0, 87.43722  (min, max)

Ahora, el segundo resultado:

aspect

## class      : RasterLayer 
## dimensions : 3399, 3422, 11631378  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 317469, 591229, 155378.5, 427298.5  (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

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 Huila (%)")

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

Recordemos qué es el objeto pendiente2:

slope2

## class      : RasterLayer 
## dimensions : 850, 856, 727600  (nrow, ncol, ncell)
## resolution : 320, 320  (x, y)
## extent     : 317469, 591389, 155298.5, 427298.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : slope 
## values     : 1.849426e-14, 84.72774  (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  : 850, 856, 1  (nrow, ncol, nlyr)
## resolution  : 320, 320  (x, y)
## extent      : 317469, 591389, 155298.5, 427298.5  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        :        slope 
## min value   : 1.849426e-14 
## max value   : 8.472774e+01

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.

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  : 850, 856, 1  (nrow, ncol, nlyr)
## resolution  : 320, 320  (x, y)
## extent      : 317469, 591389, 155298.5, 427298.5  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        :    slope 
## min value   :  1.00000 
## max value   : 84.72774

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

# Crear intervalos de clase personalizados
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)

Veamos el resultado:

slope.map

Esta guía contiene más consejos para trazar objetos raster y vectoriales. Asegúrate de probar opciones adicionales a las incluidas en este cuaderno.

Creemos un respaldo:

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

## Map saved to D:\Sexta Matricula\Geomatica\mapa DEM\huila_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

9. 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  : 37, 10  (geometries, attributes)
##  extent      : 319297.7, 565179, 171589.9, 424810.4  (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      :         41        001      NEIVA      41001  1.27e+09   2.993
##                        41        306    GIGANTE      41306 5.037e+08   2.372
##                        41        319  GUADALUPE      41319 2.496e+08   1.984
##  LONGITUD   km2   Long   Lat
##     <num> <num>  <num> <num>
##    -75.27  1270 -75.27 2.993
##    -75.52 503.7 -75.52 2.372
##     -75.7 249.6  -75.7 1.984

Convirtamos el objeto ráster en un objeto terra:

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

## class       : SpatRaster 
## dimensions  : 3399, 3422, 1  (nrow, ncol, nlyr)
## resolution  : 80, 80  (x, y)
## extent      : 317469, 591229, 155378.5, 427298.5  (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.43722

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

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

Ahora, 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 = "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 1004 by 997 cells. See tm_shape manual (argument raster.downsample)

10. Bibliografía.

Lizarazo, I., 2023. 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.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] 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-8   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.1           rlang_1.1.1             
## [13] sass_0.4.7               progress_1.2.2           tools_4.3.1             
## [16] utf8_1.2.3               yaml_2.3.7               knitr_1.44              
## [19] prettyunits_1.1.1        htmlwidgets_1.6.2        interp_1.1-4            
## [22] classInt_0.4-10          curl_5.0.2               RColorBrewer_1.1-3      
## [25] abind_1.4-5              KernSmooth_2.23-21       purrr_1.0.2             
## [28] leafsync_0.1.0           grid_4.3.1               stats4_4.3.1            
## [31] fansi_1.0.4              latticeExtra_0.6-30      e1071_1.7-13            
## [34] leafem_0.2.3             colorspace_2.1-0         progressr_0.14.0        
## [37] scales_1.2.1             dichromat_2.0-0.1        cli_3.6.1               
## [40] rmarkdown_2.25           crayon_1.5.2             generics_0.1.3          
## [43] rstudioapi_0.15.0        httr_1.4.7               tzdb_0.4.0              
## [46] tmaptools_3.1-1          DBI_1.1.3                cachem_1.0.8            
## [49] proxy_0.4-27             stars_0.6-4              slippymath_0.3.1        
## [52] parallel_4.3.1           s2_1.1.4                 base64enc_0.1-3         
## [55] vctrs_0.6.3              jsonlite_1.8.7           hms_1.1.3               
## [58] jpeg_0.1-10              crosstalk_1.2.0          jquerylib_0.1.4         
## [61] hexbin_1.28.3            units_0.8-4              lwgeom_0.2-13           
## [64] glue_1.6.2               leaflet.providers_1.13.0 codetools_0.2-19        
## [67] deldir_1.0-9             munsell_0.5.0            tibble_3.2.1            
## [70] pillar_1.9.0             htmltools_0.5.6          satellite_1.0.4         
## [73] R6_2.5.1                 wk_0.8.0                 evaluate_0.21           
## [76] png_0.1-8                bslib_0.5.1              class_7.3-22            
## [79] Rcpp_1.0.11              xfun_0.40                zoo_1.8-12              
## [82] pkgconfig_2.0.3