1. Introducción

En este cuaderno se ilustra como calcular la pendiente y el aspecto de una zona geográfica específica, permitiendo hacer inferencia acerca de las condiciones del terreno analizado, debido a que estas variables presentan reñación directa con otras como tipo de suelos.

2. Creación de Raster

Es necesario correr la librería “raster”, previo al uso de los códigos

###
rm(list=ls())
library(raster)
## Loading required package: sp
##

Para el presente caso, calcularemos las variables de penidente y el aspecto, empleando datos de 9 celdas. Por lo anterior, se crean 3 filas y 3 columnas.

La resolución espacial puede incorporarse, delimitando los ejes coordinados x y y a valores específicos. Por ejemplo, en el presente caso se trabajan valores en el eje x desde 100 hasta 115, lo que quiere decir que para 3 celdas, se presentan 15 unidades, por lo cual a cada celda le corresponden 5 unidades, siendo este el ancho de la celda. Para el caso del largo, se presenta el mismo valor, por lo que cada celda presentará una resolución espacial de 5X5.

###
dem<-raster(ncol=3, nrow=3, xmn=100, xmx=115, ymn=100, ymx=115)
ncell(dem)
## [1] 9
res(dem)
## [1] 5 5
##

Posteriormente, se incorporan los valores de elevación para cada celda específica, donde el primer valor correspondea la celda en la primera fila y primera columna, y el último, a la última celda en la última columna, organizándoce de izquierda a derecha, y de forma descendente.

###
valores <- c(50, 45, 50, 30, 30, 30, 8, 10, 10)
(values(dem) <- valores)
## [1] 50 45 50 30 30 30  8 10 10
plot(dem, main = "DEM")
text(dem)

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

2. Cálculo de pendiente y aspecto

Los valores de “neighbors” o vecinos no son al azar, solo se pueden poner 4 u 8, básicamente porque cuando se diseña el algoritmo, se plantea el uso de las 4 celdas inmediatamente cercanas a la celda de en medio, o de las 8 que lo rodean por completo

###
(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)
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     : 180.7538, 180.7538  (min, max)
plot(aspecto, main = "Aspecto")
text(aspecto)

##

3. Ensayo

###
rm(list=ls())
##
###
dem <- raster(ncol=4, nrow=4, xmn=100, xmx=120, ymn=100, ymx=120)
ncell(dem)
## [1] 16
res(dem)
## [1] 5 5
valores <- c(50, 45, 50, 48, 30, 29, 30, 29, 10, 9, 9, 10, 25, 23, 19, 21)
(values(dem) <- valores)
##  [1] 50 45 50 48 30 29 30 29 10  9  9 10 25 23 19 21
plot(dem, main = "DEM")
text(dem)

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

Debido a que los algoritmos para el cálculo de la pendiente y el aspecto utilizan información de datos vecinos, es decir, como se mencionó anteiormente, de los “neighbors”, se obtienen 4 datos, debido a que se trabaja con 16 celdas.

###
(slope = terrain(dem, '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     : 36.05503, 75.62313  (min, max)
plot(slope, main = "Pendiente")
text(slope)

(aspecto = terrain(dem, '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     : 164.0546, 181.4688  (min, max)
plot(aspecto, main = "Aspecto")
text(aspecto)

##