TITLE: Splitting a large woodland survey plot into 1 ha subdivisions
DATE: 2021-03-20
AUTHOR: John L. Godlee
====================================================================


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.

 [SEOSAW]: https://seosaw.github.io

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.

1.  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:
2.  Assign each stem to 100 m bins in the X and Y direction, using
cut().
3.  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.
4.  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_direction and xy_zero_corner function parameters.
5.  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.
6.  Name each of the subdivision polygons with the subdiv_id used
for the stems.
7.  Return a dataframe of stem measurements with subdiv_id values
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 = "")

 ![4 ha plot split, with stem
locations](https://johngodlee.xyz/img_full/split/split.png)

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)))]
     }))
   }