ElevaciónPutumayo

library(rasterVis)
library(raster)
library(rgl)
library(rgdal)
library(elevatr)
library(sf)
library(tidyverse)
library(mapview)

putumayo <- st_read("E:/Descargas en el disco duro/4to Semestre/Geomatica/Putumayo/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")

Reading layer `MGN_MPIO_POLITICO' from data source 
  `E:\Descargas en el disco duro\4to Semestre\Geomatica\Putumayo\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 13 features and 9 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -77.18681 ymin: -0.5622776 xmax: -73.84132 ymax: 1.467315
Geodetic CRS:  WGS 84

head(putumayo)

Simple feature collection with 6 features and 9 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -77.18681 ymin: 0.2391261 xmax: -76.41003 ymax: 1.467315
Geodetic CRS:  WGS 84
  DPTO_CCDGO MPIO_CCDGO             MPIO_CNMBR
1         86      86001                  MOCOA
2         86      86219                  COLÓN
3         86      86320                  ORITO
4         86      86749               SIBUNDOY
5         86      86755          SAN FRANCISCO
6         86      86757 SAN MIGUEL (La Dorada)
                            MPIO_CRSLC MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng
1                                 1958 1304.63835      2017   PUTUMAYO  2.4752281
2  Decreto 2830 de Diciembre 2 de 1989   64.27462      2017   PUTUMAYO  0.4548864
3 Decreto 2891 de Diciembre 28 de 1978 1939.39517      2017   PUTUMAYO  2.1520883
4      Decreto 1871 de Julio 1 de 1982   97.73462      2017   PUTUMAYO  0.5113819
5  Decreto 2830 de Diciembre 2 de 1989  407.35674      2017   PUTUMAYO  1.1754950
6     Ordenanza 45 de Abril 29 de 1994  379.74249      2017   PUTUMAYO  1.3275843
   Shape_Area                       geometry
1 0.105792947 POLYGON ((-76.6705 1.467315...
2 0.005209533 POLYGON ((-76.96835 1.28631...
3 0.157171417 POLYGON ((-77.07275 0.94231...
4 0.007922269 POLYGON ((-76.9043 1.299191...
5 0.033022563 POLYGON ((-76.87345 1.28986...
6 0.030777834 POLYGON ((-76.99677 0.37418...

elevation <- get_elev_raster(putumayo, z = 10)

 Accessing raster elevation [-------------------------]   0%
 Accessing raster elevation [-------------------------]   1%
 Accessing raster elevation [>------------------------]   2%
 Accessing raster elevation [>------------------------]   3%
 Accessing raster elevation [>------------------------]   5%
 Accessing raster elevation [>------------------------]   6%
 Accessing raster elevation [=>-----------------------]   7%
 Accessing raster elevation [=>-----------------------]   8%
 Accessing raster elevation [=>-----------------------]   9%
 Accessing raster elevation [==>----------------------]  10%
 Accessing raster elevation [==>----------------------]  11%
 Accessing raster elevation [==>----------------------]  12%
 Accessing raster elevation [==>----------------------]  14%
 Accessing raster elevation [===>---------------------]  15%
 Accessing raster elevation [===>---------------------]  16%
 Accessing raster elevation [===>---------------------]  17%
 Accessing raster elevation [====>--------------------]  18%
 Accessing raster elevation [====>--------------------]  19%
 Accessing raster elevation [====>--------------------]  20%
 Accessing raster elevation [====>--------------------]  22%
 Accessing raster elevation [=====>-------------------]  23%
 Accessing raster elevation [=====>-------------------]  24%
 Accessing raster elevation [=====>-------------------]  25%
 Accessing raster elevation [======>------------------]  26%
 Accessing raster elevation [======>------------------]  27%
 Accessing raster elevation [======>------------------]  28%
 Accessing raster elevation [======>------------------]  30%
 Accessing raster elevation [=======>-----------------]  31%
 Accessing raster elevation [=======>-----------------]  32%
 Accessing raster elevation [=======>-----------------]  33%
 Accessing raster elevation [========>----------------]  34%
 Accessing raster elevation [========>----------------]  35%
 Accessing raster elevation [========>----------------]  36%
 Accessing raster elevation [========>----------------]  38%
 Accessing raster elevation [=========>---------------]  39%
 Accessing raster elevation [=========>---------------]  40%
 Accessing raster elevation [=========>---------------]  41%
 Accessing raster elevation [==========>--------------]  42%
 Accessing raster elevation [==========>--------------]  43%
 Accessing raster elevation [==========>--------------]  44%
 Accessing raster elevation [==========>--------------]  45%
 Accessing raster elevation [===========>-------------]  47%
 Accessing raster elevation [===========>-------------]  48%
 Accessing raster elevation [===========>-------------]  49%
 Accessing raster elevation [===========>-------------]  50%
 Accessing raster elevation [============>------------]  51%
 Accessing raster elevation [============>------------]  52%
 Accessing raster elevation [============>------------]  53%
 Accessing raster elevation [=============>-----------]  55%
 Accessing raster elevation [=============>-----------]  56%
 Accessing raster elevation [=============>-----------]  57%
 Accessing raster elevation [=============>-----------]  58%
 Accessing raster elevation [==============>----------]  59%
 Accessing raster elevation [==============>----------]  60%
 Accessing raster elevation [==============>----------]  61%
 Accessing raster elevation [===============>---------]  62%
 Accessing raster elevation [===============>---------]  64%
 Accessing raster elevation [===============>---------]  65%
 Accessing raster elevation [===============>---------]  66%
 Accessing raster elevation [================>--------]  67%
 Accessing raster elevation [================>--------]  68%
 Accessing raster elevation [================>--------]  69%
 Accessing raster elevation [=================>-------]  70%
 Accessing raster elevation [=================>-------]  72%
 Accessing raster elevation [=================>-------]  73%
 Accessing raster elevation [=================>-------]  74%
 Accessing raster elevation [==================>------]  75%
 Accessing raster elevation [==================>------]  76%
 Accessing raster elevation [==================>------]  77%
 Accessing raster elevation [===================>-----]  78%
 Accessing raster elevation [===================>-----]  80%
 Accessing raster elevation [===================>-----]  81%
 Accessing raster elevation [===================>-----]  82%
 Accessing raster elevation [====================>----]  83%
 Accessing raster elevation [====================>----]  84%
 Accessing raster elevation [====================>----]  85%
 Accessing raster elevation [=====================>---]  86%
 Accessing raster elevation [=====================>---]  88%
 Accessing raster elevation [=====================>---]  89%
 Accessing raster elevation [=====================>---]  90%
 Accessing raster elevation [======================>--]  91%
 Accessing raster elevation [======================>--]  92%
 Accessing raster elevation [======================>--]  93%
 Accessing raster elevation [=======================>-]  94%
 Accessing raster elevation [=======================>-]  95%
 Accessing raster elevation [=======================>-]  97%
 Accessing raster elevation [=======================>-]  98%
 Accessing raster elevation [========================>]  99%
 Accessing raster elevation [=========================] 100%

Mosaicing & Projecting
Note: Elevation units are in meters.

elevation

class      : RasterLayer 
dimensions : 4096, 5633, 23072768  (nrow, ncol, ncell)
resolution : 0.0006865141, 0.0006865141  (x, y)
extent     : -77.34375, -73.47662, -1.054425, 1.757537  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : file46044f64ddc.tif 
names      : file46044f64ddc 
values     : -32768, 32767  (min, max)

munic_putu <- as(putumayo, "Spatial")

plot(elevation, main="DEM putumayo [Metros]")
plot(munic_putu,col="NA",border="black", add=TRUE)

text(coordinates(munic_putu), labels=as.character(putumayo$MPIO_CNMBR), 
     col="black", cex=0.30)

elev_crop = crop(elevation, munic_putu)
elev_crop = mask(elev_crop, munic_putu)
plot(elev_crop, main="DEM Putumayo Recortado")
plot(munic_putu, add=TRUE)

text(coordinates(munic_putu), labels=as.character(munic_putu$MPIO_CNMBR), cex=0.3)

elev_crop

class      : RasterLayer 
dimensions : 2956, 4873, 14404588  (nrow, ncol, ncell)
resolution : 0.0006865141, 0.0006865141  (x, y)
extent     : -77.18654, -73.84115, -0.5621945, 1.467141  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : file46044f64ddc 
values     : -605, 3920  (min, max)

crs(elev_crop)

CRS arguments: +proj=longlat +datum=WGS84 +no_defs 

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 : 1247, 2068, 2578796  (nrow, ncol, ncell)
resolution : 180, 180  (x, y)
extent     : 256633.6, 628873.6, -62177.54, 162282.5  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
source     : memory
names      : file46044f64ddc 
values     : -188.6136, 3908.965  (min, max)

(rep_munic = spTransform(munic_putu,spatialref))

class       : SpatialPolygonsDataFrame 
features    : 13 
extent      : 256618.4, 628938.1, -62151.84, 162252.4  (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 
min values  :         86,      86001,       COLÓN,                             1958,   64.27461649,      2017,   PUTUMAYO, 0.454886403409, 0.00520953330412 
max values  :         86,      86885, VILLAGARZÓN, Ordenanza 45 de Abril 29 de 1994, 10906.8837683,      2017,   PUTUMAYO,  7.53772787618,   0.885757928701 

plot(rep_elev, main="Reprojected digital elevation model")
plot(rep_munic, lwd=0.5, add=TRUE)

text(coordinates(rep_munic), labels=as.character(rep_munic$MPIO_CNMBR), cex=0.4
     )

writeRaster(rep_elev, filename="E:/Descargas en el disco duro/4to Semestre/Geomatica/DEM/putumayoDEM.tif", datatype='INT4S', overwrite=TRUE)

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('Promedio', 'Min', 'Max', 'DV')
Valores <- c(promedio, minimo, maximo, desviacion)

(df_estadisticas <- data.frame(Metricas, Valores))

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 de los municipios de Putumayo [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.3)

plot(slope,main="Pendiente municipal de Putumayo [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.4)

plot(aspect,main="Aspecto municipal de Putumayo [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.4)

plot(hill,
        col=grey(1:100/100),  
        legend=FALSE,         
        main="DEM municipal de Putumayo [Metros]",
        axes=FALSE)           

plot(rep_elev, 
        axes=FALSE,
        col=terrain.colors(12, alpha=0.6), add=TRUE) 

         
plot(rep_munic, add=TRUE, lwd=0.5)
text(coordinates(rep_munic), labels=as.character(rep_munic$MPIO_CNMBR), cex=0.5)

lista <- list("MOCOA", "SAN FRANCISCO","SIBUNDOY", "COLÓN", "VILLAGARZÓN", "SANTIAGO")
(algunosmunicipios <- putumayo %>% filter(MPIO_CNMBR %in% unlist(lista) ))

Simple feature collection with 6 features and 9 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -77.10145 ymin: 0.7693426 xmax: -76.41003 ymax: 1.467315
Geodetic CRS:  WGS 84
  DPTO_CCDGO MPIO_CCDGO    MPIO_CNMBR                          MPIO_CRSLC
1         86      86001         MOCOA                                1958
2         86      86219         COLÓN Decreto 2830 de Diciembre 2 de 1989
3         86      86749      SIBUNDOY     Decreto 1871 de Julio 1 de 1982
4         86      86755 SAN FRANCISCO Decreto 2830 de Diciembre 2 de 1989
5         86      86760      SANTIAGO Decreto 2830 de Diciembre 2 de 1989
6         86      86885   VILLAGARZÓN     Decreto 574 de Marzo 14 de 1977
  MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng  Shape_Area
1 1304.63835      2017   PUTUMAYO  2.4752281 0.105792947
2   64.27462      2017   PUTUMAYO  0.4548864 0.005209533
3   97.73462      2017   PUTUMAYO  0.5113819 0.007922269
4  407.35674      2017   PUTUMAYO  1.1754950 0.033022563
5  341.92382      2017   PUTUMAYO  0.9881418 0.027710081
6 1396.96644      2017   PUTUMAYO  1.9280172 0.113257175
                        geometry
1 POLYGON ((-76.6705 1.467315...
2 POLYGON ((-76.96835 1.28631...
3 POLYGON ((-76.9043 1.299191...
4 POLYGON ((-76.87345 1.28986...
5 POLYGON ((-77.0344 1.199208...
6 POLYGON ((-76.63426 1.06411...

algunosmunicipios <- as(algunosmunicipios, "Spatial")

elev_cropM = crop(elevation, algunosmunicipios)
elev_cropM = mask(elev_cropM, algunosmunicipios)
plot(elev_cropM, main=" DEM Zona Noroccidente de Putumayo ")
plot(algunosmunicipios, add=TRUE)

text(coordinates(algunosmunicipios), labels=as.character(algunosmunicipios$MPIO_CNMBR), cex=0.5)

slopeM = terrain(elev_cropM, opt='slope', unit='degrees')
aspectM = terrain(elev_cropM, opt='aspect',unit='degrees')
hillM = hillShade(slopeM, aspectM, 40, 315)

plot(slopeM,main="Pendiente Noroccidental de Putumayo [grados]", col=topo.colors(6,alpha=0.6))
plot(algunosmunicipios, add=TRUE, lwd=0.5)

text(coordinates(algunosmunicipios), labels=as.character(algunosmunicipios$MPIO_CNMBR), cex=0.6)

plot(aspectM,main="Aspecto Noroccidental de Putumayo", col=rainbow(10,alpha=0.7))
plot(algunosmunicipios, add=TRUE, lwd=0.5)

text(coordinates(algunosmunicipios), labels=as.character(algunosmunicipios$MPIO_CNMBR), cex=0.6)

Todos los codigos usados en este cuarderno de R fueron tomados del cuaderno del profesor (Lizarazo, 2021a).