💾 Archived View for republic.circumlunar.space › users › johngodlee › posts › 2021-05-10-winkelmass.… captured on 2023-04-19 at 23:49:56. Gemini links have been rewritten to link to archived content
⬅️ Previous capture (2021-12-04)
-=-=-=-=-=-=-
DATE: 2021-05-10
AUTHOR: John L. Godlee
Winkelmass means "angle measure" in German. In forestry the term has been adopted to refer to a specific metric used to describe the regularity of tree spatial distribution. Literature on the winkelmass traces mostly to Klaus von Gadow, a researcher split between the University of Göttingen and Beijing Forestry University. The winkelmass is quite simple at it's core. For each tree, a given number of nearest neighbours are identified, usually four. Then, going clockwise around the focal tree, the angle between each neighbour and the next neighbour is calculated. For each angle, if the angle is less than the 'critical angle' i.e. 360 degrees divided by the number of neighbours considered, a value of 1 is recorded, otherwise, a value of 0 is recorded. Finally, the sum of these values divided by the number of neighbours gives the winkelmass. In LaTeX math:
W_{i} = \frac{1}{k} \sum_{j=1}^{k} v_{j}
where W is the winkelmass, k is the number of neighbours considered, and v_{j} is either 0 or 1 depending on the angle between two neighbours.
To calculate the plot-level winkelmass, simply take the mean of all W_{i} values in the plot.
I wrote a function in R to calculate the Winkelmass for all trees in a plot:
#' Winkelmass (spatial regularity of trees)
#'
#' @param x vector of individual x axis coordinates
#' @param y vector of individual y axis coordinates
#' @param k number of neighbours to consider
#'
#' @references von Gadow, K., Hui, G. Y. (2001). Characterising forest spatial
#' structure and diversity. Sustainable Forestry in Temperate Regions. Proc. of
#' an international workshop organized at the University of Lund, Sweden.
#' Pages 20- 30.
#'
winkelmass <- function(x, y, k = 4) {
dat_sf <- sf::st_as_sf(data.frame(x,y), coords = c("x", "y"))
dists <- suppressMessages(nngeo::st_nn(dat_sf, dat_sf, k = k+1,
progress = FALSE))
a0 <- 360 / k
wi <- unlist(lapply(dists, function(i) {
focal_sfg <- sf::st_geometry(dat_sf[i[1],])[[1]]
nb_sfg <- sf::st_geometry(dat_sf[i[-1],])
nb_angles <- sort(unlist(lapply(nb_sfg, function(j) {
angleCalc(focal_sfg, j)
})))
aj <- nb_angles - c(NA, head(nb_angles, -1))
aj[1] <- nb_angles[k] - nb_angles[1]
aj <- ifelse(aj > 180, 360 - aj, aj)
aj <- round(aj, 1)
sum(aj < a0)
}))
out <- 1 / k * wi
return(out)
}
The function uses {nngeo} to find nearest neighbours of each tree by their x and y grid coordinates. Then it makes {sf} geometry objects and calculates the angles between them, and finally does a bit of arithmetic to find if the angles are less than the critical angle. An important aspect to the calculation is that if the angle between two neighbours is greater than 180 degrees, the opposite angle is taken. So if an angle is 275 degrees, the recorded angle is 360 - 275 = 85 degrees. In this example by taking the opposite angle the value of v_{j} changes from 0 (275 gt 90) to 1 (85 lt 90).
The angleCalc() function is:
#' Calculate angle between two sf point objects
#'
#' @param x point feature of class 'sf'
#' @param y point feature of class 'sf'
#'
#' @return azimuthal from x to y, in degrees
#'
#' @examples
#' p1 <- st_point(c(0,1))
#' p2 <- st_point(c(1,2))
#' angleCalc(p1, p2)
#'
#' @export
#'
angleCalc <- function(x, y) {
dst_diff <- as.numeric(x - y)
return((atan2(dst_diff[1], dst_diff[2]) + pi) / 0.01745329)
}
I also wrote some tests for the function to make sure it worked properly. First, packages, all for visualising results, not needed for calculations:
library(ggplot2) library(patchwork) library(knitr)
Define a function to create maps to visualise the distribution of neighbours around a focal tree:
linemake <- function(dat) {
out <- data.frame(
xmin = rep(dat$x[1], times = nrow(dat) - 1),
ymin = rep(dat$y[1], times = nrow(dat) - 1),
xmax = dat$x[-1], ymax = dat$y[-1])
return(out)
}
plotmake <- function(dat) {
lout <- linemake(dat)
dat$label <- c("F", "1", "2", "3", "4")
ggplot() +
geom_segment(data = lout,
aes(x = xmin, y = ymin, xend = xmax, yend = ymax)) +
geom_label(data = dat, aes(x = x, y = y, label = label)) +
coord_equal() +
theme_classic() +
labs(x = "X", y = "Y")
}
Create a number of test cases with different spatial distributions of neighbours. Each test case should produce a different value of Wi for the focal tree.
test_list <- list(
test1 = data.frame(x = c(0,0,1,0,-1), y = c(0,1,0,-1,0)),
test2 = data.frame(x = c(0,0,1,0,-0.5), y = c(0,1,0,-1,0.5)),
test3 = data.frame(x = c(0,0,0.5,0,-0.5), y = c(0,1,0.5,-1,0.5)),
test4 = data.frame(x = c(0,0,0.2,-0.4,-0.5), y = c(0,1,0.5,-0.1,0.5)),
test5 = data.frame(x = c(0,0,0.2,-0.5,-0.2), y = c(0,1,0.9,0.4,0.5))
)
Create the maps to visualise layout of the test cases:
plot_list <- lapply(test_list, plotmake)
wrap_plots(plot_list) +
plot_layout(nrow = 1)
Finally run the winkelmass() function and extract the first value which refers to the focal tree:
kable(
unlist(lapply(test_list, function(x) {
winkelmass(x$x, x$y)[1]
})))
┌───────┬───────┐ │ Ex │ Wi │ ╞═══════╪═══════╡ │ test1 │ 0 │ ├───────┼───────┤ │ test2 │ 0.250 │ ├───────┼───────┤ │ test3 │ 0.500 │ ├───────┼───────┤ │ test4 │ 0.750 │ ├───────┼───────┤ │ test5 │ 1 │ └───────┴───────┘
As you can see, for a given value of k, there are k+1 possible values of Wi depending on the distribution of neighbours.