In SEOSAW we have a few very large woodland monitoring plots. Six 25 ha plots in Republic of Congo, a 10 ha plot in the Democratic Republic of Congo, and a 4 ha plot in South Africa. SEOSAW advocates for using 1 ha square plots (100x100 m) in most cases, as do many other similar ecological plot networks, because they are big enough to be linked to remote sensing data, and big enough to include what is normally a representative sample of the woodland.
To increase the usefulness of the very large plots, we wanted to split them up into 1 ha square subdivisions, and attribute stems to each of these subdivisions based on stem location.
Here is the R function I wrote to split the plots:
#' Split large plots into 1 hectare squares
#'
#' @param x dataframe of stem measurements
#' @param plot_data dataframe of plot measurements
#' @param dims named character vector of length 2, specifying the relationship
#' between the xy stem grid coordinate system and the plot dimensions, i.e.
#' does the x axis run along the plot length, or the plot width measurement
#' @param xy_zero_corner either "sw", "nw", "ne", "se", specifying the corner
#' of the grid coordinate system origin (0,0)
#' @param x_direction either "ew" (east-west) or "ns" (north-south), specifying
#' the direction of the x axis on the grid coordinate system from the
#' origin corner
#' @param subplot_size the length and width of square subplots, e.g. for 1 ha,
#' 100
#'
#' @return list:
#' 1) updated dataframe of stem measurements with added \code{subdiv_id} to
#' show which subdivision each stem belongs to.
#' 2) sf dataframe of approximated subplot polygons and their centres.
#' 3) updated dataframe of plot measurements with added \code{subdiv} to
#' show whether a plot has been subdivided, and \code{subdiv_id} to
#' identify new rows which are subdivided plots
#'
#' @details This function will only split a plot if it is rectangular and the
#' area can be divided into 1 hectare squares without remainder.
#'
#' @examples
#'
#'
#' @export
#'
largePlotSplit <- function(x, plot_data, polys, stem_plot_id = "plot_id",
plot_plot_id = stem_plot_id, polys_plot_id = stem_plot_id,
dims = c("plot_length" = "x_grid", "plot_width" = "y_grid"),
xy_zero_corner = "sw", x_direction = "ew", subplot_size = 100) {
# Check parameter specified correctly
stopifnot(!any(is.na(names(dims))) & !is.null(names(dims)) & length(dims)==2)
stopifnot(xy_zero_corner %in% c("sw", "se", "nw", "ne"))
stopifnot(x_direction %in% c("ns", "ew"))
# Subset plots which are the right dimensions to split and which contain xy grid coords
plots_fil <- plot_data[
!is.na(plot_data[[names(dims)[1]]]) &
plot_data[[names(dims)[1]]] %% 100 == 0 &
plot_data[[names(dims)[1]]] %/% 100 >=2 &
!is.na(plot_data[[names(dims)[2]]]) &
plot_data[[names(dims)[2]]] %% 100 == 0 &
plot_data[[names(dims)[2]]] %/% 100 >=2 &
plot_data[[plot_plot_id]] %in%
unique(x[!is.na(x[[dims[1]]]) & !is.na(x[[dims[2]]]), stem_plot_id])
,]
# Subset stems by plot IDs
x_fil <- x[x[[stem_plot_id]] %in% plots_fil[[plot_plot_id]],]
# Split stem measurements by plot
x_split <- split(x_fil, x_fil[[stem_plot_id]])
# For each plot:
plot_out <- lapply(seq_along(x_split), function(y) {
# New object, easier to work with
x_split_iso <- as.data.frame(x_split[[y]])
# Get plot ID
plot_id <- unique(x_split_iso[[stem_plot_id]])
# Throw warning if any stems not matched
no_coords <- which(is.na(x_split_iso[[dims[1]]]) | is.na(x_split_iso[[dims[2]]]))
if (length(no_coords) > 0) {
warning(plot_id, ": discarded ", length(no_coords),
" stems with no XY grid coordinates")
}
x_split_fil <- x_split_iso[-no_coords,]
# Define bins per plot
cut_length <- plots_fil[plots_fil[[plot_plot_id]] == plot_id,names(dims)[1]] %/% subplot_size
cut_width <- plots_fil[plots_fil[[plot_plot_id]] == plot_id,names(dims)[2]] %/% subplot_size
bin_length <- seq(0, cut_length * subplot_size, by = subplot_size)
bin_width <- seq(0, cut_width * subplot_size, by = subplot_size)
# Classify each grid point by bins
x_split_fil$length_bin <- cut(x_split_fil[[dims[1]]], bin_length)
x_split_fil$width_bin <- cut(x_split_fil[[dims[2]]], bin_width)
# Get factor labels
length_bin_labels <- levels(x_split_fil$length_bin)
width_bin_labels <- levels(x_split_fil$length_bin)
# Deal with values slightly outside bins
# Length
poss_fix_length <- x_split_fil[is.na(x_split_fil$length_bin) &
!is.na(x_split_fil[[dims[1]]]), dims[1]]
if (length(poss_fix_length) > 0) {
bin_length_close <- closestMatch(bin_length, unlist(poss_fix_length))
bin_length_labs <- sapply(bin_length_close, function(i) {
if (i == min(bin_length)) {
length_bin_labels[1]
} else {
length_bin_labels[which(i == bin_length) -1]
}
})
x_split_fil[
is.na(x_split_fil$length_bin) &
!is.na(x_split_fil[[dims[1]]]), "length_bin"] <- bin_length_labs
}
# Width
poss_fix_width <- x_split_fil[is.na(x_split_fil$width_bin) &
!is.na(x_split_fil[[dims[2]]]), dims[2]]
if (length(poss_fix_width) > 0) {
bin_width_close <- closestMatch(bin_width, unlist(poss_fix_width))
bin_width_labs <- sapply(bin_width_close, function(i) {
if (i == min(bin_width)) {
width_bin_labels[1]
} else {
width_bin_labels[which(i == bin_width)-1]
}
})
x_split_fil[
is.na(x_split_fil$width_bin) &
!is.na(x_split_fil[[dims[2]]]), "width_bin"] <- bin_width_labs
}
# For each unique combination of bins, make a subplot, with name
y_split <- split(x_split_fil, list(x_split_fil$length_bin, x_split_fil$width_bin))
poss_subsets <- paste0(plot_id, "_S", seq(length(y_split)))
x_split_fil <- do.call(rbind, lapply(seq_along(y_split), function(z) {
if (nrow(y_split[[z]]) > 0) {
y_split[[z]]$subdiv_id <- poss_subsets[z]
y_split[[z]]
}
}))
# Filter polygons to current plot ID
polys_fil <- polys[polys[[polys_plot_id]] == plot_id,]
# Extract corners as dataframe
polys_points <- as.data.frame(sf::st_coordinates(polys_fil)[1:4,1:2])
# Get UTM zone of corners
utm_string <- UTMProj4(latLong2UTM(mean(polys_points[,1]), mean(polys_points[,2])))
# Convert polygons to UTM
polys_utm <- sf::st_transform(polys_fil, utm_string)
# Convert UTM polygons to points
points_utm <- sf::st_cast(polys_utm, "POINT", warn = FALSE)[1:4,]
# Extract coordinates as dataframe
coords_utm <- as.data.frame(sf::st_coordinates(points_utm)[1:4,1:2])
# Get zero and bearing point
nw_outer <- sf::st_sfc(sf::st_point(
x = c(mean(coords_utm$X) - 1000, mean(coords_utm$Y) + 1000)))
ne_outer <- sf::st_sfc(sf::st_point(
x = c(mean(coords_utm$X) + 1000, mean(coords_utm$Y) + 1000)))
sw_outer <- sf::st_sfc(sf::st_point(
x = c(mean(coords_utm$X) - 1000, mean(coords_utm$Y) - 1000)))
se_outer <- sf::st_sfc(sf::st_point(
x = c(mean(coords_utm$X) + 1000, mean(coords_utm$Y) - 1000)))
if (xy_zero_corner == "sw") {
match_point <- sw_outer
other_point <- nw_outer
} else if (xy_zero_corner == "nw") {
match_point <- nw_outer
other_point <- sw_outer
} else if (xy_zero_corner == "ne") {
match_point <- ne_outer
other_point <- se_outer
} else if (xy_zero_corner == "se") {
match_point <- se_outer
other_point <- ne_outer
}
# Set CRS
sf::st_crs(other_point) <- sf::st_crs(points_utm)
sf::st_crs(match_point) <- sf::st_crs(points_utm)
# Get xy zero and opposite corner
xy_zero <- points_utm[sf::st_nearest_feature(match_point, points_utm),]
xy_1 <- points_utm[sf::st_nearest_feature(other_point, points_utm),]
# Convert to WGS for geosphere compatibility
xy_zero_wgs <- sf::st_coordinates(sf::st_transform(xy_zero, 4326))
xy_1_wgs <- sf::st_coordinates(sf::st_transform(xy_1, 4326))
# Get bearing between points
xy_bearing <- geosphere::bearing(xy_zero_wgs, xy_1_wgs)
# Make grid of 1 ha plots
bbox_grid <- cbind(x = c(0, 0, max(bin_length), max(bin_length), 0),
y = c(0, max(bin_width), max(bin_width), 0, 0))
bbox_sf <- sf::st_polygon(list(bbox_grid))
polys_grid <- sf::st_make_grid(bbox_sf, cellsize = subplot_size)
sf::st_crs(polys_grid) <- sf::st_crs(points_utm)
# Rotate to angle and move to centre of original plot
cent <- sf::st_centroid(sf::st_combine(points_utm))
grid_cent <- sf::st_centroid(sf::st_combine(polys_grid))
angle <- NISTunits::NISTdegTOradian(xy_bearing)
polys_rot <- (polys_grid - grid_cent) * rot(angle) * 1 + cent
sf::st_crs(polys_rot) <- sf::st_crs(cent)
# Set order of subset IDs
byrow <- ifelse(x_direction == "ns", FALSE, TRUE)
if (xy_zero_corner %in% c("se", "nw")) {
ord <- cut_width:1
} else {
ord <- 1:cut_width
}
if (xy_zero_corner %in% c("sw", "nw")) {
subdiv_ids <- poss_subsets
} else {
subdiv_ids <- rev(poss_subsets)
}
# Create matrix of subset IDs
subset_mat <- matrix(subdiv_ids,
nrow = cut_width, ncol = cut_length, byrow = byrow)[ord,]
# Apply matrix to name polygons
polys_rot_sf <- sf::st_sf(geometry = polys_rot,
plot_id = plot_id,
subdiv_id = as.vector(t(subset_mat)))
polys_rot_wgs <- sf::st_transform(polys_rot_sf, 4326)
names(polys_rot_wgs)[1] <- plot_plot_id
# Remove bins
x_split_fil_clean <- x_split_fil[,!(names(x_split_fil) %in%
c("length_bin", "width_bin"))]
# Add back in stems which had no grid coords
x_split_iso$subdiv_id <- NA_character_
x_split_out <- rbind(x_split_fil_clean,
x_split_iso[no_coords,names(x_split_fil_clean)])
# Return list
return(list(x_split_out, polys_rot_wgs))
})
# Bind all stems together
stems_out <- do.call(rbind, lapply(plot_out, "[[", 1))
# Bind back in stems from plots which weren't split
if (nrow(x[!x[[stem_plot_id]] %in% plots_fil[[plot_plot_id]],]) > 0) {
stems_non <- x[!x[[stem_plot_id]] %in% plots_fil[[plot_plot_id]],]
stems_non$subdiv_id <- NA_character_
stems_all <- rbind(stems_out, stems_non[,names(stems_out)])
} else {
stems_all <- stems_out
}
# Bind all polygons
polys_out <- do.call(rbind, lapply(plot_out, "[[", 2))
# Add centres to polygons
centres <- suppressWarnings(
as.data.frame(st_coordinates(st_centroid(polys_out))))
names(centres) <- c("longitude_of_centre", "latitude_of_centre")
polys_out <- cbind(polys_out, centres)
# Add subdiv column to plots table
plot_data$subdiv <- ifelse(plot_data[[plot_plot_id]] %in%
polys_out[[plot_plot_id]], TRUE, FALSE)
# Add subdivided columns to plots data
subdiv_out <- merge(st_drop_geometry(polys_out),
plot_data[,!names(plot_data) %in% c("longitude_of_centre", "latitude_of_centre")],
by = plot_plot_id, all.x = TRUE)
subdiv_out[subdiv_out$subdiv_id %in% polys_out$subdiv_id, "plot_area"] <- subplot_size * subplot_size / 10000
subdiv_out[subdiv_out$subdiv_id %in% polys_out$subdiv_id, "plot_length"] <- subplot_size
subdiv_out[subdiv_out$subdiv_id %in% polys_out$subdiv_id, "plot_width"] <- subplot_size
subdiv_out[subdiv_out$subdiv_id %in% polys_out$subdiv_id, "plot_perimeter"] <- subplot_size * 4
subdiv_out[subdiv_out$subdiv_id %in% polys_out$subdiv_id, "subdiv"] <- FALSE
plot_data$subdiv_id <- NA_character_
plot_data_out <- rbind(plot_data, subdiv_out[,names(plot_data)])
# Return list of output
return(list(stems_all, polys_out, plot_data_out))
}
Hopefully the comments explain everything, but this function is quite complicated, and I’m proud of it, so I’ll describe it step by step.
- Find which plots are suitable for splitting. Plots must be rectangular, larger than 2 ha, the plot length and width must be divisible without remainder into 100 m lengths, and there must be at least one stem in the plot with grid coordinates. Then, for each suitable plot:
- Assign each stem to 100 m bins in the X and Y direction, using
cut(). - For stems with grid coordinates slightly outside the plot boundary (this is often a problem when converting from lat-long to grid coordinates), find the closest bin, using a separate function called
closestMatch(), see below for more details. - For each unique combination of bins, make a subdivision ID, and assign names based on the direction of stem coordinate progression and which corner the stem coordinate system originates from, referencing the
x_directionandxy_zero_cornerfunction parameters. - Make a grid of 1 ha square polygons from the plot dimensions, using
sf::st_make_grid(). Then rotate the polygon grid and move to the centre of the existing plot polygon. - Name each of the subdivision polygons with the
subdiv_idused for the stems. - Return a dataframe of stem measurements with
subdiv_idvalues for each stem, and an{sf}object with the subdivision polygons and their centres.
The plot below illustrates a 4 ha plot that has been split into four 1 ha subdivisions using the above methods.
library(ggplot2)
library(ggnewscale)
library(viridis)
ggplot() +
geom_sf(data = polys_sf, colour = "black", size = 2, fill = NA) +
geom_sf(data = out[[2]], aes(colour = subdiv_id), size = 1.5, fill = NA) +
new_scale_colour() +
geom_point(data = out[[1]],
aes(x = longitude, y = latitude, colour = x_grid, shape = subdiv_id)) +
scale_colour_viridis() +
theme_bw() +
labs(x = "", y = "")
Finally, here is the closestMatch() function:
#' Find closest match in a vector
#'
#' @param x numeric vector
#' @param y vector of numeric values, each of which be matched with the closest
#' value in \code{x}
#'
#' @return numeric vector of values of \code{x} that are closest to each
#' element of \code{y}
#'
#' @keywords internal
#' @noRd
#'
closestMatch <- function(x, y) {
unlist(lapply(y, function(i) {
x[which(abs(x-i) == min(abs(x-i)))]
}))
}