Slope and Aspect Calculation

By Valentina Valle Velasco

April 1st, 2020

This is an R Markdown Notebook in order to illustrate slop and aspect calculation.
Slop is a tool that help us to identify the steepness in a cell of a raster surface. This can be calculated in two units, degrees and percent.
Aspect is a tool that is useful for identify the direction the downhill slope faces.

Now: Raster creation

rm(list = ls())

#library(raster)

Let´s create a toy DEM.

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

If you want to know how many cells compose DEM, use the following code

ncell(dem)

[1] 9

If you want to know what is the spatial resolution of DEM, use this

res(dem)

[1] 5 5

Now, we are going to assign elevation values to 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

Let’s plot the DEM along the elevation values:

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

No, we are going to assign a coordinate reference system to DEM:

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

What does mean the EPSG code?

EPSG is the acronym for European Petroleum Survey Group, organization related to the oil industry in Europe. This body was made up of specialists in geodesy, topography and cartography applied to the exploitation area and developed a repository of geodetic parameters that contains information on ancient and modern (geocentric) reference systems (frames), cartographic projections and ellipsoids from around the world.

Where can EPSG codes be found?

That type of codes can be found on epsg.io

What is the EPSG code for MAGNA Ciudad de Bogota?

The code is EPSG:3116, MAGNA-SIRGAS/ Colombia Bogota zone.

Slope and aspect calculation

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

Exercise

We are going to calculate slope and aspect for the following DEM:

[1] 16                          
[1] 5 5

So, here are all the necessary codes:

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

ncell(exe)

[1] 16

res(exe)

[1] 5 5

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

(values(exe) <- valores)

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

plot(exe, main = "Exercise")   
text(exe)

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

Now: Calculate SLOPE

(slope = terrain(exe, '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 = "SLOPE")    
text(slope)

(aspecto = terrain(exe, '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)

Now: Calculate ASPECT

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