My fifth R Notebook
2.Introduction to elevatr
Permite traer el DEM, y tambien permite traer zonas con datos de terreno
library(rasterVis)## Loading required package: raster## Loading required package: sp## Loading required package: terra## terra version 1.3.4## Loading required package: lattice## Loading required package: latticeExtralibrary(raster)
library(rgl)
library(rgdal)## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/morea/OneDrive - Universidad Nacional de Colombia/Documentos/R/win-library/4.1/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/morea/OneDrive - Universidad Nacional de Colombia/Documentos/R/win-library/4.1/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Overwritten PROJ_LIB was C:/Users/morea/OneDrive - Universidad Nacional de Colombia/Documentos/R/win-library/4.1/rgdal/proj##
## Attaching package: 'rgdal'## The following object is masked from 'package:terra':
##
## projectlibrary(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.2 v dplyr 1.0.7
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 1.4.0 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::extract() masks terra::extract(), raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x ggplot2::layer() masks latticeExtra::layer()
## x dplyr::select() masks raster::select()
## x dplyr::src() masks terra::src()library(mapview)
library(sf)## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1library(elevatr)3. Get Raster Elevation Data
get_elev_raster() se requiere de un archivo con sistema de referencias coordenado
El nivel de “zoom” o la calidad (z) es por default 9
Creamos el simple feature a partir del shp del DANE
munic <- st_read("C:/Users/morea/OneDrive - Universidad Nacional de Colombia/Documentos/UNAL/2021 - 1S/GeomaticaBasica/MGN2017_68_SANTANDER/68_SANTANDER/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")## Reading layer `MGN_MPIO_POLITICO' from data source
## `C:\Users\morea\OneDrive - Universidad Nacional de Colombia\Documentos\UNAL\2021 - 1S\GeomaticaBasica\MGN2017_68_SANTANDER\68_SANTANDER\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 87 features and 9 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.52895 ymin: 5.707536 xmax: -72.47706 ymax: 8.14501
## Geodetic CRS: WGS 84head(munic)municelevation <- get_elev_raster(munic,z=10)## Mosaicing & Projecting## Note: Elevation units are in meters.
## Note: The coordinate reference system is:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]elevation## class : RasterLayer
## dimensions : 4598, 3603, 16566594 (nrow, ncol, ncell)
## resolution : 0.0006831156, 0.0006831156 (x, y)
## extent : -74.53125, -72.06998, 5.266202, 8.407168 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : fileb30407f906.tif
## names : fileb30407f906
## values : -32768, 32767 (min, max)munic_sp <- as(munic, "Spatial")Para plotearlo se hace una conversion en el data frame, es una visualizacion estandar con colores feos, datos que varian de 0 a 5000 con los nombres de los municipios
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)
Se toman los datos que corresponden al departamento con la funcion crop
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 : 3568, 3004, 10718272 (nrow, ncol, ncell)
## resolution : 0.0006831156, 0.0006831156 (x, y)
## extent : -74.5292, -72.47712, 5.707495, 8.144852 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : r_tmp_2021-07-17_235619_2864_10795.grd
## names : fileb30407f906
## values : -1608, 4512 (min, max)##5. Reprojectando los datos de elevacion Calculo de la pendiente y la orientacion se necesita reproyectar, una de las formas mas convencionales es con un sistema de coordenadas plano
library(sp)
spatialref <- CRS(SRS_string="EPSG:32618")pr3 <- projectExtent(elev_crop, crs=spatialref)#proyecta la zona sin datos
res(pr3) <- 180#define el tamñano de la celda
rep_elev <- projectRaster(elev_crop, pr3)#hace la proyeccionrep_elev## class : RasterLayer
## dimensions : 1502, 1264, 1898528 (nrow, ncol, ncell)
## resolution : 180, 180 (x, y)
## extent : 551863.8, 779383.8, 630819.1, 901179.1 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : fileb30407f906
## values : -374.5357, 4507.428 (min, max)(rep_munic = spTransform(munic_sp,spatialref))## class : SpatialPolygonsDataFrame
## features : 87
## extent : 552105.4, 778842.2, 630996, 900533.3 (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 : 68, 68001, AGUADA, 1539, 19.69444164, 2017, SANTANDER, 0.270018178002, 0.00160984367921
## max values : 68, 68895, ZAPATOCA, Ordenanza 33 de 1968, 3174.27867176, 2017, SANTANDER, 4.36038638488, 0.259459176725plot(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.2)
writeRaster(rep_elev, filename="C:/Users/morea/OneDrive - Universidad Nacional de Colombia/Documentos/UNAL/2021 - 1S/GeomaticaBasica/MGN2017_68_SANTANDER/68_SANTANDER/ADMINISTRATIVO/MGN_MPIO_POLITICO2.tif", datatype='INT4S', overwrite=TRUE)##6. Estadisticas basicas para los datos de elevación
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))##7.Obtencion de las variables geomorfometricas
slope = terrain(rep_elev,opt='slope', unit='degrees')
aspect = terrain(rep_elev,opt='aspect',unit='degrees')
hill = hillShade(slope,aspect,40,315) #sombra#plot de elevacion
plot(rep_elev,main="DEM for Santander municipalities [meters]", 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.2)
#Plot de la pendiente
plot(slope,main="Slope for Santander municipalities [degrees]", 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.2)
#el aspecto indica la orientacion de la pendiente, 0 es plano, si esta inclinado hacia el norte es 360 si esta orientado al oriente es de 45 y si cae al sur es de 180 si cae al oeste es**
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 DEM w/ hillshade, sombreado combinado que da sensacion de relieve
plot(hill,
col=grey(1:100/100), # create a color ramp of grey colors for hillshade
legend=FALSE, # no legend, we don't care about the grey of the hillshade
main="DEM for Santander municipalities [m]",
axes=FALSE) # makes for a cleaner plot, if the coordinates aren't necessary
plot(rep_elev,
axes=FALSE,
col=terrain.colors(12, alpha=0.6), add=TRUE) # color method
# sets how transparent the object will be (0=transparent, 1=not transparent)
plot(rep_munic, add=TRUE, lwd=0.5)
text(coordinates(rep_munic), labels=as.character(rep_munic$MPIO_CNMBR), cex=0.2)
##8. Mapeando datos de elevacion con rayshader
library(rayshader)elmat = raster_to_matrix(rep_elev)lista <- list("EL CARMEN DE CHUCURI", "SAN VICENTE DE CHUCURI", "SIMACOTA", "LANDÁZURI")
(algunos <- munic %>% filter(MPIO_CNMBR %in% unlist(lista) ))lista## [[1]]
## [1] "EL CARMEN DE CHUCURI"
##
## [[2]]
## [1] "SAN VICENTE DE CHUCURI"
##
## [[3]]
## [1] "SIMACOTA"
##
## [[4]]
## [1] "LANDÁZURI"algunos %>% mapview::mapview(zcol = "MPIO_CNMBR", legend=TRUE, selfcontained = FALSE)algunos_sp <- as(algunos, "Spatial")elev_crop = crop(elevation, algunos_sp)
plot(elev_crop, main="Cropped digital elevation model")
plot(algunos_sp, add=TRUE)
text(coordinates(algunos_sp), labels=as.character(algunos_sp$MPIO_CNMBR), cex=0.5)
# Project Raster
# to have control, we create one raster object
# using projectExtent (no values are transferred)
pr3 <- projectExtent(elev_crop, crs=spatialref)
# Adjust the cell size
res(pr3) <- 120
# now project
rep_elev <- projectRaster(elev_crop, pr3)(rep_algunos = spTransform(algunos_sp,spatialref))## class : SpatialPolygonsDataFrame
## features : 4
## extent : 599443.6, 688836.2, 681873.1, 791028 (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 : 68, 68235, EL CARMEN DE CHUCURI, 1799, 590.22625722, 2017, SANTANDER, 1.68427686489, 0.0482351082524
## max values : 68, 68745, SIMACOTA, Ordenanza 10 del 13 de Diciembre de 1974, 1116.36243904, 2017, SANTANDER, 3.11954907514, 0.0913288208216#Convert the DEM into a matrix:
elmat = raster_to_matrix(rep_elev)