Cálculo de pendientes y aspectos

Luisa Fernanda Carrión Ramírez

31/03/2019

Este cuaderno tiene como objetivo mostrar el proceso para calcular la pendiente y el aspecto empleando algunas herramientas de R. Matemáticamente hablando, la pendiente se define como el valor numérico asociado a la inclinación y dirección de un objeto con respecto a un plano. Por otro lado, el aspecto indica la orientación de la pendiente en el plano en relación con los puntos cardinales.

Creación del raster

rm(list=ls())

Antes de empezar, es necesario cargar la librería raster mediante el siguiente comando:

library(raster)

Para crear un Modelo Digital de Elevación (DEM) se debe empezar por delimitar las coordenadas en x y y de la siguiente manera:

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

Si se desea conocer cuántas celdas componen un DEM se escribe:

ncell(dem)

[1] 9

Con el fin de conocer la resolución espacial o el tamaño de la celda se emplea la siguiente función:

res(dem)

[1] 5 5

Luego, se deben asignar valores de elevación para el DEM:

valores <- c(50, 45, 50, 30, 30, 30, 8, 10, 10)

Para ejecutarlos, se emplea el siguiente comando:

(values(dem) <- valores)

[1] 50 45 50 30 30 30  8 10 10

La gráfica correspondiente al DEM se obtiene de la siguiente manera:

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

Se asigna un sistema de referencia de coordenadas:

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

Respuesta a algunas preguntas planteadas en la guía: - ¿Qué significa el código EPSG? Las siglas EPSG hacen referencia a European Petroleum Survey Group que es una organización que a partir de ciertos parámetros geodésicos determina unos códigos estándar o códigos EPSG. Estos se pueden empear para sistemas de coordenadas, esferoides, datos, entre otros. - ¿Dónde se pueden encontrar los códigos EPSG? Los códigos EPSG se encuentran compilados en la siguiente página:https://epsg.io/. Aunque, los shapefile generalmente vienen con un sistema de coordenadas establecido que puede ser leído por QGIS o ArcGIS. - ¿Cuál es el código EPSG para MAGNA Ciudad de Bogotá? El código EPSG asignado para MAGNA Ciudad de Bogotá es el 3116

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 +proj=tmerc +lat_0=4.596200416666666 +lon_0=-77.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      : slope 
values     : 75.25766, 75.25766  (min, max)

Para gráficar la pendiente:

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 +proj=tmerc +lat_0=4.596200416666666 +lon_0=-77.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      : aspect 
values     : 180.7538, 180.7538  (min, max)

Para graficar el aspecto:

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

Asignación

Para crear el DEM de la asignación se deben delimitar las coordenadas de la siguiente manera.

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

Si se quiere conocer el número de celdas se escribe:

ncell(demtar)

[1] 16

Con el fin de determinar la resolución espacial o el tamaño de la celda se escribe:

res(demtar)

[1] 5 5

Ahora, se asignan valores de elevación al DEM:

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

 [1] 50 45 50 48 30 29 30 29 10  9  9 10 25 23 19 21

Para graficar el DEM de la tarea:

plot(demtar, main = "DEM TAREA")
text(demtar)

También, se debe fijar un sistema de referencia de coordenadas:

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

Calcular la pendiente y el aspecto

(slopetar = terrain(demtar, '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 +proj=tmerc +lat_0=4.596200416666666 +lon_0=-77.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      : slope 
values     : 36.05503, 75.62313  (min, max)

Para graficar la pendiente:

plot(slopetar, main = "PENDIENTE")
text(slopetar)

(aspectotar = terrain(demtar, '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 +proj=tmerc +lat_0=4.596200416666666 +lon_0=-77.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      : aspect 
values     : 164.0546, 181.4688  (min, max)

Para graficar el aspecto:

plot(aspectotar, main="ASPECTO")
text(aspectotar)