Bird habitat use analysis 2

Introduction

In this handout we will see how to easily get a digital elevation (raster) layer for the area of interest and visualise it using mapview.

First load the libraries

library(tidyverse)
library(mapview)
library(sf)
library(tmap)
library(RColorBrewer)
library(elevatr)
library(raster)
library(leaflet.extras)

Now load the birds data as before.

birds<- st_read("Bird-data-v1.shp")

## Reading layer `Bird-data-v1' from data source 
##   `/home/rstudio/webpages/books/Bird_mapping_assignment/Bird-data-v1.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 16732 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 378500 ymin: 75500.02 xmax: 435500 ymax: 113000.1
## Projected CRS: OSGB 1936 / British National Grid

We will use some palletes later on which are defined here.

pal<- brewer.pal(12,"Paired")
pal1<-c("darkblue",terrain.colors(100))
pal2<-gray(0:100 / 100)

Using the elevatr package

The elevatr package that was loaded up with the other packages for turning R into a GIS is extremely useful. It pulls elevation data down from the web in the form of a raster layer. The boundaries of the raster match those of a vector layer given as input. The “z” parameter defines the resolution of the digital elevation model that is returned. Be careful not to set the z value too low (high resolution is low z) as the layer will be too big too handle. A z vaue of 10 is about right here.

dem<-elevatr::get_elev_raster(birds, z=10)

If you just print the dem you will get basic statistics on the raster layer.

print(dem)

## class      : RasterLayer 
## dimensions : 1550, 2075, 3216250  (nrow, ncol, ncell)
## resolution : 48.27604, 48.27604  (x, y)
## extent     : 367259.9, 467432.7, 43263.96, 118091.8  (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     : file2f795b91e39c.tif 
## names      : file2f795b91e39c 
## values     : -32768, 32767  (min, max)

Plotting the dem produces a simple static map.

plot(dem)

This might look better

dem[dem<0] <-0
plot(dem, col=pal1)

Making a hillshaded image

The raster package has a large number of functions for working with DEMs (more are also available through linking R to GRASS and SAGA).

The hillshade function requires slope and aspect layers in radians as its input.

dslope<-terrain(dem, opt='slope', unit='radians')
daspect<-terrain(dem, opt='aspect',unit='radians')
dhill<-hillShade(dslope,daspect, angle=10)

plot(dhill, col=pal2)

Quick extraction of values

Its very simple in R to extract the values that points fall on over a raster layer and add them to the point data.

birds$elevation<-raster::extract(dem,birds)

Showing the DEM as a mapview

The dem and the hillshaded layer can be added together and visualised using mapview.