El siguiente trabajo tiene como objetivo presentar la pendiente y la orientación de un DEM asignado. ## Creación del raster Antes de empezar, es necesario hacer una limpieza de la memoria de R, para que el programa trabaje más rápido y elimine datos innecesarios. Por lo tanto, el siguiente código servirá para esto.

rm(list=ls())

Ahora se va a llamar la librería necesaria, la cual será Raster.

library(raster)
## Loading required package: sp
library(sp)

Para crear el DEM,es necesario delimitar el raster con coordenadas (x, y) mínimas y máximas.

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

Si se quiere saber la cantidad de celdas que tiene el DEM se usa el siguiente código:

ncell(dem)
## [1] 16

El siguiente código funciona para saber la resolución espacial o tamaño de la celda, que por lo general es 5 metros por 5 metros.

res(dem)
## [1] 5 5

Asignación de los valores de elevación para cada celda del raster o DEM.

valores <- c(50, 45, 50, 48, 30, 29, 30, 29, 10, 9, 9, 10, 25, 23, 19, 21)

Valores de las celdas del DEM

(values(dem) <- valores)
##  [1] 50 45 50 48 30 29 30 29 10  9  9 10 25 23 19 21

A continuación se va a ploter el DEM con los valores de elevación:

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

Para asignar el sistema de coordenadas de referncia se usa el siguiente código: El sistema de coordenadas de referencia es el EPSG 3115 que hace referencia al MAGNA-SIRGAS/Colombia west zone y se encuentra este dato en (http://epsg.io/?q=colombia).

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

Cálculo de la pendiente y de la orientación del DEM

Para la pendiente es el siguiente código: Acá se pueden ver los atributos del DEM: la clase es tipo raster, las dimensiones, la resolución y el VALOR de la pendiente.

(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)

Ploteo de la pendiente, el cual muestra los valores dE la pendiente en cada celda

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

Ahora el cálculo de la orientación del DEM con sus atributos:

(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)

Ploteo de la orientación del DEM

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

FIN