In dgolicher/giscourse: Some code snippets wrapped up as functions

knitr::opts_chunk$set(
  collapse = TRUE, message = FALSE,warning = FALSE,
  comment = "#>"
)

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()
arne_lcover<-landcover("Arne, Dorset") 
arne_sssi<-sssi("Arne, Dorset") 
arne_phabitat<-phabitat("Arne, Dorset") 
arne_osm<-osm("Arne, Dorset")
arne_dem<-mkdem("Arne, Dorset")

Exploratory mapview

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

mapviewOptions(legend=FALSE)
mapview(arne_lcover,zcol="bhab",burst=TRUE,  hide = TRUE) %>% extras()

arne_dem[arne_dem<1]<-NA
sloper<-terrain(arne_dem,"slope","radians")
aspectr<-terrain(arne_dem,"aspect","radians")
hs<-hillShade(sloper,aspectr)

library(tmap)
arne_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

The 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)
str(lcover)
lcover %>% data.frame() %>% dplyr::group_by(bhab) %>% dplyr::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, 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 ..

corfe_dem<-mkdem(-2.057,50.64,dist=1000,z=14)
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()