---
title: "Spatial Statistics with ForestTools"
author: "Andrew Plowright"
date: "2024-01-19"
output: github_document
---

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

This vignette will take you through the process of generating tree statistics over polygonal areas using the outputs from `ForestTools`.

Start by loading in some test data. These are sample outputs from the `ForestTools` functions. `kootenayTrees` are individual tree points produced using `vwf`. `kootenayCrowns` are polygonal crowns generated by `mcws`. `kootenayBlocks` are some polygonal areas of interest that we'll use for generating statistics.

```{r message=FALSE, warning=FALSE}
# Load libraries
library(terra)
library(sf)
library(ForestTools)
library(dplyr)
library(magrittr)

# Load sample data
data("kootenayTrees", "kootenayBlocks", "kootenayCrowns", "kootenayCHM")

# Plot areas of interest and add trees
plot(unwrap(kootenayCHM))
plot(kootenayBlocks$geometry, add=T)
plot(kootenayTrees$geometry, pch=19, col="blue", add=T)
```

# Global Statistics

Calculating global statistics for the trees is fairly straightforward.

```{r, collapse=TRUE}
# Height
mean(kootenayTrees$height)
median(kootenayTrees$height)
max(kootenayTrees$height)

# Crown area
mean(kootenayCrowns$crownArea)
median(kootenayCrowns$crownArea)
max(kootenayCrowns$crownArea)
```

# Statistics by Polygon

To calculate statistics by polygonal areas of interest, we'll first use the `st_intersects` function to create subsets of trees for each polygon. We'll then calculate statistics for each of those subsets, and return the result to polygons. Note: this can duplicate trees that are contained in overlapping polygons.

```{r}
# Create subsets of the trees for each polygonal area of interest 
trees_by_poly <- kootenayBlocks %>% 
  st_intersects(kootenayTrees) %>%
  lapply(function(x) kootenayTrees[x,])

# Calculate statistics for each polygonal area of interest
kootenayBlocks[["mean_height"]] <- sapply(trees_by_poly, function(trees) mean(trees$height))
kootenayBlocks[["max_height"]]  <- sapply(trees_by_poly, function(trees) max(trees$height))

```

The same operation can be repeated for tree crowns. However, given that the crowns are polygonal, we don't want a single crown to be counted twice between two adjoining areas of interest. The only variation required here is that we use `st_centroid` to compute centroids for each crown before intersecting them with the areas of interest.

```{r, warning=FALSE}
# Create subsets of the crowns for each polygonal area of interest 
crowns_by_poly <- kootenayBlocks %>% 
  st_intersects(st_centroid(kootenayCrowns)) %>%
  lapply(function(x) kootenayCrowns[x,])

# Calculate statistics for each polygonal area of interest
kootenayBlocks[["mean_crownArea"]] <- sapply(crowns_by_poly, function(crowns) mean(crowns$crownArea))
kootenayBlocks[["max_crownArea"]]  <- sapply(crowns_by_poly, function(crowns) max(crowns$crownArea))

# View results
kootenayBlocks[,c("BlockID", "mean_height", "max_height", "mean_crownArea", "max_crownArea")]
```

