options(repos = c(CRAN = "https://cran.rstudio.com/"))
rm(list=ls())
# Install packages before load them
#Para Mac se necesita Xquartz (aplicaciòn).
#Verificar si están instalados los paquetes. Si no, instalar... uncomment if needed
#install.packages(c('elevatr', 'sf', 'leaflet', 'terra', 'MultiscaleDTM', 'exactextractr'))
# Cargar las librerias
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(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(leaflet)
library(terra)
## terra 1.8.15
library(MultiscaleDTM) #to use this one, it´s needed to install rlg package
library(exactextractr)
# Cargar capa
munic <- st_read("~/Desktop/cordoba/Cordoba.gpkg")
## Reading layer `cordoba' from data source
## `/Users/jsc_89/Desktop/cordoba/Cordoba.gpkg' using driver `GPKG'
## Simple feature collection with 30 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.51453 ymin: 7.347087 xmax: -74.78095 ymax: 9.447748
## Geodetic CRS: MAGNA-SIRGAS
# Conectar con Marco Geoestadistico Nacional de Referencia
(cordoba_crs = st_crs(st_geometry(munic)))
## Coordinate Reference System:
## User input: MAGNA-SIRGAS
## wkt:
## GEOGCRS["MAGNA-SIRGAS",
## DATUM["Marco Geocentrico Nacional de Referencia",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank."],
## BBOX[-4.23,-84.77,15.51,-66.87]],
## ID["EPSG",4686]]
# Area municipios de Córdoba
munic$area_m2 = sf::st_area(munic)
# el área expresada en kilómetros cuadrados
munic$area = munic$area_m2/1000000
#redondear el area a dos digitos
munic$area = signif(munic$area, digits = 6)
#llamar municipios
munic
## Simple feature collection with 30 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.51453 ymin: 7.347087 xmax: -74.78095 ymax: 9.447748
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 23 001 23001 CÓRDOBA MONTERÍA
## 2 23 068 23068 CÓRDOBA AYAPEL
## 3 23 079 23079 CÓRDOBA BUENAVISTA
## 4 23 090 23090 CÓRDOBA CANALETE
## 5 23 162 23162 CÓRDOBA CERETÉ
## 6 23 168 23168 CÓRDOBA CHIMÁ
## 7 23 182 23182 CÓRDOBA CHINÚ
## 8 23 189 23189 CÓRDOBA CIÉNAGA DE ORO
## 9 23 300 23300 CÓRDOBA COTORRA
## 10 23 350 23350 CÓRDOBA LA APARTADA
## mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1 Ordenanza 42 del 27 de Abril de 1923 MUNICIPIO 3136.8910 2023
## 2 Ordenanza 42 del 27 de Abril de 1923 MUNICIPIO 1964.8388 2023
## 3 Ordenanza 10 de Diciembre de 1969 MUNICIPIO 834.9270 2023
## 4 Ordenanza 28 de 1978 MUNICIPIO 420.0415 2023
## 5 1923 MUNICIPIO 290.8820 2023
## 6 1927 MUNICIPIO 323.9516 2023
## 7 Ordenanza 42 del 27 de Abril de 1923 MUNICIPIO 626.1762 2023
## 8 1923 MUNICIPIO 641.2572 2023
## 9 ORD 03 DE ABRIL 08 DE 1997 MUNICIPIO 88.1123 2023
## 10 ORD 07 DE MAYO 06 DE 1997 MUNICIPIO 287.0437 2023
## shape_Leng shape_Area geom area_m2
## 1 3.7233888 0.257383137 MULTIPOLYGON (((-76.12448 8... 3146708009 [m^2]
## 2 2.4096791 0.161210107 MULTIPOLYGON (((-75.00565 8... 1972548484 [m^2]
## 3 1.8888729 0.068472843 MULTIPOLYGON (((-75.35428 8... 837983885 [m^2]
## 4 1.3044687 0.034465103 MULTIPOLYGON (((-76.19951 8... 421208111 [m^2]
## 5 1.4926651 0.023889135 MULTIPOLYGON (((-75.75363 8... 291824384 [m^2]
## 6 0.9175294 0.026625648 MULTIPOLYGON (((-75.6027 9.... 325055945 [m^2]
## 7 1.7241529 0.051466850 MULTIPOLYGON (((-75.39083 9... 628471460 [m^2]
## 8 1.9769704 0.052669283 MULTIPOLYGON (((-75.59987 9... 643481366 [m^2]
## 9 0.6153369 0.007240021 MULTIPOLYGON (((-75.77023 9... 88402053 [m^2]
## 10 1.0054769 0.023535319 MULTIPOLYGON (((-75.29656 8... 288129833 [m^2]
## area
## 1 3146.7100 [m^2]
## 2 1972.5500 [m^2]
## 3 837.9840 [m^2]
## 4 421.2080 [m^2]
## 5 291.8240 [m^2]
## 6 325.0560 [m^2]
## 7 628.4710 [m^2]
## 8 643.4810 [m^2]
## 9 88.4021 [m^2]
## 10 288.1300 [m^2]
#nuevo objeto con centroides de municipios
(centers <- st_centroid(munic))
## Warning: st_centroid assumes attributes are constant over geometries
## Simple feature collection with 30 features and 13 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -76.29493 ymin: 7.711479 xmax: -75.0487 ymax: 9.35967
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 23 001 23001 CÓRDOBA MONTERÍA
## 2 23 068 23068 CÓRDOBA AYAPEL
## 3 23 079 23079 CÓRDOBA BUENAVISTA
## 4 23 090 23090 CÓRDOBA CANALETE
## 5 23 162 23162 CÓRDOBA CERETÉ
## 6 23 168 23168 CÓRDOBA CHIMÁ
## 7 23 182 23182 CÓRDOBA CHINÚ
## 8 23 189 23189 CÓRDOBA CIÉNAGA DE ORO
## 9 23 300 23300 CÓRDOBA COTORRA
## 10 23 350 23350 CÓRDOBA LA APARTADA
## mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1 Ordenanza 42 del 27 de Abril de 1923 MUNICIPIO 3136.8910 2023
## 2 Ordenanza 42 del 27 de Abril de 1923 MUNICIPIO 1964.8388 2023
## 3 Ordenanza 10 de Diciembre de 1969 MUNICIPIO 834.9270 2023
## 4 Ordenanza 28 de 1978 MUNICIPIO 420.0415 2023
## 5 1923 MUNICIPIO 290.8820 2023
## 6 1927 MUNICIPIO 323.9516 2023
## 7 Ordenanza 42 del 27 de Abril de 1923 MUNICIPIO 626.1762 2023
## 8 1923 MUNICIPIO 641.2572 2023
## 9 ORD 03 DE ABRIL 08 DE 1997 MUNICIPIO 88.1123 2023
## 10 ORD 07 DE MAYO 06 DE 1997 MUNICIPIO 287.0437 2023
## shape_Leng shape_Area geom area_m2
## 1 3.7233888 0.257383137 POINT (-75.9505 8.58464) 3146708009 [m^2]
## 2 2.4096791 0.161210107 POINT (-75.0487 8.265758) 1972548484 [m^2]
## 3 1.8888729 0.068472843 POINT (-75.43383 8.191389) 837983885 [m^2]
## 4 1.3044687 0.034465103 POINT (-76.23779 8.724093) 421208111 [m^2]
## 5 1.4926651 0.023889135 POINT (-75.83457 8.89102) 291824384 [m^2]
## 6 0.9175294 0.026625648 POINT (-75.65239 9.110253) 325055945 [m^2]
## 7 1.7241529 0.051466850 POINT (-75.35666 9.027499) 628471460 [m^2]
## 8 1.9769704 0.052669283 POINT (-75.60778 8.841785) 643481366 [m^2]
## 9 0.6153369 0.007240021 POINT (-75.77888 9.056992) 88402053 [m^2]
## 10 1.0054769 0.023535319 POINT (-75.28817 8.052216) 288129833 [m^2]
## area
## 1 3146.7100 [m^2]
## 2 1972.5500 [m^2]
## 3 837.9840 [m^2]
## 4 421.2080 [m^2]
## 5 291.8240 [m^2]
## 6 325.0560 [m^2]
## 7 628.4710 [m^2]
## 8 643.4810 [m^2]
## 9 88.4021 [m^2]
## 10 288.1300 [m^2]
#dividir cada uno en sus coordenadas x, y:
centers$x = st_coordinates(centers)[,1]
centers$y = st_coordinates(centers)[,2]
#shora se comprueba
centers
## Simple feature collection with 30 features and 15 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -76.29493 ymin: 7.711479 xmax: -75.0487 ymax: 9.35967
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 23 001 23001 CÓRDOBA MONTERÍA
## 2 23 068 23068 CÓRDOBA AYAPEL
## 3 23 079 23079 CÓRDOBA BUENAVISTA
## 4 23 090 23090 CÓRDOBA CANALETE
## 5 23 162 23162 CÓRDOBA CERETÉ
## 6 23 168 23168 CÓRDOBA CHIMÁ
## 7 23 182 23182 CÓRDOBA CHINÚ
## 8 23 189 23189 CÓRDOBA CIÉNAGA DE ORO
## 9 23 300 23300 CÓRDOBA COTORRA
## 10 23 350 23350 CÓRDOBA LA APARTADA
## mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1 Ordenanza 42 del 27 de Abril de 1923 MUNICIPIO 3136.8910 2023
## 2 Ordenanza 42 del 27 de Abril de 1923 MUNICIPIO 1964.8388 2023
## 3 Ordenanza 10 de Diciembre de 1969 MUNICIPIO 834.9270 2023
## 4 Ordenanza 28 de 1978 MUNICIPIO 420.0415 2023
## 5 1923 MUNICIPIO 290.8820 2023
## 6 1927 MUNICIPIO 323.9516 2023
## 7 Ordenanza 42 del 27 de Abril de 1923 MUNICIPIO 626.1762 2023
## 8 1923 MUNICIPIO 641.2572 2023
## 9 ORD 03 DE ABRIL 08 DE 1997 MUNICIPIO 88.1123 2023
## 10 ORD 07 DE MAYO 06 DE 1997 MUNICIPIO 287.0437 2023
## shape_Leng shape_Area geom area_m2
## 1 3.7233888 0.257383137 POINT (-75.9505 8.58464) 3146708009 [m^2]
## 2 2.4096791 0.161210107 POINT (-75.0487 8.265758) 1972548484 [m^2]
## 3 1.8888729 0.068472843 POINT (-75.43383 8.191389) 837983885 [m^2]
## 4 1.3044687 0.034465103 POINT (-76.23779 8.724093) 421208111 [m^2]
## 5 1.4926651 0.023889135 POINT (-75.83457 8.89102) 291824384 [m^2]
## 6 0.9175294 0.026625648 POINT (-75.65239 9.110253) 325055945 [m^2]
## 7 1.7241529 0.051466850 POINT (-75.35666 9.027499) 628471460 [m^2]
## 8 1.9769704 0.052669283 POINT (-75.60778 8.841785) 643481366 [m^2]
## 9 0.6153369 0.007240021 POINT (-75.77888 9.056992) 88402053 [m^2]
## 10 1.0054769 0.023535319 POINT (-75.28817 8.052216) 288129833 [m^2]
## area x y
## 1 3146.7100 [m^2] -75.95050 8.584640
## 2 1972.5500 [m^2] -75.04870 8.265758
## 3 837.9840 [m^2] -75.43383 8.191389
## 4 421.2080 [m^2] -76.23779 8.724093
## 5 291.8240 [m^2] -75.83457 8.891020
## 6 325.0560 [m^2] -75.65239 9.110253
## 7 628.4710 [m^2] -75.35666 9.027499
## 8 643.4810 [m^2] -75.60778 8.841785
## 9 88.4021 [m^2] -75.77888 9.056992
## 10 288.1300 [m^2] -75.28817 8.052216
#ahora se obtiene la DEM
elevation <- get_elev_raster(munic, z = 10)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
#Mosaicing & Projecting
#Note: Elevation units are in meters.
#objeto de elevación
elevation
## class : RasterLayer
## dimensions : 3568, 3090, 11025120 (nrow, ncol, ncell)
## resolution : 0.000682546, 0.000682546 (x, y)
## extent : -76.64062, -74.53156, 7.013738, 9.449062 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
## source : file2582500c1085.tif
## names : file2582500c1085
#escribir el DEM en su directorio local usando la función writeRaster
writeRaster(elevation, "~/Desktop/cordoba/elev_cordoba_z10.tif", overwrite=TRUE)
#convertir la elevación en un objeto ráster de Terra
(elevation2 <- terra::rast(elevation))
## class : SpatRaster
## dimensions : 3568, 3090, 1 (nrow, ncol, nlyr)
## resolution : 0.000682546, 0.000682546 (x, y)
## extent : -76.64062, -74.53156, 7.013738, 9.449062 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
## source : file2582500c1085.tif
## name : file2582500c1085
#Visualizar los datos de elevación.Primero, necesitaremos una paleta de colores
pal <- colorNumeric(c("cyan", "forestgreen","yellow","tan","orange", "brown"), values(elevation),
na.color = "transparent")
#Ahora, reduciremos la resolución espacial de los datos
(elevation3 <- terra::aggregate(elevation2, fact = 2))
## |---------|---------|---------|---------|=========================================
## class : SpatRaster
## dimensions : 1784, 1545, 1 (nrow, ncol, nlyr)
## resolution : 0.001365092, 0.001365092 (x, y)
## extent : -76.64062, -74.53156, 7.013738, 9.449062 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
## source(s) : memory
## name : file2582500c1085
## min value : -1238.5
## max value : 3711.0
#recortaremos el objeto de elevación3 para que se ajuste a los límites de cordoba
elevation4 <- terra::crop(elevation3, munic, mask=TRUE)
#Ahora es el momento de trazar - puede aparecer un mensaje de alerta sobre la capa que tiene datos inconsistentes...
leaflet(munic) %>% addTiles() %>% setView(-75.9, 8.7, 9) %>%
addPolygons(color = "white", weight = 1.0, smoothFactor = 0.5,
opacity = 0.25, fillOpacity = 0.15,
popup = paste("Municipio: ", munic$mpio_cnmbr, "<br>",
"Km2: ", munic$area, "<br>")) %>%
addLabelOnlyMarkers(data = centers,
lng = ~x, lat = ~y, label = ~mpio_cnmbr,
labelOptions = labelOptions(noHide = TRUE, direction = 'top', textOnly = TRUE, textsize = "10px")) %>%
addRasterImage(elevation4, colors = pal, opacity = 0.9) %>%
addLegend(pal = pal, values = values(elevation),
title = "Elevation data for Cordoba (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'
#DEM usando terra
(dem = terra::rast("~/Desktop/cordoba/elev_cordoba_z10.tif"))
## class : SpatRaster
## dimensions : 3568, 3090, 1 (nrow, ncol, nlyr)
## resolution : 0.000682546, 0.000682546 (x, y)
## extent : -76.64062, -74.53156, 7.013738, 9.449062 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
## source : elev_cordoba_z10.tif
## name : elev_cordoba_z10
## min value : -1241
## max value : 3724
#ahora vamos a reducir la resolución del DEM para evitar problemas de memoria:
dem2 = terra::aggregate(dem,2, "mean")
## |---------|---------|---------|---------|=========================================
#Ahora leeremos nuestros municipios usando la biblioteca
(munic <- sf::st_read("~/Desktop/cordoba/cordoba.gpkg"))
## Reading layer `cordoba' from data source
## `/Users/jsc_89/Desktop/cordoba/Cordoba.gpkg' using driver `GPKG'
## Simple feature collection with 30 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.51453 ymin: 7.347087 xmax: -74.78095 ymax: 9.447748
## Geodetic CRS: MAGNA-SIRGAS
## Simple feature collection with 30 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.51453 ymin: 7.347087 xmax: -74.78095 ymax: 9.447748
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 23 001 23001 CÓRDOBA MONTERÍA
## 2 23 068 23068 CÓRDOBA AYAPEL
## 3 23 079 23079 CÓRDOBA BUENAVISTA
## 4 23 090 23090 CÓRDOBA CANALETE
## 5 23 162 23162 CÓRDOBA CERETÉ
## 6 23 168 23168 CÓRDOBA CHIMÁ
## 7 23 182 23182 CÓRDOBA CHINÚ
## 8 23 189 23189 CÓRDOBA CIÉNAGA DE ORO
## 9 23 300 23300 CÓRDOBA COTORRA
## 10 23 350 23350 CÓRDOBA LA APARTADA
## mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1 Ordenanza 42 del 27 de Abril de 1923 MUNICIPIO 3136.8910 2023
## 2 Ordenanza 42 del 27 de Abril de 1923 MUNICIPIO 1964.8388 2023
## 3 Ordenanza 10 de Diciembre de 1969 MUNICIPIO 834.9270 2023
## 4 Ordenanza 28 de 1978 MUNICIPIO 420.0415 2023
## 5 1923 MUNICIPIO 290.8820 2023
## 6 1927 MUNICIPIO 323.9516 2023
## 7 Ordenanza 42 del 27 de Abril de 1923 MUNICIPIO 626.1762 2023
## 8 1923 MUNICIPIO 641.2572 2023
## 9 ORD 03 DE ABRIL 08 DE 1997 MUNICIPIO 88.1123 2023
## 10 ORD 07 DE MAYO 06 DE 1997 MUNICIPIO 287.0437 2023
## shape_Leng shape_Area geom
## 1 3.7233888 0.257383137 MULTIPOLYGON (((-76.12448 8...
## 2 2.4096791 0.161210107 MULTIPOLYGON (((-75.00565 8...
## 3 1.8888729 0.068472843 MULTIPOLYGON (((-75.35428 8...
## 4 1.3044687 0.034465103 MULTIPOLYGON (((-76.19951 8...
## 5 1.4926651 0.023889135 MULTIPOLYGON (((-75.75363 8...
## 6 0.9175294 0.026625648 MULTIPOLYGON (((-75.6027 9....
## 7 1.7241529 0.051466850 MULTIPOLYGON (((-75.39083 9...
## 8 1.9769704 0.052669283 MULTIPOLYGON (((-75.59987 9...
## 9 0.6153369 0.007240021 MULTIPOLYGON (((-75.77023 9...
## 10 1.0054769 0.023535319 MULTIPOLYGON (((-75.29656 8...
#Ahora, recortemos el objeto dem2 a los límites de nuestro departamento.
dem3 = terra::crop(dem2,munic, mask=TRUE)
#se Modifican los siguientes dos fragmentos de código para transformar los datos de entrada al https://rodolfofrancoweb.com/sig/proyecciones_y_sistemas_de_coordenadas/magna-sirgas-origen-nacional/
(dem_plane = project(dem3, "EPSG:9377"))
## class : SpatRaster
## dimensions : 1550, 1277, 1 (nrow, ncol, nlyr)
## resolution : 150.7599, 150.7599 (x, y)
## extent : 4611972, 4804493, 2370248, 2603926 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_cordoba_z10
## min value : -58.53565
## max value : 2688.35400
#necesitamos realizar una transformación similar para el objeto vectorial
(munic_plane = sf::st_transform(munic, "EPSG:9377"))
## Simple feature collection with 30 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4612240 ymin: 2370841 xmax: 4803924 ymax: 2603331
## Projected CRS: MAGNA-SIRGAS / Origen-Nacional
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 23 001 23001 CÓRDOBA MONTERÍA
## 2 23 068 23068 CÓRDOBA AYAPEL
## 3 23 079 23079 CÓRDOBA BUENAVISTA
## 4 23 090 23090 CÓRDOBA CANALETE
## 5 23 162 23162 CÓRDOBA CERETÉ
## 6 23 168 23168 CÓRDOBA CHIMÁ
## 7 23 182 23182 CÓRDOBA CHINÚ
## 8 23 189 23189 CÓRDOBA CIÉNAGA DE ORO
## 9 23 300 23300 CÓRDOBA COTORRA
## 10 23 350 23350 CÓRDOBA LA APARTADA
## mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1 Ordenanza 42 del 27 de Abril de 1923 MUNICIPIO 3136.8910 2023
## 2 Ordenanza 42 del 27 de Abril de 1923 MUNICIPIO 1964.8388 2023
## 3 Ordenanza 10 de Diciembre de 1969 MUNICIPIO 834.9270 2023
## 4 Ordenanza 28 de 1978 MUNICIPIO 420.0415 2023
## 5 1923 MUNICIPIO 290.8820 2023
## 6 1927 MUNICIPIO 323.9516 2023
## 7 Ordenanza 42 del 27 de Abril de 1923 MUNICIPIO 626.1762 2023
## 8 1923 MUNICIPIO 641.2572 2023
## 9 ORD 03 DE ABRIL 08 DE 1997 MUNICIPIO 88.1123 2023
## 10 ORD 07 DE MAYO 06 DE 1997 MUNICIPIO 287.0437 2023
## shape_Leng shape_Area geom
## 1 3.7233888 0.257383137 MULTIPOLYGON (((4656502 254...
## 2 2.4096791 0.161210107 MULTIPOLYGON (((4779315 250...
## 3 1.8888729 0.068472843 MULTIPOLYGON (((4740837 248...
## 4 1.3044687 0.034465103 MULTIPOLYGON (((4648131 253...
## 5 1.4926651 0.023889135 MULTIPOLYGON (((4697327 255...
## 6 0.9175294 0.026625648 MULTIPOLYGON (((4714138 257...
## 7 1.7241529 0.051466850 MULTIPOLYGON (((4737364 257...
## 8 1.9769704 0.052669283 MULTIPOLYGON (((4714292 255...
## 9 0.6153369 0.007240021 MULTIPOLYGON (((4695618 256...
## 10 1.0054769 0.023535319 MULTIPOLYGON (((4747044 245...
#Ahora, llamemos a la función SlpAsp
(slp_asp = MultiscaleDTM::SlpAsp(
dem_plane,
w = c(3, 3),
unit = "degrees",
method = "queen",
metrics = c("slope", "aspect"),
na.rm = TRUE,
include_scale = FALSE,
mask_aspect = TRUE))
## class : SpatRaster
## dimensions : 1550, 1277, 2 (nrow, ncol, nlyr)
## resolution : 150.7599, 150.7599 (x, y)
## extent : 4611972, 4804493, 2370248, 2603926 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0.0000, 1.234529e-04
## max values : 57.1558, 3.599994e+02
# Se divide el objeto slp_asp en dos partes. Capa 1
(slope = subset(slp_asp, 1))
## class : SpatRaster
## dimensions : 1550, 1277, 1 (nrow, ncol, nlyr)
## resolution : 150.7599, 150.7599 (x, y)
## extent : 4611972, 4804493, 2370248, 2603926 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 57.1558
#¿Cuál es la distribución de los valores de pendiente?
terra::hist(slope,
main = "Cordoba's slope",
xlab = "Slope (in degrees)")
## Warning: [hist] a sample of 51% of the cells was used (of which 44% was NA)
#Capa 2
(aspect =subset(slp_asp, 2))
## class : SpatRaster
## dimensions : 1550, 1277, 1 (nrow, ncol, nlyr)
## resolution : 150.7599, 150.7599 (x, y)
## extent : 4611972, 4804493, 2370248, 2603926 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 1.234529e-04
## max value : 3.599994e+02
#terra
terra::hist(aspect,
main = "Cordoba's aspect",
xlab = "Aspect (in degrees)")
## Warning: [hist] a sample of 51% of the cells was used (of which 44% was NA)
#convertimos los grados de pendiente en porcentaje de pendiente
(slope_perc = tan(slope*(pi/180))*100)
## class : SpatRaster
## dimensions : 1550, 1277, 1 (nrow, ncol, nlyr)
## resolution : 150.7599, 150.7599 (x, y)
## extent : 4611972, 4804493, 2370248, 2603926 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.000
## max value : 154.907
#terra
terra::hist(slope_perc,
main = "Cordoba's slope",
xlab = "Slope (in percentage)")
## Warning: [hist] a sample of 51% of the cells was used (of which 44% was NA)
# Reclassify the slope raster
#rc <- classify(slope_perc, c(0, 3, 7, 12, 25,50, 75), include.lowest=TRUE, brackets=TRUE)
m <- c(0, 3, 1,
3, 7, 2,
7, 12, 3,
12, 25, 4,
25, 50, 5,
50, 75, 6,
75, 160, 7)
m <- matrix(m, ncol=3, byrow = TRUE)
rc <- classify(slope_perc, m, right=TRUE)
#Ahora, calculemos estadísticas zonales
(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))
## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)
## | | | 0% | |== | 3% | |===== | 7% | |======= | 10% | |========= | 13% | |============ | 17% | |============== | 20% | |================ | 23% | |=================== | 27% | |===================== | 30% | |======================= | 33% | |========================== | 37% | |============================ | 40% | |============================== | 43% | |================================= | 47% | |=================================== | 50% | |===================================== | 53% | |======================================== | 57% | |========================================== | 60% | |============================================ | 63% | |=============================================== | 67% | |================================================= | 70% | |=================================================== | 73% | |====================================================== | 77% | |======================================================== | 80% | |========================================================== | 83% | |============================================================= | 87% | |=============================================================== | 90% | |================================================================= | 93% | |==================================================================== | 97% | |======================================================================| 100%
## [1] 2.5160389 1.2292731 1.4036171 5.2911134 0.9474660 0.5868814
## [7] 1.5853550 3.0762432 0.3953781 1.7273880 2.0587800 4.6855459
## [13] 1.4440339 10.1471510 3.8976312 3.0785139 2.0715215 3.6382058
## [19] 16.4352589 2.8243833 2.1665022 2.3823540 3.1995339 1.3972237
## [25] 3.4256206 11.3830862 1.5984520 14.8620186 2.7408032 5.6827564
#### Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)
#pendiente media de córdoba
hist(munic$mean_slope,
main = "Cordoba's mean slope",
xlab = "Slope (in percentage)")
#hist(munic$class, main = "Cordoba's reclassified slope", xlab = "Slope (as a category)")
#terra::hist(slope_perc)
#Transformemos la pendiente nuevamente a coordenadas geográficas
(rc.geo = project(rc, "EPSG:4326"))
## class : SpatRaster
## dimensions : 1556, 1291, 1 (nrow, ncol, nlyr)
## resolution : 0.001365132, 0.001365132 (x, y)
## extent : -76.53401, -74.77163, 7.33663, 9.460775 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 7
#Esta es la pendiente porcentual
(slope.geo = project(slope_perc, "EPSG:4326"))
## class : SpatRaster
## dimensions : 1556, 1291, 1 (nrow, ncol, nlyr)
## resolution : 0.001365132, 0.001365132 (x, y)
## extent : -76.53401, -74.77163, 7.33663, 9.460775 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.000
## max value : 133.009
#Primero, necesitaremos una paleta de colores para la pendiente
palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent")
#Ahora es el momento de trazar - pendiente en %
leaflet(munic) %>% addTiles() %>% setView(-75.9, 8.7, 9) %>%
addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.4, fillOpacity = 0.10,
popup = paste("Municipio: ", munic$mpio_cnmbr, "<br>",
"Slope class: ", munic$class, "<br>")) %>%
addRasterImage(slope.geo, colors = palredgreen, opacity = 0.8) %>%
addLegend(pal = palredgreen, values = values(slope.geo),
title = "Terrain slope in Cordoba (%)")
## 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'
#pendiente del terreno en grados y agregue etiquetas que muestren los nombres de los municipios junto con sus correspondientes valores de pendiente media.
## Si la pendiente está en porcentaje, conviértela a grados
slope.geo <- atan(slope.geo / 100) * (180 / pi)
# Extraer valores de pendiente media para cada municipio
valores_pendiente <- extract(slope.geo, munic, fun = mean, na.rm = TRUE)
## Warning: [extract] transforming vector data to the CRS of the raster
##Mensaje de advertencia: [extraer] transformando datos vectoriales al CRS del ráster
munic$pendiente_media <- valores_pendiente
#Ahora es el momento de trazar - pendiente en grados
leaflet(munic) %>%
addTiles() %>%
setView(-75.9, 8.7, zoom = 9) %>%
addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.4, fillOpacity = 0.10,
popup = paste("Municipio: ", munic$mpio_cnmbr, "<br>",
"Pendiente media: ", round(munic$pendiente_media, 2), "°<br>",
"Slope class: ", munic$class, "<br>")) %>%
addRasterImage(slope.geo, colors = palredgreen, opacity = 0.8) %>%
addLegend(pal = palredgreen, values = values(slope.geo),
title = "Pendiente del terreno (°)")
## 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'