library(tidyverse)
library(mapview)
library(sf)
library(tmap)
library(RColorBrewer)
load("DEM.rda")
mapviewOptions(basemaps =c("Esri.WorldImagery",
                             "OpenStreetMap",
                             "Esri.WorldStreetMap"), homebutton = FALSE)
pal<- brewer.pal(12,"Paired")

Finding the distance from the coast

The first step is to find a layer that represents the caostline The rnaturalearthhires package has a poolygon of the United Kingdom

library(rnaturalearth)
coast<-st_as_sf(rnaturalearthhires::countries10)
filter(coast, NAME=="United Kingdom") -> coast

As the result of running the chunk above produces a polygon, this has to be changed to a linestring.

coast <- st_cast(coast,"MULTILINESTRING")

The data are in latitude and longitude. So to measure the distance it needs to be transformed to British national grid (27700)

coast<-st_transform(coast,27700)
birds$coast<-st_distance(birds, coast)

Finally we can drop the units and divide by 1000 to get the distance in Km.

library(units)
birds$coast<-drop_units(birds$coast/1000)

Aggregating suburban and urban Land classes.

This is one of those operations for which there are multiple options. The simplest way might be to just rename suburban areas as urban. However this would not join the underlying geometries.

One approach is to filter the data to form two new sets. The first only contains the urban and suburban classes. The second contains the rest (! %in% reads as NOT in)

habitat %>% filter(LandCover %in% c("Urban","Suburban")) -> urban
habitat %>% filter(!LandCover %in% c("Urban","Suburban")) -> rural

Now call all the Landcover in the urban data "Urban

urban$LandCover<-"Urban"

To join the geometries we need to group the data and use at_union.

urban %>% group_by(LandCover) %>% summarise(geometry=st_union(geometry), Area=sum(Area)) ->urban

We could do the same with the rural data. However the data set is very large, and the operation does take around ten minutes to complete. If it is not run (by uncommenting the lines below) then the Landclasses stay the same.

#rural$LandCover<-"Rural"
#rural %>% group_by(LandCover) %>% summarise(geometry=st_union(geometry), Area=sum(Area)) ->rural

Now to reunite the two sets of data into a new set we can use rbind.

habitat2<-rbind(rural,urban)
birds %>% st_join(habitat2) ->birds
mapview(habitat2, zcol="LandCover") + mapview(birds)
birds %>% st_drop_geometry() -> birds
pivot_wider(birds, names_from = Species,id_cols = LandCover, values_from = Species, values_fn=length) ->bird_mat
bird_mat
## # A tibble: 8 × 9
##   LandCover  Bullfinch Chaffinch Goldfinch Greenfinch `House Sparrow` Linnet
##   <chr>          <int>     <int>     <int>      <int>           <int>  <int>
## 1 Arable            30       113        98         53              69     51
## 2 Coastal           56       262       289        211             196    239
## 3 Freshwater        59       231       200        142              33     78
## 4 Grassland        290      1077       944        670             877    375
## 5 Heathland        205       740       737        571             312    529
## 6 Saltwater         10        35        22         37              16     22
## 7 Urban            267       572       788        706             868     68
## 8 Woodland         154       576       344        203             112    252
## # … with 2 more variables: Reed Bunting <int>, Yellowhammer <int>