Geomatics
2025-07-23
Instalar y cargar las librerias
library(MultiscaleDTM)## Cargando paquete requerido: terra## terra 1.8.50library(exactextractr)
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.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUElibrary(leaflet)
library(terra)
library(dplyr)##
## Adjuntando el paquete: 'dplyr'## The following objects are masked from 'package:terra':
##
## intersect, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unioncargar archivo de polígonos
munic = st_read("Division_muncipios_completos.shp")## Reading layer `Division_muncipios_completos' from data source
## `/home/mateo/Downloads/Geom_Liza/Division_muncipios_completos.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1122 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -81.73562 ymin: -4.225938 xmax: -66.85689 ymax: 13.39473
## Geodetic CRS: WGS 84GUAINIA_SHP = munic %>%
filter(Depto == "Guainía")3.Introduction to MultiscaleDTM
Cargar el archivo .tiff del GUAINIA siguiendo la sesión anterior usando terra
dem = terra::rast("Guainia/elev_guainia_z10.tif")
dem## class : SpatRaster
## dimensions : 4605, 6147, 1 (nrow, ncol, nlyr)
## resolution : 0.0006863529, 0.0006863529 (x, y)
## extent : -71.01562, -66.79661, 1.054288, 4.214943 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : elev_guainia_z10.tif
## name : elev_guainia_z10
## min value : -799
## max value : 1153agregar por la media de los pixeles los modelos de elevación reduciendo la resolución.
dem2 = terra::aggregate(dem, 2, "mean")## |---------|---------|---------|---------|========================================= Cortar mi modelo de elevación a mi capa de polígonos en el TOLIMA para que me quede solo la información que esta dentro de esa geometría. Y lo que esté pro fuera lo deja como NA.
dem3 = terra::crop(dem2,GUAINIA_SHP, mask=TRUE)5. Transforming coordinates
El EPSG es un identificador único para sistemas de coordenadas y proyecciones geográficas. Cada código representa un Sistema de Referencia de Coordenadas (CRS) con el correspondiente datum, proyección y unidad de medida. Es muy importante para que las operaciones entre las capas se de correctamente. Si están en diferente sistema de coordenadas no se pueden operar.
Entonces, reproyectamos al sistema de coordenadas planas con origen único nacional
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## dimensions : 2098, 2989, 1 (nrow, ncol, nlyr)
## resolution : 152.66, 152.66 (x, y)
## extent : 5228299, 5684600, 1686946, 2007227 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_guainia_z10
## min value : -65.87392
## max value : 960.34509(munic_plane = sf::st_transform(GUAINIA_SHP, "EPSG:9377"))## Simple feature collection with 8 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 5228672 ymin: 1687512 xmax: 5684465 ymax: 2006550
## Projected CRS: MAGNA-SIRGAS / Origen-Nacional
## MpCodigo MpNombre Depto Cod_Mun geometry
## 1 94001 Inírida Guainía 94001 MULTIPOLYGON (((5587784 200...
## 2 94343 Barrancominas Guainía 94343 MULTIPOLYGON (((5451281 197...
## 3 94883 San Felipe Guainía 94883 MULTIPOLYGON (((5628860 183...
## 4 94884 Puerto Colombia Guainía 94884 MULTIPOLYGON (((5606955 191...
## 5 94885 La Guadalupe Guainía 94885 MULTIPOLYGON (((5667208 174...
## 6 94886 Cacahual Guainía 94886 MULTIPOLYGON (((5611930 197...
## 7 94887 Paná-Paná Guainía 94887 MULTIPOLYGON (((5373852 174...
## 8 94888 Morichal Guainía 94888 MULTIPOLYGON (((5392692 187...6. Computing terrain atributes
- Multiscale Slope and Aspect
¿Qué es w? w es el tamaño de la ventana de análisis, en este caso es en filas y columnas (3,3) es decir que la pendiente y el aspecto se calculan usando una matriz de vecinos de 3x3 pixeles
¿Qué es method? method es la forma en como se eleigen los pixeles vecinos para el cálculo, en el caso de la opción “queen” indica que tomará los 8 vecinos cercanos al rededor de la celda. Puede ser cambiado a “rook” que serian los vecinos en cruz
Calculamos la pendiente y el aspecto a partir del DEM
(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 : 2098, 2989, 2 (nrow, ncol, nlyr)
## resolution : 152.66, 152.66 (x, y)
## extent : 5228299, 5684600, 1686946, 2007227 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 4.454152e-04, 0
## max values : 5.067653e+01, 360Tomar solo la capa de pendiente
#capa1
(slope = subset(slp_asp, 1))## class : SpatRaster
## dimensions : 2098, 2989, 1 (nrow, ncol, nlyr)
## resolution : 152.66, 152.66 (x, y)
## extent : 5228299, 5684600, 1686946, 2007227 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 4.454152e-04
## max value : 5.067653e+01terra::hist(slope,
main = "GUAINIA's slope",
xlab = "Slope (in degrees)")## Warning: [hist] a sample of 16% of the cells was used (of which 51% was NA)
El histograma muestra la mayoria de los valores bajos cerca a 0, esto
indica que GUAINIA tiene pendientes muy bajas, la mayor cantidad de su
terreno es plano.
Tomar solo la capa aspecto
#capa2
(aspect =subset(slp_asp, 2))## class : SpatRaster
## dimensions : 2098, 2989, 1 (nrow, ncol, nlyr)
## resolution : 152.66, 152.66 (x, y)
## extent : 5228299, 5684600, 1686946, 2007227 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 0
## max value : 360terra::hist(aspect,
main = "GUAINIA's aspect",
xlab = "Aspect (in degrees)")## Warning: [hist] a sample of 16% of the cells was used (of which 51% was NA)
Pendiente en porcentajes
slope_perc = tan(slope*(pi/180))*100
slope_perc## class : SpatRaster
## dimensions : 2098, 2989, 1 (nrow, ncol, nlyr)
## resolution : 152.66, 152.66 (x, y)
## extent : 5228299, 5684600, 1686946, 2007227 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 7.773961e-04
## max value : 1.220741e+02terra::hist(slope_perc,
main = "GUAINIA's slope",
xlab = "Slope (in percentage)")## Warning: [hist] a sample of 16% of the cells was used (of which 51% was NA)
7. Calcular estadísticas zonales
Creamos la clasificación del IGAC para las pendientes. y clasificamos la pendiente en porcentaje con el valor superior incluido en cada categoria
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 calculamos la media de cada grupo clasificado
(GUAINIA_SHP$mean_slope <- exactextractr::exact_extract(slope_perc, GUAINIA_SHP, 'mean'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |========= | 12% | |================== | 25% | |========================== | 38% | |=================================== | 50% | |============================================ | 62% | |==================================================== | 75% | |============================================================= | 88% | |======================================================================| 100%## [1] 1.851265 2.484470 1.686722 2.279491 1.130117 1.498820 2.471356 2.574098Histograma para ver la distribución de esas medias
hist(GUAINIA_SHP$mean_slope,
main = "GUAINIA's mean slope",
xlab = "Slope (in percentage)", breaks = 18)
calcular la moda por municipio
GUAINIA_SHP$class <- exactextractr::exact_extract(rc, GUAINIA_SHP, 'mode')## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |========= | 12% | |================== | 25% | |========================== | 38% | |=================================== | 50% | |============================================ | 62% | |==================================================== | 75% | |============================================================= | 88% | |======================================================================| 100%GUAINIA_SHP$class## [1] 1 1 1 1 1 1 1 1distribución de las modas
hist(GUAINIA_SHP$class,
main = "GUAINIA's reclassified slope",
xlab = "Slope (as a category)", breaks = 12)
pendiente en coordenadas geográficas
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## dimensions : 2115, 2991, 1 (nrow, ncol, nlyr)
## resolution : 0.001372744, 0.001372744 (x, y)
## extent : -70.94753, -66.84166, 1.159429, 4.062781 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 1
## max value : 7(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## dimensions : 2115, 2991, 1 (nrow, ncol, nlyr)
## resolution : 0.001372744, 0.001372744 (x, y)
## extent : -70.94753, -66.84166, 1.159429, 4.062781 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.0085014
## max value : 103.55226148. Plotting terrain slope
Ajustar la paleta de colores con los valores de la pendiente reproyectada y que los NA sean trasnparentes
palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2","darkred"), values(slope.geo), na.color = "transparent", reverse = T)reproyectar capa de municipios para que se puedan visualizar en el mismo sistema de coordenadas
GUAINIA_SHP <- st_transform(GUAINIA_SHP, crs = 4326)graficar los polígonos y la pendiente con una letenda de los valores
leaflet(GUAINIA_SHP) %>%
addTiles() %>%
setView(-68.52293954604566, 2.7601747799884104, 8) %>%
addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.4, fillOpacity = 0.10,
popup = paste("Municipio: ", GUAINIA_SHP$MpNombre, "<br>",
"Slope class: ", GUAINIA_SHP$class, "<br>")) %>%
addRasterImage(slope.geo, opacity = 0.8) %>%
addLegend(pal = palredgreen, values = values(slope.geo),
title = "Terrain slope in GUAINIA (%)")