Se quiere ilustrar el calculo de la pendiente y aspecto. Para esto se hara uso de los Modelos digitales de elevación (DEM). Las herramientas de análisis se pueden ejecutar en los DEM para mostrar nuevas superficies como pendientes y orientaciones como se observa a continuación.
library(raster)
## Loading required package: sp
rm(list=ls())
El primer paso es crear un DEM sencillo, Para ello es importante establecer las coordenasas del cuadro delimitador:
dem<-raster(ncol=3, nrow=3, xmn=100, xmx=115, ymn=100, ymx=115)
Celdas que componen el DEM creado:
ncell(dem)
## [1] 9
Para la resolucion espacial del DEM:
res(dem)
## [1] 5 5
Y despues se asignan valores al DEM:
valores <-c(50, 45, 50,30,30, 30, 8, 10, 10)
(values(dem)<-valores)
## [1] 50 45 50 30 30 30 8 10 10
Lo siguiente es plotear el DEM:
plot(dem, main = "DEM")
text(dem)
A continuacion se asigna un sistema de referencia de coordenadas, en este caso se utilico el MAGNA zona oeste de Colombia:
crs(dem)<- CRS('+init=epsg:3115')
EPSG son las siglas, en ingles, De European Petroleum Surey Group, que es una organización cientifica relacionada a la industria petrolera. Desde el 2005 hasta la actualidad el comite de geomatica de la international Association of Oil and gas Producers (OGP) son los encargados, es por ello que para descargar los datos se debe acceder a su plataforma virtual https://epsg.org/home.html. Para MAGNA de la Ciudad de Bogotá el codigo EPSG es 3116 y si se quiere la red urbana es posible usar el codigo 6247.
Lo primero es usar la siguiente funcion para trabajar la pendiente:
(slope = terrain(dem, 'slope', unit = 'degrees', neighbors=8))
## class : RasterLayer
## dimensions : 3, 3, 9 (nrow, ncol, ncell)
## resolution : 5, 5 (x, y)
## extent : 100, 115, 100, 115 (xmin, xmax, ymin, ymax)
## crs : +init=epsg:3115
## source : memory
## names : slope
## values : 75.25766, 75.25766 (min, max)
Con esto es posible realizar el plot:
plot(slope, main="Pendiente")
text(slope)
Para el caso del aspecto sera necesario usar la siguiente función:
(aspecto=terrain(dem, 'aspect', unit='degrees', neighbors=8))
## class : RasterLayer
## dimensions : 3, 3, 9 (nrow, ncol, ncell)
## resolution : 5, 5 (x, y)
## extent : 100, 115, 100, 115 (xmin, xmax, ymin, ymax)
## crs : +init=epsg:3115
## source : memory
## names : aspect
## values : 180.7538, 180.7538 (min, max)
Y de nuevo, se realiza el plot:
plot(aspecto, main='Aspecto')
text(aspecto)
Se calculara la pendiente y el aspecto para un DEM de 16 Celdas, donde el tamaño de la celda sera de 5:
adem<- raster(ncol=4, nrow=4, xmn=100, xmx=120, ymn=100, ymx=120)
ncell(adem)
## [1] 16
res(adem)
## [1] 5 5
avalor<-c(50, 45, 50, 48, 30, 29, 30, 29, 10, 9, 9, 10, 25, 23, 19, 21)
avalor
## [1] 50 45 50 48 30 29 30 29 10 9 9 10 25 23 19 21
values(adem)<-avalor
plot(adem, main = "DEM tarea")
text(adem)
Se utiliza el codigo EPSG 3116 para la zona Ciudad de Bogota:
crs(adem)<-CRS('+init=epsg:3116')
Con el DEM creado, lo siguiente es calcular la pendiente y el aspecto, para ello se realiza lo siguiente:
(aslope = terrain(adem, 'slope', unit='degrees', neighbors=8))
## class : RasterLayer
## dimensions : 4, 4, 16 (nrow, ncol, ncell)
## resolution : 5, 5 (x, y)
## extent : 100, 120, 100, 120 (xmin, xmax, ymin, ymax)
## crs : +init=epsg:3116
## source : memory
## names : slope
## values : 36.05503, 75.62313 (min, max)
plot(aslope, main = "Pendiente tarea")
text(aslope)
Para finalizar se calcula el aspecto:
(aaspect= terrain(adem, 'aspect', unit='degrees', neighbors=8))
## class : RasterLayer
## dimensions : 4, 4, 16 (nrow, ncol, ncell)
## resolution : 5, 5 (x, y)
## extent : 100, 120, 100, 120 (xmin, xmax, ymin, ymax)
## crs : +init=epsg:3116
## source : memory
## names : aspect
## values : 164.0546, 181.4688 (min, max)
y se realiza el plot:
plot(aaspect, main = "Aspecto tarea")
text(aaspect)
Como se observa en los plot mostrados en el documento, el calcula de la pendiente muestra que el terreno no es regular, es decir que presenta una inclinacion. En cambio el aspecto nos muestra que el modelo se orienta principalmente hacia el sur.