Skip to content

Latest commit

 

History

History
107 lines (88 loc) · 3.66 KB

spatial_statistics.md

File metadata and controls

107 lines (88 loc) · 3.66 KB
Raw

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.

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

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