My 5th Notebook DEM Antioquia
Juan Sebastian Camelo Garcia
library(lattice)
library(rgdal)## Loading required package: sp## 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/User/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/User/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/User/Documents/R/win-library/4.0/rgdal/projlibrary(raster)
library(elevatr)## Warning: package 'elevatr' was built under R version 4.0.5library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.3 v purrr 0.3.4
## v tibble 3.0.6 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 tidyr::extract() masks raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks raster::select()library(dplyr)
library(ggplot2)library(sf)## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1munic <- st_read("../GEOMATICA/05_ANTIOQUIA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")## Reading layer `MGN_MPIO_POLITICO' from data source `C:\Users\User\Desktop\GEOMATICA\05_ANTIOQUIA\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 125 features and 9 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -77.12783 ymin: 5.418558 xmax: -73.88128 ymax: 8.873974
## Geodetic CRS: WGS 84class(munic)## [1] "sf" "data.frame"head(munic)elevation <- get_elev_raster(munic, z = 10)## Note: Your request will download approximately 118.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(elevation, filename=file.path("../GEOMATICA/05_ANTIOQUIA/ADMINISTRATIVO/antioquia_dem_z10.tif"), format="GTiff", overwrite=TRUE)
}elevation## class : RasterLayer
## dimensions : 6134, 5667, 34761378 (nrow, ncol, ncell)
## resolution : 0.0006824976, 0.0006824976 (x, y)
## extent : -77.34375, -73.47604, 4.915657, 9.102097 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : C:/Users/User/AppData/Local/Temp/Rtmp0EKXvQ/file3650578f6c81.tif
## names : file3650578f6c81
## 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)
Realizamos un recorte
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 : 5063, 4757, 24084691 (nrow, ncol, ncell)
## resolution : 0.0006824976, 0.0006824976 (x, y)
## extent : -77.12808, -73.88144, 5.418657, 8.874143 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : file3650578f6c81
## values : -1609, 4758 (min, max)Reproyectamos para calcular la elevacion.
crs(elev_crop)## CRS arguments: +proj=longlat +datum=WGS84 +no_defslibrary(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 : 2126, 1998, 4247748 (nrow, ncol, ncell)
## resolution : 180, 180 (x, y)
## extent : 264195.4, 623835.4, 598928.5, 981608.5 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : file3650578f6c81
## values : -692.154, 4490.655 (min, max)(rep_munic = spTransform(munic_sp,spatialref))## class : SpatialPolygonsDataFrame
## features : 125
## extent : 265363.8, 623567.1, 598938.4, 981215 (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 : 05, 05001, ABEJORRAL, 1541, 15.83597275, 2017, ANTIOQUIA, 0.172962481259, 0.0012928296672
## max values : 05, 05895, ZARAGOZA, Ordenanza 9 de Diciembre 15 de 1983, 2959.36315089, 2017, ANTIOQUIA, 6.61177822771, 0.23894871119plot(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)
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 Antioquia 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(slope,main="Slope for Antioquia 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)
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 Antioquia 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)
library(rayshader)## Warning: package 'rayshader' was built under R version 4.0.5elmat = raster_to_matrix(rep_elev)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()
lista <- list("05045", "05490", "05837", "05147", "05172")
(algunos <- munic %>% filter(MPIO_CCDGO %in% unlist(lista) ))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_CCDGO), cex=0.5)
Debido a que algunos nombres contienen simbolos que no reconoce el software se realizo la lectura del listado con los codigos de los municipios de Apartado, Carepa, Chigorodo, Necocli y Turbo.
pr3 <- projectExtent(elev_crop, crs=spatialref)
res(pr3) <- 120
rep_elev <- projectRaster(elev_crop, pr3)(rep_algunos = spTransform(algunos_sp,spatialref))## class : SpatialPolygonsDataFrame
## features : 5
## extent : 265363.8, 346991.8, 824404.7, 964326.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 : 05, 05045, APARTADÓ, 1840, 387.27064143, 2017, ANTIOQUIA, 1.28765230797, 0.0316837676599
## max values : 05, 05837, TURBO, Ordenanza 7 de Noviembre 30 de 1967, 2959.36315089, 2017, ANTIOQUIA, 6.61177822771, 0.23894871119elmat = raster_to_matrix(rep_elev)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 18362)
##
## 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.19.2 sf_0.9-8 forcats_0.5.1 stringr_1.4.0
## [5] dplyr_1.0.5 purrr_0.3.4 readr_1.4.0 tidyr_1.1.3
## [9] tibble_3.0.6 ggplot2_3.3.3 tidyverse_1.3.0 elevatr_0.3.4
## [13] raster_3.4-5 rgdal_1.5-23 sp_1.4-5 lattice_0.20-41
##
## loaded via a namespace (and not attached):
## [1] fs_1.5.0 lubridate_1.7.10 webshot_0.5.2
## [4] doParallel_1.0.16 progress_1.2.2 httr_1.4.2
## [7] tools_4.0.4 backports_1.2.1 bslib_0.2.4
## [10] utf8_1.1.4 R6_2.5.0 KernSmooth_2.23-18
## [13] DBI_1.1.1 colorspace_2.0-0 manipulateWidget_0.10.1
## [16] withr_2.4.1 tidyselect_1.1.0 prettyunits_1.1.1
## [19] curl_4.3 compiler_4.0.4 cli_2.4.0
## [22] rvest_1.0.0 xml2_1.3.2 sass_0.3.1
## [25] scales_1.1.1 classInt_0.4-3 proxy_0.4-25
## [28] digest_0.6.27 rmarkdown_2.7 pkgconfig_2.0.3
## [31] htmltools_0.5.1.1 dbplyr_2.1.0 fastmap_1.1.0
## [34] highr_0.8 htmlwidgets_1.5.3 rlang_0.4.10
## [37] readxl_1.3.1 rstudioapi_0.13 shiny_1.6.0
## [40] jquerylib_0.1.3 generics_0.1.0 jsonlite_1.7.2
## [43] crosstalk_1.1.1 rayimage_0.5.1 magrittr_2.0.1
## [46] Rcpp_1.0.6 munsell_0.5.0 fansi_0.4.2
## [49] lifecycle_1.0.0 stringi_1.5.3 yaml_2.2.1
## [52] grid_4.0.4 parallel_4.0.4 promises_1.2.0.1
## [55] crayon_1.4.1 miniUI_0.1.1.1 haven_2.3.1
## [58] hms_1.0.0 knitr_1.31 pillar_1.5.0
## [61] codetools_0.2-18 reprex_2.0.0 glue_1.4.2
## [64] evaluate_0.14 modelr_0.1.8 vctrs_0.3.6
## [67] httpuv_1.5.5 foreach_1.5.1 cellranger_1.1.0
## [70] gtable_0.3.0 assertthat_0.2.1 xfun_0.21
## [73] mime_0.10 xtable_1.8-4 broom_0.7.5
## [76] e1071_1.7-6 later_1.1.0.1 class_7.3-18
## [79] iterators_1.0.13 units_0.7-1 rgl_0.105.22
## [82] ellipsis_0.3.1plot(cars)
Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.