Slope and Aspect Calculation

Introduction

This is an R Markdown (http://rmarkdown.rstudio.com) Notebook to illustrate slope and aspect calculation. In mathematics, the slope or gradient of a plane is a number that describes both the direction and the steepness of the plane. On the another hand, aspect refers to the orientation of the plane slope in relation to the cardinal points. In short, aspect is the compass direction that a slope faces. Before going ahead, read the following two topics: (i) How slope works? (http://desktop.arcgis.com/en/arcmap/10.6/tools/spatial-analyst-toolbox/how-slope-works.htm) (ii) How aspect works? (http://desktop.arcgis.com/en/arcmap/10.6/tools/spatial-analyst-toolbox/how-aspect-works.htm) Take your time to review the links explaining slope and aspect algorithms. Then, take notes on your notebook. Once you understand how the two algorithms works, proceed to the next section.

Raster creation

Let’s create a toy DEM. Note that is necessary to pass bounding box coordinates (minimum and maximum along the two axes, x and y):

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

How many cells compose DEM?

ncell(dem)

[1] 9

What is the spatial resolution (i.e. cell size) of DEM?

res(dem)

[1] 5 5

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)

Let’s assign a coordinate reference system to DEM:

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

What does the EPSG code mean?

EPSG stands for European Petroleum Survey Group and is an organization that maintains a geodetic parameter database with standard codes, the EPSG codes, for coordinate systems, datums, spheroids, units and such alike. Additionally, this database contains the parameters for these objects or - if they cannot easily be expressed as values - at least references to where such parameters can be found.

Where can EPSG codes be found?

We can find EPSG codes at https://epsg.io/

What is the EPSG code for MAGNA Ciudad de Bogota?

The EPSG cofe for MGNA Ciudad Bogota is: 7458

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)

Assignment

Write R code to calculate slope and aspect for the following DEM:

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

ncell(ndem)

[1] 16

res(ndem)

[1] 5 5

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

(values(ndem) <- valores)

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

plot(ndem, main = "NDEM")
text(ndem)

Let’s assign a coordinate reference system to NDEM:

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

(slope = terrain(ndem, '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(ndem, '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)