02-SpatialSubsets.Rmd

---
title: "02-SpatialSubsets"
author: "Dimitrios Markou"
date: "`r Sys.Date()`"
output: html_document
---

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

# Chapter 2: Spatial Subsetting: KBA's and Priority Places

##### Authors: Dimitrios Markou, Danielle Ethier

| In [**Chapter 1: Spatial Data Exploration**](01-SpatialDataExploration.Rmd), you distinguished between spatial data types (vector and raster) and explored NatureCounts data using a spatio-temporal map. Now you are ready to focus your data exploration to a specific geographic area. Analyzing NatureCounts data within specific boundaries is relevant to many research applications. Here we use Key Biodiversity Areas (KBA) and Priority Places as examples of research areas of interest.

# 2.0 Learning Objectives

By the end of **Chapter 2 - Spatial Subsetting**, users will know how to:

-   Import and map polygon boundary files based on attributes
-   Reproject NatureCounts data to the coordinate reference system of your spatial layer
-   Spatially filter and map NatureCounts data for an area of interest
-   Read, process, and visualize spatial vector data within areas of significant conservation potential: [Key Biodiversity Areas (KBAs)](https://kbacanada.org/about/) and [Priority Places for Species at Risk](https://environmental-maps.canada.ca/CWS_Storylines/index-ca-en.html#/en/priority_places-lieux_prioritaires).

The data used in this tutorial are downloaded from NatureCounts, the [KBA Canada Map Viewer](https://kbacanada.org/explore/map-viewer/), and the [Priority Places - Open Government Portal](https://open.canada.ca/data/en/dataset/91219d24-e877-4c8a-8bd2-b2b662e573e0). You will save the downloaded shapefiles to your working directory.

To view your working directory.

```{r}
getwd()
```

This R tutorial requires the following **packages**:

```{r, message = FALSE}
library(naturecounts)
library(sf)
library(tidyverse)
library(mapview)
library(leaflet)
library(leaflet.extras)
```

# 2.1 Key Biodiversity Areas

#### *Example 2a* - You are interested in assessing the spatial distribution of Wood Ducks (Ontario Breeding Bird Atlas data) across the KBAs found within the province of Ontario.

Navigate the the [KBA Canada Map Viewer](https://kbacanada.org/explore/map-viewer/) and filter the data for Ontario using the left hand `Province/Territory` filter. Then select 'Download'. You will want to select both `csv` and `shp` for this example.

We can read in our KBA polygons using the `sf` package once it is in your working directory. The downloaded files was renamed for this example so you will need to change the code to match your file name.

```{r}
ontario_kba <- sf::st_read("ontario_kba.shp")
```

`sf` objects are stored in R as a spatial dataframe which contains the attribute table of the vector along with the geometry type. When we examine the dataframe, it looks like there are many duplicate entries including duplicate geometries (vertices). To clean this up, we can apply the `st_make_valid()` and `distinct()` functions to our spatial dataframe:

```{r}
ontario_kba <- ontario_kba %>% st_make_valid() %>% distinct()
```

Our spatial data is also accompanied by a CSV file that contains additional useful attributes (landcover, species, etc) concerning our KBAs. Let's read in the accompanying CSV file for our KBA layer.

```{r}
kba_attributes <- read.csv("ontario_kba.csv")
```

Great! We can now join these dataframes using the handy `tidyverse` package. However, we'll want to select for specific columns first to avoid redundancies before performing our join:

```{r}
kba_attributes <- kba_attributes %>%
  select("SiteCode",
         "DateAssessed",
         "PercentProtected",
         "BoundaryGeneralized",
         "Level",
         "CriteriaMet",
         "ConservationActions",
         "Landcover",
         "Province",
         "Species")
```

Both dataframes now contain unique columns, after our selection. We apply the `full_join()` function to hold all attributes within one dataframe.

```{r}
ontario_kba <- full_join(ontario_kba, kba_attributes, by = "SiteCode")
```

We can visualize the **ontario_kba** data with an interactive map, using the `leaflet` package.

```{r}
leaflet(width = "100%") %>%
  addTiles() %>%
  addPolygons(data=ontario_kba,color = "black", weight = 2, smoothFactor = 1,opacity = 1.0, fillOpacity = 0.5, fillColor = "red") %>% addFullscreenControl() %>%
  addLegend(colors = c("red"),labels = c("Ontario KBAs"),position = "bottomright")
```

Similarly, the package `mapview` (based on leaflet) can also be used to make interactive plots. We can represent specific attributes like so:

```{r}
mapview::mapview(ontario_kba, zcol = "PercentProtected")
```

In this example, were interested in all the KBA polygons of Ontario. However, if you were working with a larger data set, it is possible to filter your dataframe to retrieve only specific polygons that meet certain criteria relevant to your research. To do so, we can apply filters based on a variable condition. For example, say we only wanted KBA's greater than 100km\^2 in size:

> To execute this code chunk, remove the \#

```{r}
# kba_name <- ontario_kba %>% 
  # dplyr::filter(Area > 100) 
```

Let's search NatureCounts for the Ontario Breeding Bird Atlas point count dataset using `meta_collections()` and the Wood Duck species ID using `search_species()`.

```{r}
collections <- meta_collections()
View(meta_collections())
```

```{r}
search_species("wood duck")
```

Now we can download the NatureCounts data. Remember to change `testuser` to your personal username.

```{r}
atlas_on <- nc_data_dl(collections = "OBBA2PC", species = 360, username = "testuser", info = "spatial_data_tutorial")
```

We can then convert our NatureCounts data into a spatial object. To do so, we deploy the `st_as_sf` function and specify the coordinate reference system (CRS).

The CRS of our KBA sf object can be returned with `st_crs()`.

```{r}
sf::st_crs(ontario_kba)
```

Our KBA sf object is stored with World Geodetic System 1984 (WGS 84) coordinates, EPSG = 4326. Now we can convert our **atlas_on** dataframe to an sf object using the same CRS.

```{r}
atlas_on_sf <- sf::st_as_sf(atlas_on,
                        coords = c("longitude", "latitude"), crs = 4326)
```

Now let's ensure that the conversion was successful. You'll notice a new geometry column where each observation is a point.

```{r}
str(atlas_on_sf) # view the sf object
```

The `st_transform()` function can be applied to project our spatial object using a different CRS like NAD83 / UTM zone 16N (EPSG = 26916).

```{r}
ontario_kba <- sf::st_transform(ontario_kba, crs = 26916) 
```

It can also be used to ensure that the CRS of our spatial objects match.

```{r}
atlas_on_sf <- st_transform(atlas_on_sf, crs = st_crs(ontario_kba))
```

To identify Wood Ducks observed across Ontario KBA's, we can apply the `st_intersection()` function.

```{r}
wood_ducks_kba <- sf::st_intersection(ontario_kba, atlas_on_sf)
```

You will get a `warning` message, which you can safely ignore.

```{r}
# If need be, transform your spatial data back to EPSG:4326 to visualize with leaflet
ontario_kba <- st_transform(ontario_kba, crs = 4326)
wood_ducks_kba <- st_transform(wood_ducks_kba, crs = 4326)
```

Apply leaflet to visualize our polygon data and point data, using the `addPolygons()` `addCircleMarkers()` arguments.

```{r}
leaflet(width = "100%") %>%
  addTiles() %>%
  addPolygons(data = ontario_kba, color = "black", weight = 2, smoothFactor = 1, 
  opacity = 1.0, fillOpacity = 0.5, fillColor = "violet") %>%
  addCircleMarkers(data = wood_ducks_kba, radius = 3, color = "orange", 
                   stroke = FALSE, fillOpacity = 0.8) %>%
  addFullscreenControl() %>%
  addLegend(colors = c("violet", "orange"), labels = c("Ontario KBA", "Wood Duck Observations"), position = "bottomright")
```

After geoprocessing our data, we can write out any sf objects to Shapefiles on a disk, where the argument delete_layer = TRUE is used to overwrite an existing file.

> To execute this code chunk, remove the \#

```{r}
# st_write(wood_ducks_kba,"wood_ducks_kba.shp", driver = "ESRI Shapefile", delete_layer = TRUE)
```

# 2.2 Priority Places

#### *Example 2b* - You want to assess the spatial distribution of Wood Ducks (Ontario Breeding Bird Atlas data) across the Long Point Walsingham Forest Priority Place.

Navigate to [Priority Places - Open Government Portal](https://open.canada.ca/data/en/dataset/91219d24-e877-4c8a-8bd2-b2b662e573e0). Scroll down to the **Data and Resources** section and select the Priority Places file labeled `English`, `dataset`, and `FGDB/GDB`.

First, lets create a path to our downloaded Priority Place file.

```{r}
gdb_path <- "PriorityPlaces.gdb"
```

Then, let's inspect the spatial data.

```{r}
gdb_layers <- st_layers(gdb_path)
print(gdb_layers)
```

To read in our spatial data object, we apply the `st_read` function and specify our desired data layer.

```{r}
priori_place_polygons <- st_read(dsn = gdb_path, layer = "PriorityPlacesBoundary")
```

Were interested in the spatial distribution of Wood Ducks across the Long Point Walsingham Forest Priority Place. We'll filter based on a variable condition.

```{r}
long_point_polygon <- priori_place_polygons %>%
  filter(Name == "Long Point Walsingham Forest") # filters based on multipolygon name 
```

Then reproject the Wood Duck data to match our Priority Place using `st_transform`.

```{r}
atlas_on_sf <- st_transform(atlas_on_sf, crs = st_crs(long_point_polygon))
```

Next, we'll apply our geoprocessing function to find the Wood Duck observations that intersect with our chosen Priority Place.

```{r}
wood_ducks_longpoint <- sf::st_intersection(long_point_polygon, atlas_on_sf)
```

Finally, we'll transform our spatial objects one more time before visualizing them with `leaflet`.

```{r}
long_point_polygon <- st_transform(long_point_polygon, crs = 4326)
wood_ducks_longpoint <- st_transform(wood_ducks_longpoint, crs = 4326)
```

Apply leaflet to visualize our polygon data and point data, using the `addPolygons()` `addCircleMarkers()` arguments.

```{r}
leaflet(width = "100%") %>%
  addTiles() %>%
  addPolygons(data = long_point_polygon, color = "black", weight = 2, smoothFactor = 1, 
  opacity = 1.0, fillOpacity = 0.5, fillColor = "red") %>%
  addCircleMarkers(data = wood_ducks_longpoint, radius = 5, color = "green", 
                   stroke = FALSE, fillOpacity = 0.8) %>%
  addFullscreenControl() %>%
  addLegend(colors = c("red", "green"), labels = c("Long Point Walsingham Forest", "Wood Duck Observations"), position = "bottomright")
```

After geoprocessing our data, we can write out any sf objects to Shapefiles on a disk, where the argument delete_layer = TRUE is used to overwrite an existing file.

> To execute this code chunk, remove the \#

```{r}
# st_write(wood_ducks_longpoint,"wood_ducks_longpoint.shp", driver = "ESRI Shapefile", delete_layer = TRUE)
```

**Congratulations**! You completed **Chapter 2 - Spatial Subsetting: KBA's and Priority Places**. Here, you spatially filtered and visualized NatureCounts data. In [Chapter 3](03-ClimateData.Rmd), you can explore more spatial data visualization while linking climate data to NatureCounts observations.