Chapter 17 Raster data and terrain analysis
Terrain analysis involves the analysis of topographic features (geomorphology). The operations required typically take raster layers as input,although three dimensional point cloud data can also be used directly for some analyses. The input to most forms of terrain analyses are digital elevation models (DEMs). These can take various forms.Traditionally they were produced through interpolation between contours from a traditional map. DEM’s are now mainly derived from remote sensing involving some form of radar based sensor. For example
- SRTM (Shuttle Radar Topography mission): Global extent 90m resolution
- Lidar (Light detection and ranging): Local extent, varying resolution
Raw lidar data consists of a point cloud of returns with some additional information regarding the intensity of the return. The point cloud has two key reference points. The first return represents the top of vegetation or buildings and the last return should typically be the ground.
Lidar point clouds can be processed to form digital terrain models (DTMs) by taking the last returns and using smoothing and interpolation. Digital surface models use the first return. This processing can be achieved using a range of different software. For example QGIS can be linked to LASTools to process point clouds directly.
Terrain analysis uses algorithms to calculate a potentially large number of different derived layers from the digital terrsin models and the digital surface models. These typically include slope, aspect, hillshading and elements such as direct beam insolation and topographic wetness index.
There are many uses for terrain analysis. Hydrologists can use it to delimit watersheds and analyse flows. The layers derived from terrain analysis can be queried to produce information at points, along lines or within polygons. We can also carry out raster alegebra operations on the layers.
https://rpubs.com/dgolicher/lidar
To see some point clouds open http://plas.io/ and load the lidar data from the course web site.
Here are some point clouds of Arne. Download them locally then open with plas.io.
http://r.bournemouth.ac.uk:82/Quantitative_and_Spatial_Analysis/Week_2/point_clouds/
Processing the point cloud (not conducted in this class) leads to the production of a Digital surface model (DSM) and a Digital terrain model (DTM)
The digital surface model is obtained by smoothing the first returns in order to produce values for a two dimensional (x and y) raster grid placed on the point cloud. The third dimension is now the value of the pixel. The resolution can be chosen when producing the model. it can’t be finer than the typical gaps between points, but it can be coarser.
The digital terrain model is obtained by smoothing the last returns with some interpolation techniques included to smooth around hard surfaces in order to produce values that should correspond to the ground surface. In cities and towns this will be more challenging than in areas with vegetation, as returns from the actual ground may be difficult.
Archaeologists and geomorphologists are often interested in the visual aspect of terrain analysis. Techniques such as analytical hill shading can reveal subtle surface features that would be otherwise hidden from view. As Lidar penetrates through many types pf vegetation analysis of Lidar derived digital surface models have been used to find Mayan ruins beneath tropical forest.
In many ecological studies the main purpose behind running terrain analysis is to derive variables that can be combined with others within a quantiative analysis. For example if we place a set of quadrats down within a study area in order to look at the relationships between vegetation types and other variables we may want to know variables such as slope and aspect. Rather than try to measure this on the ground, we can easily add the values to the analysis if have derived them from a dtm and we have the coordinates of the quadrats using a gps. In many studies it is preferable to derive more informative variables than slope and aspect. For example slope and aspect combined determine how much direct beam insolation the surface recieves on a given day. This could be useful in a study of butterflies or small reptiles that may seek warm micro climates. Rather than just looking at slope, the position on a slope may be useful. There are a wide range of options that can be explored.
Raster layers can represent the same features that are represented by vector layers, but the information is held in a different format. For example land use classes can be converted from vector to raster.
17.1 Obtaining SRTM data
Before we go onto to look at Lidar data at the Arne we can look at data derived from the SRTM. The data from lidar or from SRTM point clouds are converted into raster layers.
The simplest derived layer is the digital terrain model. This is sometimes called a digital elevation model, but digital elevations can also be derived from the first return.
When working with raster data it is very important to be aware of scale. The resolution of the raster layer is determined by the size of the grid cell. As high resolution layers quickly become very large, analysis at regional scales uses lower resolution data. This leads to results that look good when viewed at the right scale, but are less suitable for local, site specific analysis. The solution to this problem is to run analyses at a scale which matches the questions being asked while being prepared to run similar code on local high resolution data sets to look at specific sites.
The code shown below will produce maps if all the chunks are run in the correct order. The maps have not been shown in the handout. The exercise associated with this chapter is to work through the code and investigate the maps which it produces yourself. You may want to set up a project to do this. The analysis is not part of the assignment, but will help you understand how the Arne data has been processed.
library(tidyverse)
library(raster)
library(giscourse)
library(mapview)
library(sf)
library(rnaturalearth)
library(dplyr)
library(RColorBrewer)
pal1<-terrain.colors(100)
pal2<-gray(0:100 / 100)## The code below will produce the data without the server. However it takes a while to run.
## The data has been added to the aqm package
## Removing the # comments on the lines would load the same data off the server
# counties<-ne_states(country = 'United Kingdom', returnclass = c("sf"))
# counties %>% filter (name=='Dorset') %>% as("Spatial")->dorset
# dem<-elevatr::get_elev_raster(dorset,z=11)
# dem<- raster::mask(dem,dorset)
# dem<-raster::projectRaster(dem,crs=27700)
## Load: Only these lines are needed on the server.
library(aqm)
data("dorset_dem")17.2 Attributes of a raster layer
You can check the attributes of a raster layer just by typing its name.
dem## class : RasterLayer
## dimensions : 2736, 5582, 15272352 (nrow, ncol, ncell)
## resolution : 20.2, 31.9 (x, y)
## extent : 329697.5, 442453.9, 55579.94, 142858.3 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs
## source : memory
## names : file7b6021cb1386
## values : -1.960342, 274 (min, max)Notice the number of pixels, the resolution in meters and coordinate reference system. This is British National grid.
The area of each pixel in square meters can be calculated by multiplying the two elements of the resolution. We can use this later for calculating total areas.
parea<-res(dem)[1]*res(dem)[2]
parea## [1] 644.3817.3 Quick plot of a raster layer
Rasters can be plotted simply.
plot(dem)Its also simple to get a histogram of the raster values.
hist(dem)If we wanted to find the total area of Dorset over 100m in square km we can do this.
vals <- getValues(dem)
sum(vals>100, na.rm=TRUE) * parea/1000000## [1] 864.106517.4 Building a hillshaded layer
In order to build a hillshaded layer we first need to derive the slope and aspect from the DTM in radians, then run the hillsheade algorithm.
dslope<-terrain(dem, opt='slope', unit='radians')
daspect<-terrain(dem, opt='aspect',unit='radians')
dhill<-hillShade(dslope,daspect, angle=10)We can now look at the result. Run the code yourself to see it.
plot(dhill, col = gray(0:100 / 100), legend = FALSE)
plot(dem, col = pal1, alpha = 0.4, legend = TRUE, add = TRUE)
plot(dorset,add=TRUE)This can be shown as a dynamic mapview.
map<-mapview(dem,col.regions=pal1, legend=FALSE, na.color="transparent") +
mapview(dhill,col.regions=pal2,na.color="transparent", legend=FALSE)
## Add full screen control from leaflet.extras
library(leaflet.extras)
map@map <- map@map %>%
addFullscreenControl() This will show the results. Not run here. Run the code yourself making sure that all lines above have been run.
map17.5 Finding contours
contours<-rasterToContour(dem,levels=c(50,100,200,300))A simple map can be plotted and the contours added.
plot(dem)
plot(contours, add= TRUE)The contours can be added to the mapview map.
map + mapview(contours, color="black", legend=FALSE)17.6 Raster algebra.
When you have a set of raster layers to work with it is possible to combine them to form new layers using either sophisticated algorithms based on physical rules, such as those used to calculate insolation and topographic wetness index, or simpler “bespoke” rules that you can design. For example, say we want to find the steep south facing slopes in Dorset.
We can first derive layers representing slope and aspect in degrees.
slope<-terrain(dem, opt='slope', unit='degrees')
aspect<-terrain(dem, opt='aspect', unit='degrees')Now we can design a rule. Let’s assume our definition of steep is a slope over 5 degrees and aspect lies between 130 and 230 degrees.
## Find slopes and set other pixels to NA
steep_south_slope<-(slope>5 & aspect> 130 & aspect< 230)
steep_south_slope[steep_south_slope<1]<-NAThis will show the results.
map + mapview(steep_south_slope,legend=FALSE, na.color ='transparent') ->map
mapIf we wanted to only look at calcareous grassland on south facing slopes we could include another layer.
data(calc_dorset)
steep_south_calc<-(steep_south_slope & calc==1)
steep_south_calc[steep_south_calc<1]<-NAmap +
mapview(steep_south_calc,legend=FALSE, na.color ='transparent')17.7 Additional algorithms
Over and above the raster processing included in the raster package R provides an interface to two major pieces of software that can run complex algorithms used for terrain analysis. These are SAGA and GRASS. The code to run the algorithms is rather complex and it is not shown here. Some care is required to use algorithms on large rasters as processing times can be long and lead to the server freezing. However once derived they can be used for further analysis.
SAGA can be installed as a standalone GIS on any laptop. It also provides plug in algorithms for QGIS. The version installed on the server is SAGA 7.3. This provides all the algorithms used by the desktop program, but they are called through a command line in R. Thus the SAGA GUI is not used for visualisation.
http://www.saga-gis.org/en/index.html
GRASS (version 7.8) is also installed on the server and can be called through R. GRASS is one of the oldest and most powerful sets of algorithms for raster processing and is used in many scientific papers.
https://grass.osgeo.org/grass78/manuals/
GRASS can also be used from QGIS and also has its own standalone GUI.
# Adds layers from GRASS and SAGA processing
## Code used to run the algorithms not shown
data(dorset_tm)GRASS can be used to analyse the digital terrain model and find the most likely paths that would be taken by water flowing over it. This can be used to delimit watersheds and to find unmapped stream paths. The accuracy of this depends on the reliability of the DEM. It does not work perfectly for coarse resolution rasters, such as this one, but it can still be useful as a rough guide.
map+mapview(streams) Direct beam radition recieved over a given time period by a pixel under clear sky conditions can be calculated using SAGA. This can be useful in many more advanced analyses, particularly when finding “warm spots” on a landscape that may be attractive to butterflies and reptiles.
mapview(insol )17.8 Exercise
Work through the code lines and investigate the maps produced.
17.9 Extra examples
These are just added to illustrate that code can be applied to other areas easily once designed. You do not need to run the chunks except as an illustration.
17.9.1 Snowdon
## Just for fun. Get the DEM for Snowdon
g<-tmaptools::geocode_OSM("Snowdon") # Find the place
d<-data.frame(id=1,x=g$coords[1],y=g$coords[2]) # Make a data frame
d<-st_as_sf(d,coords=c("x","y"),crs=4326) # Turn spatial
dem<-elevatr::get_elev_raster(d, z=14) # Get the DEM
dslope<-terrain(dem, opt='slope', unit='radians')
daspect<-terrain(dem, opt='aspect',unit='radians')
dhill<-hillShade(dslope,daspect, angle=10)
map<-mapview(dem,col.regions=pal1, legend=FALSE, na.color="transparent") +
mapview(dhill,col.regions=pal2,na.color="transparent", legend=FALSE)
map@map <- map@map %>%
addFullscreenControl()
map17.9.2 Corfe castle
g<-tmaptools::geocode_OSM("Corfe Castle") # Find the place
d<-data.frame(id=1,x=g$coords[1],y=g$coords[2]) # Make a data frame
d<-st_as_sf(d,coords=c("x","y"),crs=4326) # Turn spatial
dem<-elevatr::get_elev_raster(d, z=14) # Get the DEM
dslope<-terrain(dem, opt='slope', unit='radians')
daspect<-terrain(dem, opt='aspect',unit='radians')
dhill<-hillShade(dslope,daspect, angle=10)
map<-mapview(dem,col.regions=pal1, legend=FALSE, na.color="transparent") +
mapview(dhill,col.regions=pal2,na.color="transparent", legend=FALSE)
map@map <- map@map %>%
addFullscreenControl() map