Loading layers from the data base

The giscourse package has four composite high level functions that load in data for a site.

You can either use a string of characters to geocode the site, or the coordinates as latitude and longitude. If a distance in meters is provided the functions download all features within this radius. If no distance is given then the default is to download all the features inside the OS 5km grid square for the site.

Example 1: Geocode Arne and get 5km grid square

library(giscourse)
conn<-connect()
lcover<-landcover("Arne, Dorset") 
sssi<-sssi("Arne, Dorset") 
phabitat<-phabitat("Arne, Dorset") 
osm<-osm("Arne, Dorset")
dem<-mkdem("Arne, Dorset")

This code has loaded land-cover, sssi boundaries, priority habitat, open street map paths, tracks and roads and the 12m resolution digital elevation model for the 5km grid square that contains R.

Higher resolution DTMs are available by changing the zoom value. The highest resolution provided by elevatr is 3.5 meters. The server contains some Lidar data at higher resolution for chosen sites such as R. The functions given here will work for any site in the South of England.

Exploratory mapview

Use burst =TRUE to plot each land use separately. Useful for exploration.

mapviewOptions(legend=FALSE)
mapview(lcover,zcol="bhab",burst=TRUE,  hide = TRUE) %>% extras()
mapview(osm) +mapview(dem)

Terrain analysis

dem[dem<1]<-NA
sloper<-terrain(dem,"slope","radians")
aspectr<-terrain(dem,"aspect","radians")
hs<-hillShade(sloper,aspectr)
library(tmap)
lcover %>% filter(bhab !="Saltwater") %>% 

qtm("bhab") + tmap_style("classic") + tmap_options(legend.outside=TRUE) + tm_scale_bar() + tm_compass(type="arrow") + tm_shape(hs) +
 tm_raster(alpha=0.5,palette = gray(0:100 / 100), n = 100, legend.show = FALSE)

Example 2: Distance around a point

These coordinates are for a study point in the New forest

lon<--1.5
lat<-50.86
landcover(lon,lat,dist=2000) ->lcover
sssi(lon,lat,dist=2000) ->sssi
phabitat(lon,lat,dist=2000) -> phabitat
mapview(lcover,zcol="bhab",burst=TRUE,  hide = TRUE) %>% extras()

lcover %>% 

qtm("bhab") + tmap_style("classic") + tmap_options(legend.outside=TRUE) + tm_scale_bar() + tm_compass(type="arrow")

It’s easy to calculate total areas of each habitat.


lcover$bhab<-as.factor(lcover$bhab)

lcover %>% data.frame() %>% group_by(bhab) %>% summarise(n=n(), area_ha=round(sum(area)/10000,1)) -> areas
library(ggplot2)

areas %>% arrange(area_ha) %>%  mutate(bhab = factor(bhab, bhab)) %>%
ggplot(aes(x=bhab,y=area_ha,lab=area_ha)) +geom_bar(stat="identity",fill="grey") + geom_label(aes(x = bhab, y = area_ha, label = paste(area_ha,"ha"))) + coord_flip() + theme_bw()

Saving the results to a geodata file

All the layers can be written to a single geodata package file. The geodata package can then be exported from the server and dropped into a local QGIS project for work offline.

st_write(arne_lcover, dsn="geodata.gpkg", layer="arne_lcover")
st_write(arne_sssi, dsn="geodata.gpkg", layer="arne_sssi")
## Etc, etc ..

Insolation and terrain wetness analysis

The package allows the SAGA functions to be run on a clipped digital terrain model.

corfe_dem<-mkdem(-2.057,50.64,dist=1000,z=13)
sol<-insol(corfe_dem)
twi<-twi(corfe_dem)
mapview(twi)
sloper<-terrain(corfe_dem,"slope","radians")
aspectr<-terrain(corfe_dem,"aspect","radians")
hs<-hillShade(sloper,aspectr)
landcover(-2.057,50.64,dist=1000)->lc

tm_shape(hs) +
 tm_raster(palette = gray(0:100 / 100), n = 100, legend.show = FALSE) + tm_shape(lc) + tm_fill("bhab",alpha=0.8) +tmap_style("classic") +  tmap_options(legend.outside=TRUE) ->mp

tmap_mode("view")
mp + tm_minimap()

Deforestation analysis

We can use the get_hansen function to cut out raw data for any reasonably sized study area from Hansen’s global mosaic. For example, to obtain a 4 km by 4km grid square use a distance of 2km around a point showing a high rate of deforestation.


fr<-get_hansen(x=101,y=0.01, dist=2)

The loss layer holds the year at which the loss was detected. The WMS loss layer simply shows all loss during the period from 2001 to 2018 in red.

fr$loss[fr$loss<1]<-NA
mapview(fr$loss,legend=TRUE) %>% hansen_wms()
d<-data.frame(table(round(getValues(fr$loss),0))) # Tabulate the results
pix_area<-28*28
cell_area<- 40*40
d %>% 
  mutate(area_ha=Freq*pix_area/10000) %>%
  mutate(percent=round(100*area_ha/cell_area,1)) %>%
  mutate(year=as.numeric(as.character(Var1))+2000) ->d
DT::datatable(d)

So we can see that the deforestation intensified in 2008.

d %>% ggplot(aes(x=year,y=percent)) +geom_line() +ylab("Annual deforestation as percentage of total area")

The pattern of loss can be calculated easily by forming a raster stack containing the forest remaining for each year.

forest<-fr$forest2000>50
ls<-function(yr=1){
 forest * (1 -round((fr$loss<yr)*1,0))
}

list<-7:18

loss<-raster::stack(lapply(list,ls))
names(loss)<-paste("Forest", 2007:2018)
plot(loss,legend=FALSE)