rm(list=ls())
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.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(leaflet)
library(terra)
## terra 1.8.29
list.files()
## [1] "codigoR_elevacion.html" "codigoR_elevacion.nb.html"
## [3] "codigoR_elevacion.Rmd" "putumayo"
## [5] "putumayo.gpkg"
munic = sf::st_read("putumayo.gpkg")
## Reading layer `putumayo' from data source `D:\GB2\P4\R_elevacion\putumayo.gpkg' using driver `GPKG'
## Simple feature collection with 13 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -77.18681 ymin: -0.5622776 xmax: -73.84132 ymax: 1.467315
## Geodetic CRS: MAGNA-SIRGAS
(putumayo_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]]
munic$area_m2 = sf::st_area(munic)
munic$area = munic$area_m2/1000000
munic$area = signif(munic$area, digits = 6)
munic
## Simple feature collection with 13 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -77.18681 ymin: -0.5622776 xmax: -73.84132 ymax: 1.467315
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 86 001 86001 PUTUMAYO MOCOA
## 2 86 219 86219 PUTUMAYO COLÓN
## 3 86 320 86320 PUTUMAYO ORITO
## 4 86 568 86568 PUTUMAYO PUERTO ASÍS
## 5 86 569 86569 PUTUMAYO PUERTO CAICEDO
## 6 86 571 86571 PUTUMAYO PUERTO GUZMÁN
## 7 86 573 86573 PUTUMAYO PUERTO LEGUÍZAMO
## 8 86 749 86749 PUTUMAYO SIBUNDOY
## 9 86 755 86755 PUTUMAYO SAN FRANCISCO
## 10 86 757 86757 PUTUMAYO SAN MIGUEL
## mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1 1958 MUNICIPIO 1304.63788 2023
## 2 Decreto 2830 de Diciembre 2 de 1989 MUNICIPIO 64.27462 2023
## 3 Decreto 2891 de Diciembre 28 de 1978 MUNICIPIO 1939.35399 2023
## 4 Decreto 1951 de Octubre 24 de 1967 MUNICIPIO 2819.15478 2023
## 5 Ordenanza 12 de Noviembre 24 de 1992 MUNICIPIO 926.50246 2023
## 6 Ordenanza 13 de Noviembre 24 de 1992 MUNICIPIO 4576.59156 2023
## 7 Decreto 698 de Noviembre 13 de 1953 MUNICIPIO 10906.87681 2023
## 8 Decreto 1871 de Julio 1 de 1982 MUNICIPIO 97.73462 2023
## 9 Decreto 2830 de Diciembre 2 de 1989 MUNICIPIO 407.35674 2023
## 10 Ordenanza 45 de Abril 29 de 1994 MUNICIPIO 379.74213 2023
## shape_Leng shape_Area geom area_m2
## 1 2.4752281 0.105792947 MULTIPOLYGON (((-76.6705 1.... 1307757014 [m^2]
## 2 0.4548864 0.005209533 MULTIPOLYGON (((-76.96835 1... 64397806 [m^2]
## 3 2.1523001 0.157168461 MULTIPOLYGON (((-76.76607 0... 1943148182 [m^2]
## 4 3.6444021 0.228679089 MULTIPOLYGON (((-76.2263 0.... 2827365906 [m^2]
## 5 1.7199295 0.075142214 MULTIPOLYGON (((-76.46575 0... 929009823 [m^2]
## 6 4.7346323 0.371458749 MULTIPOLYGON (((-75.94095 1... 4592442615 [m^2]
## 7 7.5377279 0.885757929 MULTIPOLYGON (((-75.20032 0... 10951684883 [m^2]
## 8 0.5113819 0.007922269 MULTIPOLYGON (((-76.9043 1.... 97931160 [m^2]
## 9 1.1754950 0.033022563 MULTIPOLYGON (((-76.87345 1... 408222057 [m^2]
## 10 1.3275843 0.030777834 MULTIPOLYGON (((-76.99677 0... 380542476 [m^2]
## area
## 1 1307.7600 [m^2]
## 2 64.3978 [m^2]
## 3 1943.1500 [m^2]
## 4 2827.3700 [m^2]
## 5 929.0100 [m^2]
## 6 4592.4400 [m^2]
## 7 10951.7000 [m^2]
## 8 97.9312 [m^2]
## 9 408.2220 [m^2]
## 10 380.5420 [m^2]
4 Obtener centroide
(centers <- st_centroid(munic))
## Warning: st_centroid assumes attributes are constant over geometries
## Simple feature collection with 13 features and 13 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -77.00563 ymin: 0.05947393 xmax: -75.0648 ymax: 1.228746
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 86 001 86001 PUTUMAYO MOCOA
## 2 86 219 86219 PUTUMAYO COLÓN
## 3 86 320 86320 PUTUMAYO ORITO
## 4 86 568 86568 PUTUMAYO PUERTO ASÍS
## 5 86 569 86569 PUTUMAYO PUERTO CAICEDO
## 6 86 571 86571 PUTUMAYO PUERTO GUZMÁN
## 7 86 573 86573 PUTUMAYO PUERTO LEGUÍZAMO
## 8 86 749 86749 PUTUMAYO SIBUNDOY
## 9 86 755 86755 PUTUMAYO SAN FRANCISCO
## 10 86 757 86757 PUTUMAYO SAN MIGUEL
## mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1 1958 MUNICIPIO 1304.63788 2023
## 2 Decreto 2830 de Diciembre 2 de 1989 MUNICIPIO 64.27462 2023
## 3 Decreto 2891 de Diciembre 28 de 1978 MUNICIPIO 1939.35399 2023
## 4 Decreto 1951 de Octubre 24 de 1967 MUNICIPIO 2819.15478 2023
## 5 Ordenanza 12 de Noviembre 24 de 1992 MUNICIPIO 926.50246 2023
## 6 Ordenanza 13 de Noviembre 24 de 1992 MUNICIPIO 4576.59156 2023
## 7 Decreto 698 de Noviembre 13 de 1953 MUNICIPIO 10906.87681 2023
## 8 Decreto 1871 de Julio 1 de 1982 MUNICIPIO 97.73462 2023
## 9 Decreto 2830 de Diciembre 2 de 1989 MUNICIPIO 407.35674 2023
## 10 Ordenanza 45 de Abril 29 de 1994 MUNICIPIO 379.74213 2023
## shape_Leng shape_Area geom area_m2
## 1 2.4752281 0.105792947 POINT (-76.68436 1.228571) 1307757014 [m^2]
## 2 0.4548864 0.005209533 POINT (-76.96673 1.223217) 64397806 [m^2]
## 3 2.1523001 0.157168461 POINT (-76.94314 0.6737547) 1943148182 [m^2]
## 4 3.6444021 0.228679089 POINT (-76.3282 0.4786145) 2827365906 [m^2]
## 5 1.7199295 0.075142214 POINT (-76.481 0.7250721) 929009823 [m^2]
## 6 4.7346323 0.371458749 POINT (-75.88022 0.7464368) 4592442615 [m^2]
## 7 7.5377279 0.885757929 POINT (-75.0648 0.05947393) 10951684883 [m^2]
## 8 0.5113819 0.007922269 POINT (-76.91308 1.228746) 97931160 [m^2]
## 9 1.1754950 0.033022563 POINT (-76.84727 1.13607) 408222057 [m^2]
## 10 1.3275843 0.030777834 POINT (-76.85776 0.3058743) 380542476 [m^2]
## area
## 1 1307.7600 [m^2]
## 2 64.3978 [m^2]
## 3 1943.1500 [m^2]
## 4 2827.3700 [m^2]
## 5 929.0100 [m^2]
## 6 4592.4400 [m^2]
## 7 10951.7000 [m^2]
## 8 97.9312 [m^2]
## 9 408.2220 [m^2]
## 10 380.5420 [m^2]
centers$x = st_coordinates(centers)[,1]
centers$y = st_coordinates(centers)[,2]
centers
## Simple feature collection with 13 features and 15 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -77.00563 ymin: 0.05947393 xmax: -75.0648 ymax: 1.228746
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 86 001 86001 PUTUMAYO MOCOA
## 2 86 219 86219 PUTUMAYO COLÓN
## 3 86 320 86320 PUTUMAYO ORITO
## 4 86 568 86568 PUTUMAYO PUERTO ASÍS
## 5 86 569 86569 PUTUMAYO PUERTO CAICEDO
## 6 86 571 86571 PUTUMAYO PUERTO GUZMÁN
## 7 86 573 86573 PUTUMAYO PUERTO LEGUÍZAMO
## 8 86 749 86749 PUTUMAYO SIBUNDOY
## 9 86 755 86755 PUTUMAYO SAN FRANCISCO
## 10 86 757 86757 PUTUMAYO SAN MIGUEL
## mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1 1958 MUNICIPIO 1304.63788 2023
## 2 Decreto 2830 de Diciembre 2 de 1989 MUNICIPIO 64.27462 2023
## 3 Decreto 2891 de Diciembre 28 de 1978 MUNICIPIO 1939.35399 2023
## 4 Decreto 1951 de Octubre 24 de 1967 MUNICIPIO 2819.15478 2023
## 5 Ordenanza 12 de Noviembre 24 de 1992 MUNICIPIO 926.50246 2023
## 6 Ordenanza 13 de Noviembre 24 de 1992 MUNICIPIO 4576.59156 2023
## 7 Decreto 698 de Noviembre 13 de 1953 MUNICIPIO 10906.87681 2023
## 8 Decreto 1871 de Julio 1 de 1982 MUNICIPIO 97.73462 2023
## 9 Decreto 2830 de Diciembre 2 de 1989 MUNICIPIO 407.35674 2023
## 10 Ordenanza 45 de Abril 29 de 1994 MUNICIPIO 379.74213 2023
## shape_Leng shape_Area geom area_m2
## 1 2.4752281 0.105792947 POINT (-76.68436 1.228571) 1307757014 [m^2]
## 2 0.4548864 0.005209533 POINT (-76.96673 1.223217) 64397806 [m^2]
## 3 2.1523001 0.157168461 POINT (-76.94314 0.6737547) 1943148182 [m^2]
## 4 3.6444021 0.228679089 POINT (-76.3282 0.4786145) 2827365906 [m^2]
## 5 1.7199295 0.075142214 POINT (-76.481 0.7250721) 929009823 [m^2]
## 6 4.7346323 0.371458749 POINT (-75.88022 0.7464368) 4592442615 [m^2]
## 7 7.5377279 0.885757929 POINT (-75.0648 0.05947393) 10951684883 [m^2]
## 8 0.5113819 0.007922269 POINT (-76.91308 1.228746) 97931160 [m^2]
## 9 1.1754950 0.033022563 POINT (-76.84727 1.13607) 408222057 [m^2]
## 10 1.3275843 0.030777834 POINT (-76.85776 0.3058743) 380542476 [m^2]
## area x y
## 1 1307.7600 [m^2] -76.68436 1.22857059
## 2 64.3978 [m^2] -76.96673 1.22321748
## 3 1943.1500 [m^2] -76.94314 0.67375467
## 4 2827.3700 [m^2] -76.32820 0.47861446
## 5 929.0100 [m^2] -76.48100 0.72507212
## 6 4592.4400 [m^2] -75.88022 0.74643679
## 7 10951.7000 [m^2] -75.06480 0.05947393
## 8 97.9312 [m^2] -76.91308 1.22874637
## 9 408.2220 [m^2] -76.84727 1.13607034
## 10 380.5420 [m^2] -76.85776 0.30587427
elevation <- get_elev_raster(munic, z = 10)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
elevation
## class : RasterLayer
## dimensions : 3584, 5120, 18350080 (nrow, ncol, ncell)
## resolution : 0.0006866186, 0.0006866186 (x, y)
## extent : -77.34375, -73.82826, -0.7033042, 1.757537 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
## source : file292c492e3f85.tif
## names : file292c492e3f85
writeRaster(elevation, "putumayo/elev_putumayo_z10.tif", overwrite = TRUE)
(elevation2 <- terra::rast(elevation))
## class : SpatRaster
## dimensions : 3584, 5120, 1 (nrow, ncol, nlyr)
## resolution : 0.0006866186, 0.0006866186 (x, y)
## extent : -77.34375, -73.82826, -0.7033042, 1.757537 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
## source : file292c492e3f85.tif
## name : file292c492e3f85
5 visualización del mapa
pal <- colorNumeric(c("cyan", "forestgreen","yellow","tan","orange", "brown"), values(elevation),
na.color = "transparent")
(elevation3 <- terra::aggregate(elevation2, fact = 2))
##
|---------|---------|---------|---------|
=========================================
## class : SpatRaster
## dimensions : 1792, 2560, 1 (nrow, ncol, nlyr)
## resolution : 0.001373237, 0.001373237 (x, y)
## extent : -77.34375, -73.82826, -0.7033042, 1.757537 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
## source(s) : memory
## name : file292c492e3f85
## min value : -146.00
## max value : 4170.75
elevation4 <- terra::crop(elevation3, munic, mask=TRUE)
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 = "Elevación Putumayo (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'