cuaderno5-2 elevacion
Juan Guillermo Narvaez Diaz
29/4/2021
library(rasterVis)## Warning: package 'rasterVis' was built under R version 4.0.5## Loading required package: raster## Warning: package 'raster' was built under R version 4.0.5## Loading required package: sp## Loading required package: terra## Warning: package 'terra' was built under R version 4.0.5## terra version 1.1.4## Loading required package: lattice## Loading required package: latticeExtra## Warning: package 'latticeExtra' was built under R version 4.0.5library(raster)
library(rgl)## Warning: package 'rgl' was built under R version 4.0.5library(rgdal)## Warning: package 'rgdal' was built under R version 4.0.5## 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/Juan Narvaez/Documents/R/win-library/4.0/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/Juan Narvaez/Documents/R/win-library/4.0/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/Juan Narvaez/Documents/R/win-library/4.0/rgdal/proj##
## Attaching package: 'rgdal'## The following object is masked from 'package:terra':
##
## projectlibrary(elevatr)## Warning: package 'elevatr' was built under R version 4.0.5library(sf)## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.3 v purrr 0.3.4
## v tibble 3.1.0 v dplyr 1.0.5
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 1.4.0 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::collapse() masks terra::collapse()
## x dplyr::desc() masks terra::desc()
## x tidyr::expand() masks terra::expand()
## x tidyr::extract() masks terra::extract(), raster::extract()
## x tidyr::fill() masks terra::fill()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x ggplot2::layer() masks latticeExtra::layer()
## x dplyr::near() masks terra::near()
## x tidyr::pack() masks terra::pack()
## x dplyr::select() masks terra::select(), raster::select()
## x tidyr::separate() masks terra::separate()
## x purrr::transpose() masks terra::transpose()library(mapview)## Warning: package 'mapview' was built under R version 4.0.5(munici <- sf::st_read("geomatica/94_GUAINIA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))## Reading layer `MGN_MPIO_POLITICO' from data source `C:\Users\Juan Narvaez\Documents\geomatica\94_GUAINIA\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 9 features and 9 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -70.94249 ymin: 1.165633 xmax: -66.84722 ymax: 4.045026
## geographic CRS: WGS 84## Simple feature collection with 9 features and 9 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -70.94249 ymin: 1.165633 xmax: -66.84722 ymax: 4.045026
## geographic CRS: WGS 84
## DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR
## 1 94 94001 INÍRIDA
## 2 94 94343 BARRANCO MINA
## 3 94 94663 MAPIRIPANA
## 4 94 94883 SAN FELIPE
## 5 94 94884 PUERTO COLOMBIA
## 6 94 94885 LA GUADALUPE
## 7 94 94886 CACAHUAL
## 8 94 94887 PANÁ-PANÁ (Campo Alegre)
## 9 94 94888 MORICHAL (Morichal Nuevo)
## MPIO_CRSLC MPIO_NAREA
## 1 Decreto 1593 de Agosto 5 de 1974 15970.268
## 2 Resolución 83 del 1 de Febrero de 1988 que aprueba con modi 9467.804
## 3 ACUERDO COMISARIAL 018 DEL 25 DE AGOSTO DE 1990 4927.971
## 4 Resolución 83 2926.365
## 5 Resolución 83 15700.545
## 6 Resolución 83 1178.868
## 7 Resolución 83 2335.323
## 8 Resolución 83 de 1988 10227.139
## 9 1988 8554.996
## MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area geometry
## 1 2017 GUAINÍA 9.256685 1.2869895 POLYGON ((-67.67638 3.91228...
## 2 2017 GUAINÍA 6.825914 0.7652557 POLYGON ((-68.91332 3.68215...
## 3 2017 GUAINÍA 5.010789 0.3990354 POLYGON ((-70.10453 3.38436...
## 4 2017 GUAINÍA 4.029719 0.2346116 POLYGON ((-67.34976 2.50451...
## 5 2017 GUAINÍA 9.166487 1.2632740 POLYGON ((-67.5385 3.177568...
## 6 2017 GUAINÍA 1.750810 0.0943385 POLYGON ((-66.96243 1.66858...
## 7 2017 GUAINÍA 2.364174 0.1876089 POLYGON ((-67.49529 3.75681...
## 8 2017 GUAINÍA 7.805305 0.8250362 POLYGON ((-68.77132 2.42182...
## 9 2017 GUAINÍA 6.916544 0.6917881 POLYGON ((-69.46493 2.88960...elevacion<- get_elev_raster(munici,z=10)## Note: Your request will download approximately 124.2Mb.## Mosaicing & Projecting## Note: Elevation units are in meters.
## Note: The coordinate reference system is:
## GEOGCRS["WGS 84 (with axis order normalized for visualization)",
## 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["geodetic longitude (Lon)",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]if (require(rgdal)){
rf <- writeRaster(elevacion, filename=file.path("C:/Users/Juan Narvaez/Documents/geomatica/Guainia_Dem_z10.tif"), format="GTiff", overwrite=TRUE)
}munici_sp <- as(munici, "Spatial")
spatialref <- CRS(SRS_string="EPSG:32618")
elev_crop = crop(elevacion, munici_sp)
pr3 <- projectExtent(elev_crop, crs=spatialref)
res(pr3) <- 180
rep_elev <- projectRaster(elev_crop, pr3)
(rep_munic = spTransform(munici_sp,spatialref))## class : SpatialPolygonsDataFrame
## features : 9
## extent : 951413.4, 1410054, 130083.7, 450762.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 : 94, 94001, BARRANCO MINA, 1988, 1178.86832903, 2017, GUAINÍA, 1.75081007547, 0.0943384971244
## max values : 94, 94888, SAN FELIPE, Resolución 83 del 1 de Febrero de 1988 que aprueba con modi, 15970.2680471, 2017, GUAINÍA, 9.25668481051, 1.28698950156###datos basicos estadisticos de elevacion generamos un histograma para explorar los datos del DEM
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))## metricas valores
## 1 mean 150.94975
## 2 min -268.23602
## 3 max 960.60735
## 4 std 42.12976obtencion de variables geomorfometricos
slope = terrain(rep_elev,opt='slope', unit='degrees')
aspect = terrain(rep_elev,opt='aspect',unit='degrees')
hill = hillShade(slope,aspect,40,315)obtenemos el Dem de los municipios
plot(rep_elev,main="DEM de los municios de Guainia [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.2)
pendiente de la zona
plot(slope,main="pendiente de los municipios de Guainia [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.2)
aspecto de la zona podremos ver donde se encuentran los picos y en que municipios
plot(aspect,main="aspectos de grados municipios [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.2)
se combinan los mapas para tener uno solo
plot(hill,
col=grey(1:100/100),
legend=FALSE,
main="DEM de municipios de Guainia [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)
Utilizaremos la libreria rayshader para hacer un mapeo de los datos de elevacion
library(rayshader)## Warning: package 'rayshader' was built under R version 4.0.5#Convert the DEM into a matrix:
elmat = raster_to_matrix(rep_elev)
#And we can add a raytraced layer from that sun direction as well:
elmat %>%
sphere_shade(texture = "desert") %>%
add_water(detect_water(elmat), color = "desert") %>%
add_shadow(ray_shade(elmat), 0.5) %>%
add_shadow(ambient_shade(elmat), 0) %>%
plot_map()
visualizaremos una zona mas pequeña , la zona norte y central del departamento
lista <- list("MORICHAL (Morichal Nuevo)", "INÍRIDA", "BARRANCO MINA", "MAPIRIPANA", "PANÁ-PANÁ (Campo Alegre)")
(algunos <- munici %>% filter(MPIO_CNMBR %in% unlist(lista) ))## Simple feature collection with 5 features and 9 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -70.94249 ymin: 1.692744 xmax: -67.53842 ymax: 4.045026
## geographic CRS: WGS 84
## DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR
## 1 94 94001 INÍRIDA
## 2 94 94343 BARRANCO MINA
## 3 94 94663 MAPIRIPANA
## 4 94 94887 PANÁ-PANÁ (Campo Alegre)
## 5 94 94888 MORICHAL (Morichal Nuevo)
## MPIO_CRSLC MPIO_NAREA
## 1 Decreto 1593 de Agosto 5 de 1974 15970.268
## 2 Resolución 83 del 1 de Febrero de 1988 que aprueba con modi 9467.804
## 3 ACUERDO COMISARIAL 018 DEL 25 DE AGOSTO DE 1990 4927.971
## 4 Resolución 83 de 1988 10227.139
## 5 1988 8554.996
## MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area geometry
## 1 2017 GUAINÍA 9.256685 1.2869895 POLYGON ((-67.67638 3.91228...
## 2 2017 GUAINÍA 6.825914 0.7652557 POLYGON ((-68.91332 3.68215...
## 3 2017 GUAINÍA 5.010789 0.3990354 POLYGON ((-70.10453 3.38436...
## 4 2017 GUAINÍA 7.805305 0.8250362 POLYGON ((-68.77132 2.42182...
## 5 2017 GUAINÍA 6.916544 0.6917881 POLYGON ((-69.46493 2.88960...algunos_sp <- as(algunos, "Spatial")
elev_crop = crop(elevacion, algunos_sp)
plot(elev_crop, main=" modelo digital recortado ")
plot(algunos_sp, add=TRUE)
text(coordinates(algunos_sp), labels=as.character(algunos_sp$MPIO_CNMBR), cex=0.5)
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 : 5
## extent : 951413.4, 1330940, 187882.8, 450762.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 : 94, 94001, BARRANCO MINA, 1988, 4927.97124104, 2017, GUAINÍA, 5.01078939041, 0.399035434995
## max values : 94, 94888, PANÁ-PANÁ (Campo Alegre), Resolución 83 del 1 de Febrero de 1988 que aprueba con modi, 15970.2680471, 2017, GUAINÍA, 9.25668481051, 1.28698950156lo veremos de una zona 3D
#Convert the DEM into a matrix:
elmat = raster_to_matrix(rep_elev)
#detect_water and add_water adds a water layer to the map:
elmat %>%
sphere_shade(texture = "desert") %>%
add_water(detect_water(elmat), color = "desert") %>%
add_shadow(ray_shade(elmat), 0.5) %>%
add_shadow(ambient_shade(elmat), 0) %>%
plot_map()
sessionInfo()## R version 4.0.4 (2021-02-15)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Spanish_Colombia.1252 LC_CTYPE=Spanish_Colombia.1252
## [3] LC_MONETARY=Spanish_Colombia.1252 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Colombia.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rayshader_0.24.5 mapview_2.9.0 forcats_0.5.1
## [4] stringr_1.4.0 dplyr_1.0.5 purrr_0.3.4
## [7] readr_1.4.0 tidyr_1.1.3 tibble_3.1.0
## [10] ggplot2_3.3.3 tidyverse_1.3.0 sf_0.9-7
## [13] elevatr_0.3.4 rgdal_1.5-23 rgl_0.105.22
## [16] rasterVis_0.50.1 latticeExtra_0.6-29 lattice_0.20-41
## [19] terra_1.1-4 raster_3.4-5 sp_1.4-5
##
## loaded via a namespace (and not attached):
## [1] fs_1.5.0 satellite_1.0.2 lubridate_1.7.10
## [4] doParallel_1.0.16 progress_1.2.2 webshot_0.5.2
## [7] RColorBrewer_1.1-2 httr_1.4.2 tools_4.0.4
## [10] backports_1.2.1 bslib_0.2.4 utf8_1.1.4
## [13] R6_2.5.0 KernSmooth_2.23-18 DBI_1.1.1
## [16] colorspace_2.0-0 manipulateWidget_0.10.1 withr_2.4.1
## [19] prettyunits_1.1.1 tidyselect_1.1.0 leaflet_2.0.4.1
## [22] curl_4.3 compiler_4.0.4 leafem_0.1.3
## [25] cli_2.3.1 rvest_1.0.0 xml2_1.3.2
## [28] sass_0.3.1 scales_1.1.1 classInt_0.4-3
## [31] hexbin_1.28.2 proxy_0.4-25 digest_0.6.27
## [34] rmarkdown_2.7 base64enc_0.1-3 jpeg_0.1-8.1
## [37] pkgconfig_2.0.3 htmltools_0.5.1.1 highr_0.8
## [40] dbplyr_2.1.0 fastmap_1.1.0 htmlwidgets_1.5.3
## [43] rlang_0.4.10 readxl_1.3.1 rstudioapi_0.13
## [46] shiny_1.6.0 jquerylib_0.1.3 generics_0.1.0
## [49] zoo_1.8-9 jsonlite_1.7.2 crosstalk_1.1.1
## [52] rayimage_0.5.1 magrittr_2.0.1 Rcpp_1.0.6
## [55] munsell_0.5.0 fansi_0.4.2 lifecycle_1.0.0
## [58] stringi_1.5.3 yaml_2.2.1 grid_4.0.4
## [61] parallel_4.0.4 promises_1.2.0.1 crayon_1.4.1
## [64] miniUI_0.1.1.1 haven_2.3.1 hms_1.0.0
## [67] knitr_1.31 pillar_1.5.1 stats4_4.0.4
## [70] codetools_0.2-18 reprex_1.0.0 glue_1.4.2
## [73] evaluate_0.14 modelr_0.1.8 foreach_1.5.1
## [76] png_0.1-7 vctrs_0.3.6 httpuv_1.5.5
## [79] cellranger_1.1.0 gtable_0.3.0 assertthat_0.2.1
## [82] xfun_0.22 mime_0.10 xtable_1.8-4
## [85] broom_0.7.5 e1071_1.7-5 later_1.1.0.1
## [88] class_7.3-18 viridisLite_0.3.0 iterators_1.0.13
## [91] units_0.7-0 ellipsis_0.3.1