---
title: "Bird habitat use analysis"
author: "Duncan Golicher"
date: "11/24/2023"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message=FALSE, warning = FALSE)
```

```{r}
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

```{r}
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.


```{r}
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)

```{r}
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.

```{r}
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)

```{r}
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


```{r}
urban$LandCover<-"Urban"
```

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

```{r}
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.



```{r}
#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.

```{r}
habitat2<-rbind(rural,urban)
birds %>% st_join(habitat2) ->birds
```



```{r}
mapview(habitat2, zcol="LandCover") + mapview(birds)
```


```{r}
birds %>% st_drop_geometry() -> birds
```

```{r}
pivot_wider(birds, names_from = Species,id_cols = LandCover, values_from = Species, values_fn=length) ->bird_mat
bird_mat
```

