Spatial Statistics with ForestTools"
In ForestTools: Tools for Analyzing Remote Sensing Forest Data

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.

# 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.

# 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.

# 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.

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