Spatial Statistics with ForestTools
================
Andrew Plowright
2024-01-19

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
# 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)
```

![](spatial_statistics_files/figure-gfm/unnamed-chunk-1-1.png)<!-- -->

# Global Statistics

Calculating global statistics for the trees is fairly straightforward.

``` r
# Height
mean(kootenayTrees$height)
## [1] 5.251959
median(kootenayTrees$height)
## [1] 4.256785
max(kootenayTrees$height)
## [1] 13.49121

# Crown area
mean(kootenayCrowns$crownArea)
## [1] 9.006453
median(kootenayCrowns$crownArea)
## [1] 7
max(kootenayCrowns$crownArea)
## [1] 54.25
```

# 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
# 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")]
```

    ## Simple feature collection with 3 features and 5 fields
    ## Geometry type: MULTIPOLYGON
    ## Dimension:     XY
    ## Bounding box:  xmin: 439689.1 ymin: 5526454 xmax: 439832.5 ymax: 5526563
    ## Projected CRS: WGS 84 / UTM zone 11N
    ##   BlockID mean_height max_height mean_crownArea max_crownArea
    ## 0     101    7.338588  13.491207      13.683476         51.25
    ## 1    3308    3.208017   7.125149       4.847403         23.50
    ## 2     113    7.685308  12.583441      13.430990         54.25
    ##                         geometry
    ## 0 MULTIPOLYGON (((439707.9 55...
    ## 1 MULTIPOLYGON (((439832.5 55...
    ## 2 MULTIPOLYGON (((439832.4 55...
