library(rasterVis)
## Loading required package: lattice
library(raster)
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.1.2
library(rgl)
library(rgdal)
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-27, (SVN revision 1148)
## 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/ADMIN/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/ADMIN/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 sp or rgdal.
## Overwritten PROJ_LIB was C:/Users/ADMIN/Documentos/R/win-library/4.1/rgdal/proj
library(elevatr)
library(sf)
## Warning: package 'sf' was built under R version 4.1.2
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.4 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.1.0 v forcats 0.5.1
## Warning: package 'readr' was built under R version 4.1.2
## -- 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(mapview)
library(ggplot2)
Importar los archivos
- ShapeFile (*.shp) de municipios de Arauca
- Raster (*.tif) DEM de Arauca
mpio = st_read('../../arauca/81_ARAUCA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp')
## Reading layer `MGN_MPIO_POLITICO' from data source
## `C:\Users\ADMIN\Documentos\CARLOS\UNAL\Asignaturas\Geomatica\arauca\81_ARAUCA\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 10 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -72.36662 ymin: 6.036228 xmax: -69.42756 ymax: 7.104381
## Geodetic CRS: WGS 84
elev = raster("arauca_dem_z10.tif")
Mapa: superposición de Arauca sobre DEM
mpio_sp = as(mpio, 'Spatial')
plot(elev, main="DEM [metros]")
plot(mpio_sp, col=NA, border="black", add=TRUE)
text(coordinates(mpio_sp), col="black", cex=0.6,
labels=as.character(mpio$MPIO_CNMBR))
Mapa: recorte del DEM para la region Arauca
elev_crop = crop(elev, mpio_sp)
elev_crop = mask(elev_crop, mpio_sp)
plot(elev_crop, main="Elevación recortada")
plot(mpio_sp, add=TRUE)
text(coordinates(mpio_sp), col="black", cex=0.6,
labels=as.character(mpio$MPIO_CNMBR))
Pendiente
slope = terrain(elev_crop,opt='slope', unit='degrees')
aspect = terrain(elev_crop,opt='aspect',unit='degrees')
hill = hillShade(slope,aspect,40,315)
Mapa: pendientes del departamento de Arauca
plot(slope, main="Pendiente para municipos de Arauca [grados]",
col=topo.colors(6,alpha=0.6))
plot(mpio_sp, add=TRUE, lwd=0.5)
text(coordinates(mpio_sp), col="black", cex=0.6,
labels=as.character(mpio$MPIO_CNMBR))
Mapa: Aspecto del relieve departamento de Arauca
plot(aspect, main="Pendiente para municipos de Arauca [grados]",
col=rainbow(10, alpha = 0.6))
plot(mpio_sp, add=TRUE, lwd=0.5)
text(coordinates(mpio_sp), col="black", cex=0.6,
labels=as.character(mpio$MPIO_CNMBR))
Pendiente y Aspecto unicamente para el municipio de Tame, Arauca
- Filtrado del shapefile y recortado del raster para el municipio
- Calculo de pendiente para el municipio
mpio_tame = filter(mpio, MPIO_CNMBR == 'TAME')
mpio_tame_sp = as(mpio_tame, 'Spatial')
elev_tame_crop = crop(elev, mpio_tame_sp)
elev_tame_crop = mask(elev_tame_crop, mpio_tame_sp)
slope_tame = terrain(elev_tame_crop, opt='slope', unit='degrees')
aspect_tame = terrain(elev_tame_crop, opt='aspect',unit='degrees')
hill_tame = hillShade(slope_tame, aspect_tame, 40, 315)
- Conversion de raster a tabla para graficar
elev_tame_crop_tbl = rasterToPoints(elev_tame_crop) %>%
as.tibble()
## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
head(elev_tame_crop_tbl)
## # A tibble: 6 x 3
## x y arauca_dem_z10
## <dbl> <dbl> <dbl>
## 1 -71.8 6.70 270
## 2 -71.7 6.70 258
## 3 -71.7 6.70 256
## 4 -71.7 6.70 256
## 5 -71.7 6.70 252
## 6 -71.7 6.70 252
aspect_tame_crop_tbl = rasterToPoints(aspect_tame) %>%
as.tibble()
head(aspect_tame_crop_tbl)
## # A tibble: 6 x 3
## x y aspect
## <dbl> <dbl> <dbl>
## 1 -71.7 6.70 71.6
## 2 -71.7 6.70 207.
## 3 -71.8 6.70 45.0
## 4 -71.8 6.70 71.6
## 5 -71.8 6.70 71.6
## 6 -71.8 6.70 71.6
- Mapa: DEM de Tame usando ggplot2
ggplot()+
geom_sf(data = st_as_sf(mpio_tame), fill = NA)+
geom_tile(data = elev_tame_crop_tbl,
aes(x, y, fill = arauca_dem_z10))+
scale_fill_gradientn(colors = terrain.colors(15))+
labs(x = 'Longitud', y = 'Latitud', fill = 'Altitud (msnm)')+
theme_bw()+
theme(legend.position = 'bottom',
legend.key.width = unit(3, 'line'))
- Mapa: Aspecto del relieve de Tame usando ggplot2
ggplot()+
geom_sf(data = st_as_sf(mpio_tame), fill = NA)+
geom_tile(data = aspect_tame_crop_tbl,
aes(x, y, fill = aspect))+
scale_fill_gradientn(colors = rainbow(15, alpha = 0.5))+
labs(x = 'Longitud', y = 'Latitud',
title = 'Pendiente municipio de Tame, Arauca [grados]')+
theme_bw()+
theme(legend.position = 'bottom',
legend.key.width = unit(3, 'line'))
