library(giscourse)
library(dplyr)
library(stringr)
library(elevatr)
conn<-connect()
mapviewOptions(basemaps =c("OpenStreetMap",
"Stamen.TonerLite", "Esri.WorldImagery", 
"Esri.WorldShadedRelief", 
"Stamen.Terrain"))

Data sources

The data from the following sources have been transfered to the local postgis data base for easy access.

Land cover 2015: https://www.ceh.ac.uk/services/land-cover-map-2015

Priority habitat inventory: https://naturalengland-defra.opendata.arcgis.com/datasets/priority-habitat-inventory-south-england

Sites of special Scientific interest: https://data.gov.uk/dataset/5b632bd7-9838-4ef2-9101-ea9384421b0d/sites-of-special-scientific-interest-england

GBIF: https://www.gbif.org/en/

Loading the data

The data are loaded into R through querying postgis to find all the features within a 10 km radius of Corfe Common.

corfe_common_sssi<-giscourse::sssi(x=-2.0595,y=50.6337,dist=1000)
phabitat<-giscourse::phabitat(x=-2.0595,y=50.6337,dist=10000)
lcover<-giscourse::landcover(x=-2.0595,y=50.6337,dist=10000)
write_sf(corfe_common_sssi,"data/Corfe_common.shp")
write_sf(phabitat,"data/phabitat.shp")
write_sf(lcover,"data/lcover.shp")

National distributions

The data from GBIF counted by 5km grid cell.

dd<-st_read(conn,"bees_gbif_grid")

dd %>% filter(species =='Succisa pratensis') -> devils_bit
dd %>% filter(species =='Knautia arvensis') ->  field_scabious
dd %>% filter(species =='Scabiosa columbaria') -> small_scabious
dd %>% filter(species =='Andrena marginata') -> andrena_marginata
dd %>% filter(species =='Nomada argentata') -> nomada_argentata


mapview(andrena_marginata,zcol="count",legend=FALSE) +
mapview(nomada_argentata,zcol="count",legend=FALSE) +mapview(devils_bit,zcol="count",legend=FALSE) +
mapview(field_scabious,zcol="count",legend=FALSE) -> mp

mp@map  %>% leaflet.extras::addFullscreenControl()

Map centred around Corfe common

devils_bit_scabious<-st_read(conn,query="select * from bees_gbif where species = 'Succisa pratensis'")
andrena_marginata<-st_read(conn,query="select * from bees_gbif where species = 'Andrena marginata'")
devils_bit_scabious<-st_crop(devils_bit_scabious,st_bbox(lcover))
andrena_marginata<-st_crop(andrena_marginata,st_bbox(lcover))

library(stringr)
phabitat %>% filter (str_detect(main_habit, 'grass')) -> priority_grassland
lcover %>% filter (str_detect(bhab, 'grass')) -> all_grassland

mapview(corfe_common_sssi,col.regions="darkgreen") + mapview(devils_bit_scabious,zcol="year", legend=FALSE) +
  mapview(andrena_marginata,zcol="year", legend=FALSE) +
  mapview(priority_grassland,zcol="main_habit",burst=TRUE,legend=FALSE) + mapview(all_grassland,zcol='bhab',legend=FALSE) -> mp
mp@map %>% leaflet.extras::addFullscreenControl() %>% 
  addMeasure(primaryLengthUnit ='meters', secondaryLengthUnit='kilometers',primaryAreaUnit='hectares', secondaryAreaUnit='acres',position="topright")

Lidar data

# path<- "~/geoserver/data_dir/rasters/lidar/dsm/composite_1m"
# system(sprintf("gdalbuildvrt dsm1m.vrt %s/*.asc",path))
# path<- "~/geoserver/data_dir/rasters/lidar/dtm/composite_1m"
# system(sprintf("gdalbuildvrt dtm1m.vrt %s/*.asc",path))

There is a problem with the openly available Lidar data. It does not cover the whole of Corfe Common.

library(raster)
dsm<-raster("dsm1m.vrt")
crs(dsm) <- st_crs(corfe_common_sssi)$proj4string
dsm<-raster::crop(dsm,corfe_common_sssi)
mapview(dsm) + mapview(corfe_common_sssi)

Seamless open digital elevation model

The same problem affects the data obtained using the elevatr package in R. This is because the same lidar data source is being used.

dem <- get_elev_raster(corfe_common_sssi, z = 14,clip="bbox")

contours<-rasterToContour(dem,nlevels=20)
sloper<-terrain(dem, opt='slope', unit='radians')
aspectr<-terrain(dem, opt='aspect',unit='radians')
hillshade<-hillShade(sloper,aspectr,angle=5)
  plot(hillshade, col = gray(0:100 / 100), legend = FALSE)
plot(dem, col = terrain.colors(25), alpha = 0.2, legend = TRUE, add = TRUE)

require(RColorBrewer)
pal<-terrain.colors(100)
pal <- colorRampPalette(pal)


mapview(dem,col.regions=pal, legend=FALSE) + mapview(contours,legend=FALSE)
# ## Write out the raster in saga format first.
# writeRaster(dem,"data/saga_dem.sgrd",format="SAGA", overwrite=TRUE)
# 
# ## Run the commands (WILL ONLY WORK IF THE LIBRARIES ARE INSTALLED) 
# 
# system ("saga_cmd ta_lighting 2 -GRD_DEM data/saga_dem.sgrd  -DAY '13/08/19' -GRD_DIRECT   data/saga_sol")
# system ("saga_cmd ta_hydrology 15 -DEM  data/saga_dem.sgrd  -TWI   data/saga_TWI")
# 
# ## Read the results back to R
# 
# sol<-raster("data/saga_sol.sdat")
# twi<-raster("data/saga_TWI.sdat")
# 
# ## Re set the coordinate reference systems (which get lost in the process of translation to SAGA and back)
# sol@crs<-dem@crs
# twi@crs<-dem@crs
# 
# ## Write out the results 
# writeRaster(sol,"data/insol.tif",format="GTiff", overwrite=TRUE)
# writeRaster(twi,"data/twi.tif",format="GTiff", overwrite=TRUE)
# 
# ## Remove the intermediate SAGA files
# system("rm data/saga*")

Insolation

sol<-raster("data/insol.tif")
sol[sol<3]<-NA
mapview(sol) 
# system ("zip -r data.zip . data")

Download the data files for ARC GIS.

http://r.bournemouth.ac.uk:82/examples/corfe_common/data/data.zip