Cálculo de pendiente y aspecto

Introducción

Este es un cuaderno R Markdown (http://rmarkdown.rstudio.com) para ilustrar el cálculo de pendiente y orientación. En matemáticas, la pendiente o pendiente de un plano es un número que describe tanto la dirección como la inclinación del plano. Por otro lado, el aspecto se refiere a la orientación de la pendiente del plano en relación con los puntos cardinales. En resumen, el aspecto es la dirección de la brújula a la que se enfrenta una pendiente. Antes de continuar, lea los dos temas siguientes:

1. ¿Cómo funciona la pendiente? (http://desktop.arcgis.com/en/arcmap/10.6/tools/spatial-analyst-toolbox/how-slope-works.htm)

2. ¿Cómo funciona el aspecto? (http://desktop.arcgis.com/en/arcmap/10.6/tools/spatial-analyst-toolbox/how-aspect-works.htm)

rm(list=ls())

library(raster)

Loading required package: sp

Creación de ráster

Creemos un DEM de juguete. Tenga en cuenta que es necesario pasar las coordenadas del cuadro delimitador (mínimo y máximo a lo largolos dos ejes, “x” y “y”):

dem <- raster(ncol=3, nrow=3, xmn=100, xmx=115, ymn=100, ymx=115)

¿Cuántas células componen DEM?

ncell(dem)

[1] 9

¿Cuál es la resolución espacial (es decir, el tamaño de celda) de DEM?

res(dem)

[1] 5 5

Asignar valores de elevación a DEM:

valores <- c(50, 49, 53, 35, 36, 35, 9, 15, 16)

(values(dem) <- valores)

[1] 50 49 53 35 36 35  9 15 16

Tracemos el DEM a lo largo de los valores de elevación:

plot(dem, main = "DEM")
text(dem)

Asignemos un sistema de referencia de coordenadas a DEM:

crs(dem) <- CRS('+init=epsg:3115')

Cálculo de pendiente y aspecto

(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     : 74.71263, 74.71263  (min, max)

plot(slope, main = "Pendiente")
text(slope)

(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     : 183.9182, 183.9182  (min, max)

plot(aspecto, main = "Aspecto")
text(aspecto)

Asignación

creamos un nuevo dem:

dem2 <- raster(ncol=4, nrow=4, xmn=100, xmx=120, ymn=100, ymx=120)

ncell(dem2)

[1] 16

res(dem2)

[1] 5 5

valores2 <- c(50, 49, 53,50, 35, 36, 35, 30, 9, 15, 16, 12, 20, 24, 28, 22)
(values(dem2) <- valores2)

 [1] 50 49 53 50 35 36 35 30  9 15 16 12 20 24 28 22

plot(dem2, main="Calcule slope & aspect for this DEM")
text(dem2)

crs(dem2) <- CRS('+init=epsg:3115')
(slope2 = terrain(dem2, '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:3115 
source     : memory
names      : slope 
values     : 42.59043, 74.74518  (min, max)

(aspecto2 = terrain(dem2, '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:3115 
source     : memory
names      : aspect 
values     : 157.6199, 205.56  (min, max)

plot(aspecto2, main = "Aspecto")
text(aspecto2)

plot(slope2, main = "slope")
text(slope2)