1. Introducción.

En este cuaderno se desarrolla varios códigos, los cuales son útiles para poder analizar la elevación de un territorio y sus accidentes geográfiocos, también ayuda a determinar otros factores, que se verán a continuación.

2. Código inicial.

library(rasterVis)
## Loading required package: lattice
library(raster)
## Loading required package: sp
library(rgl)
library(rgdal)
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.1, released 2021/12/27
## Path to GDAL shared files: C:/Users/GERARDO/Documents/R/win-library/4.1/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/GERARDO/Documents/R/win-library/4.1/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-7
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(elevatr)
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.6     v dplyr   1.0.8
## v tidyr   1.2.0     v stringr 1.4.0
## v readr   2.1.2     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::extract() masks raster::extract()
## x dplyr::filter()  masks stats::filter()
## x dplyr::lag()     masks stats::lag()
## x dplyr::select()  masks raster::select()
library(mapview)

3. Obtener datos de elevación raster

(munic <-  st_read("C:/Users/GERARDO/Desktop/GB2022/datos/11_BOGOTA/11_BOGOTA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))
## Reading layer `MGN_MPIO_POLITICO' from data source 
##   `C:\Users\GERARDO\Desktop\GB2022\Datos\11_BOGOTA\11_BOGOTA\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.45085 ymin: 3.730633 xmax: -73.98613 ymax: 4.836827
## Geodetic CRS:  WGS 84
## Simple feature collection with 1 feature and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.45085 ymin: 3.730633 xmax: -73.98613 ymax: 4.836827
## Geodetic CRS:  WGS 84
##   DPTO_CCDGO MPIO_CCDGO   MPIO_CNMBR MPIO_CRSLC MPIO_NAREA MPIO_NANO
## 1         11      11001 BOGOTÁ, D.C.       1538   1622.853      2017
##     DPTO_CNMBR Shape_Leng Shape_Area                       geometry
## 1 BOGOTÁ, D.C.   3.760453  0.1322079 POLYGON ((-74.07419 4.83654...
head(munic)
## Simple feature collection with 1 feature and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.45085 ymin: 3.730633 xmax: -73.98613 ymax: 4.836827
## Geodetic CRS:  WGS 84
##   DPTO_CCDGO MPIO_CCDGO   MPIO_CNMBR MPIO_CRSLC MPIO_NAREA MPIO_NANO
## 1         11      11001 BOGOTÁ, D.C.       1538   1622.853      2017
##     DPTO_CNMBR Shape_Leng Shape_Area                       geometry
## 1 BOGOTÁ, D.C.   3.760453  0.1322079 POLYGON ((-74.07419 4.83654...
head(munic)
## Simple feature collection with 1 feature and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.45085 ymin: 3.730633 xmax: -73.98613 ymax: 4.836827
## Geodetic CRS:  WGS 84
##   DPTO_CCDGO MPIO_CCDGO   MPIO_CNMBR MPIO_CRSLC MPIO_NAREA MPIO_NANO
## 1         11      11001 BOGOTÁ, D.C.       1538   1622.853      2017
##     DPTO_CNMBR Shape_Leng Shape_Area                       geometry
## 1 BOGOTÁ, D.C.   3.760453  0.1322079 POLYGON ((-74.07419 4.83654...
elevation <- get_elev_raster(munic, z = 10)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
if (require(rgdal)) {
 rf <- writeRaster(elevation, filename=file.path("C:/Users/GERARDO/Desktop/GB2022/datos/bogota_dem_z10.tif"), format="GTiff", overwrite=TRUE)
}
elevation
## class      : RasterLayer 
## dimensions : 2047, 1026, 2100222  (nrow, ncol, ncell)
## resolution : 0.0006851466, 0.0006851466  (x, y)
## extent     : -74.53125, -73.82829, 3.513338, 4.915833  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : file14a444ea2865.tif 
## names      : file14a444ea2865 
## values     : -32768, 32767  (min, max)
munic_sp <- as(munic, "Spatial")
plot(elevation, main="Este es el DEM descargado [metros]")
plot(munic_sp,col="NA",border="black", add=TRUE)
text(coordinates(munic_sp), labels=as.character(munic$MPIO_CNMBR), 
     col="black", cex=0.05)

4. Recolectar los datos de elevación para que coincidan con la extensión del área de estudio

elev_crop = crop(elevation, munic_sp)
plot(elev_crop, main="Modelo de elevación digital recortado")
plot(munic_sp, add=TRUE)
text(coordinates(munic_sp), labels=as.character(munic_sp$MPIO_CNMBR), cex=0.2)

elev_crop
## class      : RasterLayer 
## dimensions : 1615, 679, 1096585  (nrow, ncol, ncell)
## resolution : 0.0006851466, 0.0006851466  (x, y)
## extent     : -74.45109, -73.98587, 3.730529, 4.837041  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : file14a444ea2865 
## values     : -565, 4255  (min, max)

5. Reproyectar los datos de elevación.

crs(elev_crop)
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["WGS 84 (with axis order normalized for visualization)",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]
library(sp)
spatialref <- CRS(SRS_string="EPSG:32618")
pr3 <- projectExtent(elev_crop, crs=spatialref)
res(pr3) <- 180
rep_elev <- projectRaster(elev_crop, pr3)
rep_elev
## class      : RasterLayer 
## dimensions : 680, 288, 195840  (nrow, ncol, ncell)
## resolution : 180, 180  (x, y)
## extent     : 560865, 612705, 412335, 534735  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : file14a444ea2865 
## values     : -258.1412, 4200.682  (min, max)
(rep_munic = spTransform(munic_sp,spatialref))
## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 560977.7, 612459.1, 412372.7, 534697.2  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## variables   : 9
## names       : DPTO_CCDGO, MPIO_CCDGO,   MPIO_CNMBR, MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO,   DPTO_CNMBR,    Shape_Leng,     Shape_Area 
## value       :         11,      11001, BOGOTÁ, D.C.,       1538, 1622.85261422,      2017, BOGOTÁ, D.C., 3.76045332576, 0.132207854177
plot(rep_elev, main="Modelo de elevación digital reproyectado")
plot(rep_munic, lwd=0.5, add=TRUE)
text(coordinates(rep_munic), labels=as.character(rep_munic$MPIO_CNMBR), cex=0.2)

writeRaster(rep_elev, filename="C:/Users/GERARDO/Desktop/GB2022/datos/rep_boyaca_dem.tif", datatype='INT4S', overwrite=TRUE)

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

hist(rep_elev)

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)
(df_estadisticas <- data.frame(metricas, valores))
##   metricas   valores
## 1     mean 2863.2863
## 2      min -258.1412
## 3      max 4200.6818
## 4      std  671.8537

7. Obtención de variables geomorfométricas.

slope = terrain(rep_elev,opt='slope', unit='degrees')
aspect = terrain(rep_elev,opt='aspect',unit='degrees')
hill = hillShade(slope,aspect,40,315)
plot(rep_elev,main="DEM para Bogotá [metros]", col=terrain.colors(25,alpha=0.8))
plot(rep_munic, add=TRUE, lwd=0.5)
text(coordinates(rep_munic), labels=as.character(rep_munic$MPIO_CNMBR), cex=0.2)

plot(slope,main="Pendiente para Bogtotá [grados]", col=topo.colors(6,alpha=0.6))
plot(rep_munic, add=TRUE, lwd=0.5)
text(coordinates(rep_munic), labels=as.character(rep_munic$MPIO_CNMBR), cex=0.2)

plot(aspect,main="Aspecto para la extendión de Bogotá [grados]", col=rainbow(10,alpha=0.7))
plot(rep_munic, add=TRUE, lwd=0.5)
text(coordinates(rep_munic), labels=as.character(rep_munic$MPIO_CNMBR), cex=0.2)

plot(hill,
        col=grey(1:100/100),  
        legend=FALSE,         
        main="DEM for Bogota [m]",
        axes=FALSE)           

plot(rep_elev, 
        axes=FALSE,
        col=terrain.colors(12, alpha=0.6), add=TRUE) # color method
         # sets how transparent the object will be (0=transparent, 1=not transparent)
plot(rep_munic, add=TRUE, lwd=0.5)
text(coordinates(rep_munic), labels=as.character(rep_munic$MPIO_CNMBR), cex=0.2)

8. Mapeo de datos de elevación con Rayshader.

library(rayshader)
elmat = raster_to_matrix(rep_elev)
#summary(elmat)
elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  #add_shadow(ray_shade(elmat), 0.5) %>%
  #add_shadow(ambient_shade(elmat), 0) %>%
  plot_map()

9. Bibliografía

“Lizarazo, I., 2022. My fifth notebook: Elevation data in R. Available at https://rpubs.com/ials2un/boyaca_dem