Datos de elevaciĆ³n en R
# install.packages("rgdal")
# install.packages("raster")
# install.packages("elevatr")
# install.packages("rasterVis")
# install.packages("rgl")library(rasterVis)## Loading required package: raster## Loading required package: sp## Loading required package: lattice## Loading required package: latticeExtralibrary(raster)
library(rgl)
library(rgdal)## rgdal: version: 1.4-8, (SVN revision 845)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/yisus/Documents/R/win-library/3.6/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/yisus/Documents/R/win-library/3.6/rgdal/proj
## Linking to sp version: 1.4-1library(elevatr)OBTENER DATOS RASTER DE ELEVACIĆN
Primero se carga un archivo shapefile, en este caso el archivo utilizado en este Notebook ha sido tomado del Marco GeoestadĆstico Nacional para el departamento de Antioquia del aƱo 2017 disponible en el siguiente enlace: https://geoportal.dane.gov.co/servicios/descarga-y-metadatos/descarga-mgn-marco-geoestadistico-nacional/
list.files("C:/Users/yisus/Documents/Antioquia-Geomatica/Antioquia")## [1] "19-03-2020.qgz" "Anexo_Antioquia.csv"
## [3] "Anexo_Antioquia.xls" "Anexo_Antioquia2.csv"
## [5] "Anexo_Antioquia3.csv" "Dem.tif"
## [7] "Departamento.cpg" "Departamento.dbf"
## [9] "Departamento.prj" "Departamento.qpj"
## [11] "Departamento.shp" "Departamento.shx"
## [13] "Departamento.xlsx" "MGN_MPIO_POLITICO.cpg"
## [15] "MGN_MPIO_POLITICO.dbf" "MGN_MPIO_POLITICO.prj"
## [17] "MGN_MPIO_POLITICO.sbn" "MGN_MPIO_POLITICO.sbx"
## [19] "MGN_MPIO_POLITICO.shp" "MGN_MPIO_POLITICO.shp.xml"
## [21] "MGN_MPIO_POLITICO.shx" "mun_3116.gpkg"
## [23] "Municipios.cpg" "Municipios.dbf"
## [25] "Municipios.prj" "Municipios.qpj"
## [27] "Municipios.shp" "Municipios.shx"
## [29] "potw1726a.jpg" "Proyecto 12 de Marzo.qgz"
## [31] "Vias.dbf" "Vias.prj"
## [33] "Vias.qpj" "Vias.shp"
## [35] "Vias.shx"(munic <- shapefile("C:/Users/yisus/Documents/Antioquia-Geomatica/Antioquia/MGN_MPIO_POLITICO.shp"))## class : SpatialPolygonsDataFrame
## features : 125
## extent : -77.12783, -73.88128, 5.418558, 8.873974 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## 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.23894871119Atributos del objeto āmunicā:
head(munic)Para este caso seleccionĆ© la ciudad llamada āCiudad BolĆvarā con el fin de no repetir el Notebook de ejemplo.
CiudadBolivar <- munic[munic$MPIO_CNMBR=="CIUDAD BOLĆ\215VAR",]
plot(CiudadBolivar, main="Ciudad BolĆvar", axes=TRUE)
plot(munic, add=TRUE)
invisible(text(coordinates(munic), labels=as.character(munic$MPIO_CNMBR), cex=0.5))
Se descargĆ³ los datos de elevaciĆ³n.
elevation <- get_elev_raster(CiudadBolivar, z = 10)## Merging DEMs## Reprojecting DEM to original projection## Note: Elevation units are in meters.
## Note: The coordinate reference system is:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0Atributos del objeto āelevationā:
elevation## class : RasterLayer
## dimensions : 1035, 1545, 1599075 (nrow, ncol, ncell)
## resolution : 0.000687, 0.000683 (x, y)
## extent : -76.2925, -75.23108, 5.262264, 5.969169 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 161.7008, 4685.116 (min, max)plot(elevation, main="DEM descargado [metros]")
plot(CiudadBolivar, add=TRUE)
RECORTAR LOS DATOS PARA QUE COINCIDAN CON LA REGIĆN SELECCIONADA
Se escoge los datos correspondientes a Ciudad BolĆvar.
elev_crop = crop(elevation, CiudadBolivar)
plot(elev_crop, main="Modelo de elevaciĆ³n digital recortado")
plot(CiudadBolivar, add=TRUE)
VerificaciĆ³n del nuevo objeto:
elev_crop## class : RasterLayer
## dimensions : 261, 275, 71775 (nrow, ncol, ncell)
## resolution : 0.000687, 0.000683 (x, y)
## extent : -76.09945, -75.91053, 5.753341, 5.931604 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 648.3034, 3937.323 (min, max)REPROYECTAR LOS DATOS DE ELEVACIĆN
Se reproyecta con el fin de usar coordenadas de mapa en lugar de coordenadas geogrĆ”ficas, para esto busca la definiciĆ³n para esta referencia espacial in formato PROJ.4.
spatialref <- "+proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" Se procede reproyectar.
pr3 <- projectExtent(elev_crop, spatialref)res(pr3) <- 100rep_elev <- projectRaster(elev_crop, pr3)VerificaciĆ³n de la reproyecciĆ³n:
rep_elev## class : RasterLayer
## dimensions : 198, 210, 41580 (nrow, ncol, ncell)
## resolution : 100, 100 (x, y)
## extent : 775998.8, 796998.8, 1128282, 1148082 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## source : memory
## names : layer
## values : 650.9104, 3936.504 (min, max)Se reproyecta el SpatialPolygonsDataFrame que representa a Ciudad BolĆvar.
(rep_CiudadBolivar = spTransform(CiudadBolivar,spatialref))## class : SpatialPolygonsDataFrame
## features : 1
## extent : 776023.9, 796960.5, 1128331, 1148060 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +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
## value : 05, 05101, CIUDAD BOLĆĀVAR, 1869, 260.44610771, 2017, ANTIOQUIA, 0.708503874151, 0.0212422352127VisualizaciĆ³n:
plot(rep_elev, main="Reprojected Digital elevation model")
plot(rep_CiudadBolivar, add=TRUE)
ESTADĆSTICAS BĆSICAS DE 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))OBTENCIĆN DE VARIABLES GEOMORFOMĆTRICAS
Calcular la pendiente, el aspect y el sombreado.
slope = terrain(rep_elev,opt='slope', unit='degrees')
aspect = terrain(rep_elev,opt='aspect',unit='degrees')
hill = hillShade(slope,aspect,40,315)VisualizaciĆ³n de la elevaciĆ³n:
plot(rep_elev,main="DEM para Ciudad BolĆvar [metros]", col=terrain.colors(25,alpha=0.7))
VisualizaciĆ³n de la pendiente:
plot(slope,main="Pendiente para Ciudad BolĆvar [grados]", col=topo.colors(25,alpha=0.7))
VisualizaciĆ³n del aspect:
plot(aspect,main="Aspect para Ciudad BolĆvar [grados]", col=rainbow(25,alpha=0.7))
MAPEO DE DATOS DE ELEVACIĆN CON RAYSHADER
# install.packages("rayshader")library(rayshader)elmat = raster_to_matrix(rep_elev)elmat %>%
sphere_shade(texture = "imhof2") %>%
plot_map()
elmat %>%
sphere_shade(texture = "desert") %>%
add_water(detect_water(elmat), color = "desert") %>%
plot_map()
elmat %>%
sphere_shade(texture = "desert") %>%
add_water(detect_water(elmat), color = "desert") %>%
add_shadow(ray_shade(elmat), 0.5) %>%
plot_map()
OTRA FORMA DE VISUALIZACIĆN
# install.packages("jpeg")library(jpeg)getv=function(i,a,s){
ct = dim(i)[1:2]/2
sx = values(s)/90 * ct[1]
sy = values(s)/90 * ct[2]
a = values(a) * 0.01745
px = floor(ct[1] + sx * -sin(a))
py = floor(ct[2] + sy * cos(a))
template = brick(s,s,s)
values(template)=NA
cellr = px + py * ct[1]*2
cellg = px + py * ct[1]*2 + (ct[1]*2*ct[2]*2)
cellb = px + py * ct[1]*2 + 2*(ct[1]*2*ct[2]*2)
template[[1]] = i[cellr]
template[[2]] = i[cellg]
template[[3]] = i[cellb]
template = template * 256
template
}map=readJPEG("C:/Users/yisus/Documents/Antioquia-Geomatica/Antioquia/potw1726a.jpg")out = getv(map, aspect, slope)plotRGB(out, main = "Supposedly pretty mountains in Ciudad BolĆvar")