process-landscape.Rmd

---
title: "Process Landscape Data"
subtitle: "Download and summarise landscape data at point locations"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Process Landscape Data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

This article details various ways to download and process landscape data to build predictive models of animal diversity. For example, the landscape data may be derived from (1) remotely sensed satellite imagery, (2) open-source maps, or (3) manually generated sources (e.g., from on-site mapping, design scenarios). We give examples of each in this article.

<br>

```{r out.width = "70%", fig.align='center', dpi = 300, echo = FALSE}
knitr::include_graphics("framework_process-landscape.png", dpi = 300)
```
<center><b> Figure: Broad overview of the data workflow for a chosen animal group </b></center>

<br>

First, load the necessary packages to run the analysis:

```{r load required libraries, message = FALSE, warning = FALSE, eval = FALSE}
library("biodivercity")
library("dplyr") # to process/wrangle data
library("sf") # to process vector data
library("tmap") # for visualisation
```

```{r load biodivercity while in dev, include = FALSE}
devtools::load_all()
library("sf")
library("dplyr")
library("tmap")
```


Next, define an output directory for raw and processed data to be exported to:

```{r define output directories, warning = FALSE, message = FALSE, eval = FALSE}
output_dir <- "</FILEPATH/TO/DIRECTORY>"
```

<br>

In this article, the Punggol (`PG`) area in Singapore from the example dataset `sampling_areas` will be used. Landscape data can be summarised at point locations where animal data are available (e.g., animal surveys have been conducted); to demonstrate this, we randomly generate three points within Punggol and assign it to the object '`points`':

```{r get example data ready}
data(sampling_areas)

sampling_areas <- sampling_areas %>%
   filter(area == "PG")

# randomly generate 3 points
set.seed(30)
points <- st_sample(sampling_areas, 3) %>% 
  st_as_sf() %>% 
  mutate(point_id = row.names(.)) # unique identifier
```

```{r eval = FALSE}
tmap_mode("view")
tm_basemap(c("CartoDB.Positron")) +
  tm_shape(sampling_areas) + tm_borders() +
  tm_shape(points) + tm_dots()
```


```{r plot punggol area, warning = FALSE, message = FALSE, fig.height = 1.0, dpi = 300, out.width="100%", fig.align='center', echo = FALSE}
tmap_mode("view")
tm_basemap(c("CartoDB.Positron")) +
  tm_shape(sampling_areas) + tm_borders() +
  tm_shape(points) + tm_dots()
```

<br>

---

## 1. Remotely sensed land cover

Land cover can be a useful predictor of animal diversity. For example, publicly-available data from Sentinel-2 may be downloaded using the package [`sen2r`](http://sen2r.ranghetti.info/index.html). The package includes data pre-processing steps such as cloud masking, atmospheric correction and the calculation of spectral indices.

```{r sen2r gui, eval = FALSE}
library(sen2r)
sen2r() # view GUI
```


```{r out.width = "80%", fig.align='center', dpi = 300, echo = FALSE, fig.cap = "**Figure: A screenshot of the `sen2r` graphical user interface accessible by running `sen2r::sen2r()`.**"}
knitr::include_graphics("sen2r-gui.png", dpi = 300)
```

<br>

Data processing parameters can be specified as arguments within `sen2r()`, or loaded from a JSON file. In the example below, we download data from 1 Jan to 30 Jun 2021, with not more than 15% cloud cover within `sampling_areas`. For each image, the spectral indices [`NDVI`](https://www.indexdatabase.de/db/i-single.php?id=58) and [`NDWI2`](https://www.indexdatabase.de/db/i-single.php?id=546) are processed; these will be used to map vegetation and water cover, respectively. Do note that you will need to enter your account credentials for the [ESA Sentinel Hub](https://scihub.copernicus.eu/dhus/#/home). Detailed steps and the use of other servers (e.g. Google Cloud) can be found in the [package website](http://sen2r.ranghetti.info/index.html).

```{r sen2r, eval = FALSE}

# create sub-directories
dir.create(paste0(output_dir, "/sen2r")) # for processed output files
dir.create(paste0(output_dir, "/sen2r/raw")) # for raw files

# execute workflow
out_paths <- sen2r(
  timewindow = as.Date(c("2021-01-01", "2021-06-30")),
  max_cloud_safe = 15, # max 15% cloud cover in raw tile
  max_mask = 15, # max 15% cloud cover over sampling_areas
  extent = geojsonsf::sf_geojson(sampling_areas) # needs geojson as input
  clip_on_extent = TRUE,
  list_indices = c("NDVI", "NDWI2") # to map vegetation and water
  list_rgb = "RGB432B", # output an RGB image as well
  path_l2a = paste0(output_dir, "/sen2r/raw"), 
  path_l1c = paste0(output_dir, "/sen2r/raw"),
  path_out = paste0(output_dir, "/sen2r"),
  path_indices = paste0(output_dir, "/sen2r"))
```

For each spectral index, we can combine all raster images captured within the specified date range by averaging the pixel values across files, thus forming an image mosaic. This allows us to avoid relying on any one image for our analysis, and to deal with missing data (e.g. due to high cloud cover) during the period of interest. The function `mosaic_sen2r()` creates the mosaic by averaging the pixel values across the multiple rasters. It is a convenience function that processes the images based on the output folder structure of `sen2r()`; ensure that the function argument `parent_dir` in `mosaic_sen2r()` is similar to the function argument `path_out` in `sen2r()`.

```{r mosaic generation, eval = FALSE}
mosaic_sen2r(parent_dir = paste0(output_dir, "/sen2r"))
```

<br>

At this point, we have a single image mosaic (continuous raster) for each spectral index, for a given area/period of interest. While the continuous values from these rasters may be used directly in analyses, there may be instances were we want to work with discrete classes of land cover. To do so, we can classify pixels into one of two classes (e.g. vegetated or non-vegetated; water or non-water), based on an adaptively derived threshold value. For example, we can use Otsu's thresholding ([Otsu, 1979](https://cw.fel.cvut.cz/wiki/_media/courses/a6m33bio/otsu.pdf)), which tends to outperform other techniques in terms of stability of results and processing speed, even with the presence of > 2 peaks in the histogram of pixel values ([Bouhennache et al., 2019](https://www.tandfonline.com/doi/abs/10.1080/10106049.2018.1497094)). This can be implemented using the function `threshold_otsu()`. As an example, let's load the NDVI raster for the Punggol area in Singapore for subsequent classification. The NDVI is a measure of healthy green vegetation, based on the tendency of plants to reflect NIR & absorb red light. It ranges from -1 (non-vegetated) to 1 (densely vegetated). The figure below shows the histogram of NDVI values within Punggol, as well as the adaptively derived threshold value (vertical line).

```{r}
ndvi_mosaic <- 
  system.file("extdata", "ndvi_mosaic.tif", package="biodivercity") %>% 
  terra::rast()

# get Otsu's threshold
threshold <- threshold_otsu(ndvi_mosaic)
threshold
```

```{r include = FALSE}
histogram <- terra::hist(ndvi_mosaic, breaks = 70, yaxt = 'n')
```

```{r echo = FALSE, fig.align='center', fig.height=4, fig.width = 5, warning = FALSE, message = FALSE, dpi = 300, out.width="60%", fig.cap="**Figure: Distribution of NDVI values across the Punggol area in Singapore, based on a mosaic of images collected from 2019-07-01 to 2021-06-30.** Pixel values to the right of the dashed line may be classified as 'vegetated', while those to the left are 'non-vegetated'."}
plot(histogram, 
     xlab = "NDVI", xlim = c(-1, 1), main = "", yaxt = 'n',
      col = rgb(178/255,223/255,138/255, 0.5))
abline(v = threshold, col="black", lwd=2, lty = 2)
```

<br>

The function `classify_image_binary()` can be used to directly classify a continuous raster (e.g., image mosaic). Note that it uses `threshold_otsu()` internally to define the threshold value. See the map below to examine both the original NDVI and classified rasters within the Punggol area in Singapore (toggle the map layers to view each one separately).

```{r mosaic classification, warning = FALSE, message = FALSE, fig.width=2.6, fig.height = 2.0, dpi = 300, out.width="100%"}
veg_classified <- classify_image_binary(ndvi_mosaic, threshold = "otsu")

tm_basemap(c("CartoDB.Positron")) +
  tm_shape(ndvi_mosaic) +
    tm_raster(palette = c("Greens")) +
  tm_shape(veg_classified) +
    tm_raster(style = "cat",
            palette = c("grey", "darkgreen")) +
    tm_shape(points) + tm_dots()
```

<br>

'Landscape metrics' can then be used to represent the composition and spatial patterns of land cover within the area of interest. For each classified raster, metrics can be summarised at specific point locations and at multiple buffer radii. Metrics can also be calculated at various 'levels', or unit of analysis — _patch_, _class_, and _landscape_ (in order of widening spatial scale). Refer to the package [landscapemetrics](https://r-spatialecology.github.io/landscapemetrics/index.html) for more information.

```{r info on landscape metrics, echo = FALSE}

lsm <- landscapemetrics::list_lsm() %>%
  dplyr::select(-function_name) %>%
  filter(level != "patch") %>%
  pivot_wider(names_from = "level", values_from = "name") %>%
  filter(metric != "ta", # remove unnecessary metrics
         #metric != "pr",
         #metric != "prd",
         metric != "rpr") %>% 
  mutate(name = coalesce(class, landscape)) %>% # new col for name
  relocate(name, .after = metric) %>%
  mutate(across(.cols = c("class", "landscape"), # simplify to Yes/NA 
                ~ifelse(is.na(.), NA, "Yes"))) %>%
  mutate(name = case_when( # add details to name
    grepl("_cv", metric) ~ paste(name, "(coef. of variation)"), 
    grepl("_mn", metric) ~ paste(name, "(mean)"),
    grepl("_sd", metric) ~ paste(name, "(st. dev.)"),
    TRUE ~ name
  )) %>%
  arrange(class, landscape, metric, type) %>%
  mutate(class = ifelse(is.na(class), "No", paste0("[", class, "](https://r-spatialecology.github.io/landscapemetrics/reference/lsm_c_", metric, ".html)"))) %>% 
  mutate(landscape = ifelse(is.na(landscape), "No", paste0("[", landscape, "](https://r-spatialecology.github.io/landscapemetrics/reference/lsm_l_", metric, ".html)"))) 


knitr::kable(lsm, format = "html", caption = "**Table: Examples of landscape metrics and their corresponding level(s) they may be calculated at (e.g., class, landscape)**. Details in hyperlinked text.") %>%
  kableExtra::kable_styling(full_width = T, position = "left")  %>%
    kableExtra::kable_styling("striped") %>% kableExtra::scroll_box(width = "100%", height = "600px")

rm(lsm)
```

<br>

A wide variety of landscape metrics can be summarised at point locations of interest, using the function `calc_lsm()`. The following example summarises all 'class-level' metrics for vegetation cover (`veg_classified`), up to `200` and `400` metres away from `points`. This returns a list, with each element corresponding to a particular buffer radius. Column names include a prefix indicating that values are landscape metrics derived from a classified raster (`lsm_`).

```{r calc_lsm, warning = FALSE, message = FALSE}
lsm_perpoint <- 
  calc_lsm(raster = veg_classified,
           points = points,
           buffer_sizes = c(200, 400), # in metres
           level = c("class"),
           class_names = c("veg"), # land cover class name
           class_values = c(1)) # pixel value of class

lsm_perpoint
```

This list object can then be converted into a dataframe, with a new column indicating the radii in which these landscape metrics have been summarised.

```{r warning = FALSE, message = FALSE}
lsm_perpoint %>% 
  bind_rows(.id = "radius")
```


<br>

---

## 2. OpenStreetMap landscape elements

Landscape elements can be retrieved from the open-source maps such as OpenStreetMap (OSM). For instance, built elements such as building polygons and road lines can be extracted from `sampling_areas` using the functions `get_buildings_osm()` and `get_roads_osm()`, respectively. The `date` argument, if specified, must be a snapshot date present within the [Geofabrik database](https://download.geofabrik.de). If not, the latest data will be downloaded. In this example, we download these urban elements for the date 1 Jan 2021 and view them as an interactive map:

```{r download urban elements, eval = FALSE}

dir.create(paste0(output_dir, "/osm")) # create sub-directory for raw data (avoid re-downloading each time)

buildings_osm <- get_buildings_osm(place = sampling_areas,
                               date = as.Date("2021-01-01"),
                               dir_raw = paste0(output_dir, "/osm"))

roads_osm <- get_roads_osm(sampling_areas,
                       date = as.Date("2021-01-01"),
                       dir_raw = paste0(output_dir, "/osm"))

# alternatively, load example dataset from package
filepath <- system.file("extdata", "osm_data.Rdata", package = "biodivercity")
load(filepath)

# plot
tm_basemap(c("CartoDB.Positron")) +
  tm_shape(sampling_areas) + tm_borders() +
  tm_shape(buildings_osm) + tm_polygons(col = "levels") +
  tm_shape(roads_osm) + tm_lines(col = "lanes", palette = "YlOrRd")+
  tm_shape(points) + tm_dots()
```

```{r load and view osm elements, echo = FALSE, warning = FALSE, message = FALSE, fig.width=2.6, fig.height = 2.0, dpi = 300, out.width="100%"}

filepath <- system.file("extdata", "osm_data.Rdata", package = "biodivercity")
load(filepath)

tmap_mode("view")
tm_basemap(c("CartoDB.Positron")) +
  tm_shape(sampling_areas) + tm_borders() +
  tm_shape(buildings_osm) + tm_polygons(col = "levels") +
  tm_shape(roads_osm) + tm_lines(col = "lanes", palette = "YlOrRd")+
  tm_shape(points) + tm_dots()
```

<br>

The function `calc_osm()` can be be used to summarise the following metrics of buildings and roads, at multiple buffer radii. 

- Buildings (polygons):
    - Building footprint area (`buildingArea_m2`)
    - Building gross floor area (`buildingGFA_m2`) - _Column for the no. of levels must be present_
    - Average number of levels (`buildingAvgLvl`) - _Column for the no. of levels must be present_
    - Building Volume (`buildingVol_m3`) - _Column for height must be present_
    - Floor area ratio (`buildingFA_ratio`) - _Column for the no. of levels must be present_    

- Roads (lines):
    - Total lane length (`lanelength`) - _Column for the no. of lanes must be present_
    - Lane density (`laneDensity`) - _Column for the no. of lanes must be present_

<br>

For example, we can calculate building metrics up to `200` and `400` metres away from `points`. This returns a list, with each element corresponding to a particular buffer radius. Column names include a prefix indicating that values are from OpenStreetMap (`osm_`).

```{r calc_osm, warning = FALSE, message = FALSE}

osm_perpoint <- 
  calc_osm(vector = buildings_osm, 
           points = points,
           name = "buildings",
           buffer_sizes = c(200, 400),
           building_levels = "levels") 

osm_perpoint
```

This list object can then be converted into a dataframe, with a new column indicating the radii in which these building metrics have been summarised.

```{r warning = FALSE, message = FALSE}
do.call(rbind, # alternative way to convert list to dataframe
       mapply(transform, osm_perpoint, 
              radius = names(osm_perpoint), # new column
              SIMPLIFY = FALSE))
```

<br>

---

## 3. Manually generated landscape elements

There may be cases where you intend to summarise data that have been generated manually, for instance, from on-site mapping or artificial design scenarios (see `vignette("use-cases")`). Load example data for landscape vectors (points, polygons, lines) within the Punggol (PG) area in Singapore (Chong et al., 2014, 2019), visualised in the interactive map below. Note that there are no water bodies mapped within the target area.

```{r import manually mapped vector data, message=FALSE}
filepath <- system.file("extdata", "landscape-vectors_mapped.Rdata", package = "biodivercity") # load example data
load(filepath) 
```

```{r plot interative map for landscape data, echo = FALSE, warning = FALSE, message = FALSE, fig.width=2.6, fig.height = 2.0, dpi = 300, out.width="100%"}

tmap_mode("view")
tmap_options(check.and.fix = TRUE)

tm_basemap(c("CartoDB.Positron")) +
  tm_shape(bound) +
    tm_borders() +
  tm_shape(natveg) +
    tm_polygons(title = "Natural vegetation",
                group = "natveg",
                col = "darkgreen",
                alpha = 0.6,
                border.col = "transparent") +
  tm_shape(trees %>% relocate(species)) +
    tm_dots(title = "Trees",
            group = "trees",
            col = "brown",
            size = 0.001,
            border.col = "transparent") +
  tm_shape(shrubs %>% relocate(species) %>% st_make_valid()) +
    tm_polygons(title = "Shrubs",
                group = "shrubs",
                col = "#ffd92f",
                alpha = 0.6,
                border.col = "transparent") +
  tm_shape(turf) +
    tm_polygons(title = "Turf",
                group = "turf",
                col = "#a6d854",
                alpha = 0.6,
                border.col = "transparent") +
  # tm_shape(water) +
  #   tm_polygons(title = "Water",
  #               group = "water",
  #               col = "blue",
  #               alpha = 0.6,
  #               border.col = "transparent") +
  tm_shape(buildings) +
    tm_polygons(title = "Building levels",
                group = "buildings",
                col = "levels",
                palette = viridis::magma(5),
                style = "jenks",
                border.col = "transparent") + 
  tm_shape(roads) +
    tm_lines(title.col = "Road lanes",
             group = "roads",
             col = "lanes",
             palette = viridis::magma(5),
             style = "fixed",
             breaks = c(2, 4, 6, 8))

```

The function `calc_manual()` is provided for interested users, to calculate the following landscape predictors at multiple buffer radii. Column names include a prefix indicating that values were manually generated (`man_`).
    
- Buildings (polygons):
    - Building footprint area (`buildingArea_m2`)
    - Building gross floor area (`buildingGFA_m2`) - _Column for the no. of levels must be present_
    - Average number of levels (`buildingAvgLvl`) - _Column for the no. of levels must be present_
    - Floor area ratio (`buildingFA_ratio`) - _Column for the no. of levels must be present_    

- Roads (lines):
    - Total lane length (`lanelength`) - _Column for the no. of lanes must be present_
    - Lane density (`laneDensity`) - _Column for the no. of lanes must be present_
    
- Trees (points):
    - Count of trees (`tree_count`)
    - Species richness of trees (`tree_sprich`) - _Column with species name must be present_

- Shrubs (polygons):
    - Percentage of landscape area with shrubs (`shrub_pland`)
    - Species richness of shrubs (`shrub_sprich`) - _Column with species name must be present_

- Turf (polygons):
    - Percentage of landscape area with turf (`turf_pland`)

- Natural vegetation (polygons):
    - Percentage of landscape area with natural vegetation (`natveg_pland`)

- Water (polygons):
    - Percentage of landscape area with water (`water_pland`)

<br>

For example, we can summarise the species and count of trees up to `50` metres away from `points`. This returns a list, with each element corresponding to a particular buffer radius.

```{r calculate metrics for manually generated data, message = FALSE}
manual_perpoint <-
  calc_manual(vector = trees, name = "trees",
              points = points, buffer_sizes = 50,
              plant_species = "species")

manual_perpoint
```

This list object can then be converted into a dataframe, with a new column indicating the radii in which these metrics have been summarised.

```{r warning = FALSE, message = FALSE}
manual_perpoint %>% 
  bind_rows(.id = "radius")
```


---

## References

Chong KY, Teo S, Kurukulasuriya B, Chung YF, Giam X & Tan HTW (2019) The effects of landscape scale on greenery and traffic relationships with urban birds and butterflies. _Urban Ecosystems_, _22_(5): 917–926.

Chong KY, Teo S, Kurukulasuriya B, Chung YF, Rajathurai S & Tan HTW (2014) Not all green is as good: Different effects of the natural and cultivated components of urban vegetation on bird and butterfly diversity. _Biological Conservation_, _171_: 299–309.