Terrain analysis

Load the packages

There are a large number of packages used here, as they convert R into a GIS. Use this chunk whenever you are conducting spatial analysis in R.

library(sf)
library(tidyverse)
library(aqm)
library(mapview)
library(mgcv)
library(raster)
library(tmap)
library(giscourse)

Get the data

The data in this case are already loaded on the server within the aqm package. The full data includes quadrats from the “restoration” site that we will not use here.

data("arne_pines")
arne_pines %>% filter(site=="Heath" ) -> arne_pines
quadrats<-st_as_sf(arne_pines,coords = c("lon","lat"),crs=4326)

Load spatial data

Again, the spatial data is on the server. If you were assembling the layers from data used in a GIS there would be more steps needed to load each layer.

data(arne_lidar)

Plot the digital terran model

A simple, but effective map can be made by first converting the DTM into contours then plotting the results.

contours<-rasterToContour(dtm, levels=seq(0,60,5)) ## Note the sequence
qtm(dtm) + qtm(contours) +qtm(quadrats) ## Qtm is a quick way of forming a thematic map with deafault settings

Study area boundaries

The study area is included in the loaded data as a polygon.

mapview(study_area,alpha.regions = 0,lwd=4)

Cropping

A very common GIS operation is to crop a large raster layer to the bounding box of a smaller area. This is very easy in R. We can use the study area polygon for cropping.

pol<-as(study_area,"Spatial") ## Convert to the older spatial classes
dsm<-raster::crop(dsm,pol)
dtm<-raster::crop(dtm,pol)

Mapping results

We can now build a mapview to look at the cropped dsm and dtm quite easily. The default colours don’t look very good, so the code also includes a palette.

contours<-rasterToContour(dtm, nlevels=10)
map2<- mapview(dsm,col.regions=terrain.colors(1000),legend=TRUE) + mapview(dtm,col.regions=terrain.colors(1000)) + mapview(contours) + mapview(quadrats)
map2 %>% extras()
## Loading required package: leaflet.extras
## Loading required package: leaflet

Slope and aspect

There are a large number of algorithms built into R that can be used to conduct terrain analysis. These produce the same results as the comparable options in QGIS. Working in R has the advantage of being quicker once a script has been built, and also reproducible.

For example deriving slope and aspect in degrees just involve running the following lines. Really quick and easy.

slope<-terrain(dtm, opt='slope', unit='degrees')
aspect<-terrain(dtm, opt='aspect',unit='degrees')

A quick call to plot shows the results with a default palette.

plot(slope)

Same for aspect.

plot(aspect)

Hillshade

Making a hillshade layer requires slope and aspect to be calculated in radians, rather than degrees

sloper<-terrain(dtm, opt='slope', unit='radians')
aspectr<-terrain(dtm, opt='aspect',unit='radians')
hillshade<-hillShade(sloper,aspectr,angle=5)
plot(hillshade,col=grey.colors(1000),legend=FALSE)

The dtm can be draped over the plot by making it semi transparent by setting an alpha level.

plot(hillshade,col=grey.colors(1000),legend=FALSE) 
plot(dtm, add=TRUE, col=terrain.colors(1000),alpha=0.4)

Insolation and topographic wetness index

R can also call some more advanced algorithms, such as those which calculate topographic wetness index and insolation. These algorithms effectively combine elements such as slope and aspect into more interpretable variables.

twi<-giscourse::twi(dtm)
sol<- giscourse::insol(dtm, day = "11/06/2030")

The topographic wetness index measures the upslope area draining through a certain point per unit contour length. In other words values of the topographic wetmess index (which is unitless) can be used to compare where a quadrat lies on a slope. Water drains from low values to high values, so on average the soil will be moister although we may not know by how much) when the index is low.

plot(twi)

The insolation layer is calculated precisely using mechanistic preinciples. When a date is provided for the algorithm the routine steps through each hour of the day and calculates the solar energy that each pixel will receive on a totally clear day through direct beam insolation. All other things being equal, pixels with high values will warm up more than pixels with low values on a sunny day. The pattern will change over the year.

plot(sol)

Raster algebra

Carrying out raster algebra (i.e. sums) operations in R is very simple and inuitive.

To obtain the height of vegetation (and maybe some buildings ) we just need to subract the digital terrain model from the digital surface model. This is a canopy height model (chm)

chm<-dsm-dtm
plot(chm)

## Adding the values to the quadrats

In the context of the assignment the reason for producing the raster maps is to find the values for the variables that coincide with the quadrats that we placed on the ground. This is like placing your finger on each of the maps at the position of the quadrat, reading off the value, then adding the value to the original spreadheet of results. Rather than go through this very tedious process we can run some lines of R code to do the same.

q<-as(quadrats,"Spatial")
quadrats$twi<-raster::extract(twi,q)
quadrats$sol<-raster::extract(sol,q)
quadrats$dtm<-raster::extract(dtm,q)
quadrats$slope<-raster::extract(slope,q)
quadrats$chm <-raster::extract(chm,q)

Distances to trees

To calculate the distance between the quadrats and the nearest tree we can use the canopy height model. Assuming that the heights represent vegetation we can first make polygons that represent groups of pixels over a given height.

chm[chm<2]<-NA ## Set below 2m as blanks (NA)
chm[chm>=2]<-1 ## Set above 2m to a single value
trees<-rasterToPolygons(chm, dissolve=TRUE)
trees<-st_as_sf(trees)
st_crs(trees)<-27700
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform
## for that
trees<-st_cast(trees$geometry,"POLYGON")
mapview(trees)

Now once we have the objects that represent the “trees” we can make a distance matrix. This contains all the distances measured between the quadrats and the trees.

quadrats<-st_transform(quadrats,crs=27700)

distances<-st_distance(quadrats,trees)
dim(distances)
## [1]  184 1666
quadrats$min_dist<-apply(distances,1,min)
mapview(quadrats, zcol='min_dist') + mapview(trees) ->map4
map4 %>% extras()
quadrats %>% st_set_geometry(NULL) -> quads
save(quads,file="/home/rstudio/rstudio/aqm/arne_quads.rda") ## Make available for independent analysis
dt(quads)