library(tidyverse)
library(aqm)
library(giscourse)
library(sf)
library(mapview)
library(raster)
library(leaflet.extras)

Introduction

In this session you wil practice using some of the GIS capabilities of R in the context of the data collected at Arne. If you have already analysed the pine growth data and the pine density data with respect to the derived variables there is not too much to add in order to complete the assignment. However the way the variables that might be related to the pine densities was actually obtained may still be rather mysterious. This class will show you how you can work with lidar derived digital elevation rasters (DEMs and DTMs)

Obtaining Lidar data

The repository for openly available lidar data online has changed several times in the last few years. The most comprehensive source is now held by DEFRA.

https://environment.data.gov.uk/DefraDataDownload/?Mode=survey

Downloading the data is fairly straightforward, but is still something of a hassle for use in class.The DEFRA interface can produce a large number of files for a single site. You should be aware that this data is available and can be imported into R if you are working at a specific site for your own project, but for the purposes of this unit we will use data that is already loaded. This data is a little out of date, but corresponds to the situation when the field data was collected.

Quadrats

The pine density data was obtained by students through the use of circular quadrats. The quadrats were georeferenced by matching the ids to waypoints from a GPS. Pine densities were calculated by dividing the number of pines found in the circular areas by the area of the circles.

data("arne_pines")
d<-arne_pines
str(d)
## 'data.frame':    198 obs. of  6 variables:
##  $ lon         : num  -2.04 -2.04 -2.04 -2.04 -2.05 ...
##  $ lat         : num  50.7 50.7 50.7 50.7 50.7 ...
##  $ Year        : num  2019 2019 2019 2019 2019 ...
##  $ site        : chr  "Heath" "Heath" "Heath" "Heath" ...
##  $ npines      : num  11 7 12 3 28 18 51 54 4 8 ...
##  $ pine_density: num  0.875 0.557 0.955 0.239 2.228 ...

The actual areas varied between years as a result of the students using different lengths of rope to measure out the circles. So use the densities, not the raw counts. Notice that these data simply form a data frame. They could be loaded into R from Excel as a CSV.

Pine densities in heath and restauration area.

The data is not very representative of the two areas due to very limited sampling. However you can apply a simple t-test based approach to look at any differences between densities in the two areas. Strong caveats apply.

g0<- ggplot(d,aes(x=site,y=pine_density)) 
g0+geom_boxplot()

aqm::ci(g0)

What would a t test show?

Making a spatial object

In order to turn them into a spatial object we need to tell R something about them.

We know that the x coordinates are the longitude (lon) and the y coordinates are the latitude (lat). These are not projected so the coordinate reference system (CRS) is 4326. If the points had been collected in OS national grid form then the CRS would be 27700.

d<-st_as_sf(arne_pines,coords = c("lon","lat"),crs=4326)

Now we can map the positions. If you want a full screen control this is added using the leaflet.extras package. By making year the zcol and “bursting” the data the map will show a layer for each year and a legend if legend is true, or without in this case.

map<-mapview(d, zcol="Year", burst=TRUE, legend=FALSE)
map@map %>% leaflet.extras::addFullscreenControl()

Adding a minimap also uses the leaflet.extras package. The lower the negative value for the zoomLevelOffset the greater the difference in scale between the two maps. You can also select different base maps. This can be useful for showing the position of a study site.

map@map %>% addFullscreenControl() %>% 
  addMiniMap(position = "bottomleft",tiles="OpenStreetMap",zoomLevelOffset = -6, toggleDisplay=TRUE)