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