Quinto cuaderno: DEM de Valle del Cauca

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

(munic <-  st_read("C:/Users/LUISA CARRION/Downloads/76_VALLE_DEL_CAUCA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

elevation <- get_elev_raster(munic, z = 10)

if (require(rgdal)) {
 rf <- writeRaster(elevation, filename=file.path("/Users/LUISA CARRION/Downloads/76_VALLE_DEL_CAUCA/ADMINISTRATIVO/valle_dem_z10.tif"), format="GTiff", overwrite=TRUE)
}

elevation

class      : RasterLayer 
dimensions : 4095, 3592, 14709240  (nrow, ncol, ncell)
resolution : 0.0006852919, 0.0006852919  (x, y)
extent     : -77.69531, -75.23374, 2.459737, 5.266008  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : C:/Users/LUISA CARRION/AppData/Local/Temp/RtmpQBuZ5Y/file5fc4cb2b22.tif 
names      : file5fc4cb2b22 
values     : -32768, 32767  (min, max)

munic_sp <- as(munic, "Spatial")

plot(elevation, main="This the downloaded DEM [meters]")
plot(munic_sp,col="NA",border="black", add=TRUE)

text(coordinates(munic_sp), labels=as.character(munic$MPIO_CNMBR), 
     col="black", cex=0.20)

elev_crop = crop(elevation, munic_sp)
plot(elev_crop, main="Cropped digital elevation model")
plot(munic_sp, add=TRUE)

text(coordinates(munic_sp), labels=as.character(munic_sp$MPIO_CNMBR), cex=0.2)

elev_crop

class      : RasterLayer 
dimensions : 2854, 2689, 7674406  (nrow, ncol, ncell)
resolution : 0.0006852919, 0.0006852919  (x, y)
extent     : -77.55003, -75.70728, 3.091577, 5.0474  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : file29182d35f46 
values     : -1393, 4574  (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 : 1204, 1139, 1371356  (nrow, ncol, ncell)
resolution : 180, 180  (x, y)
extent     : 216562, 421582, 341737.7, 558457.7  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
source     : memory
names      : file5fc4cb2b22 
values     : -1392.295, 4552.94  (min, max)

(rep_munic = spTransform(munic_sp,spatialref))

class       : SpatialPolygonsDataFrame 
features    : 42 
extent      : 216628.4, 421558.5, 341813.7, 558001.9  (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  :         76,      76001,     ALCALÁ,                          1536,   41.86090736,      2017, VALLE DEL CAUCA, 0.453826056161, 0.00340901158098 
max values  :         76,      76895,     ZARZAL, Ordenanza 9 de Diciembre 1954, 6292.50083741,      2017, VALLE DEL CAUCA,  6.59527269127,   0.510778519244 

plot(rep_elev, main="DEM de Valle del Cauca")
plot(rep_munic, lwd=0.5, add=TRUE)

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

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))

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 for Valle municipalities [meters]", col=heat.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 los municipios de Valle [grados]", col=cm.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="Aspect for several municipalities [degrees]", 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 Valle municipalities [m]",
        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.2)