💾 Archived View for republic.circumlunar.space › users › johngodlee › posts › 2021-11-14-minimum_rec… captured on 2023-09-08 at 16:19:32. Gemini links have been rewritten to link to archived content

-=-=-=-=-=-=-

Minimum bounding rectangles around sf polygons in R

DATE: 2021-11-14

AUTHOR: John L. Godlee

I wanted to calculate the minimum bounding rectangle around various polygons in a dataset I was working on. The minimum bounding rectangle is the smallest rectangle that totally encompasses the given polygon.

Example of a polygon and a minimum bounding rectangle.

I have been using R and the sf package to work with the polygons. The polygons represent projected tree crown areas. The minimum bounding rectangle is a useful thing to calculate as it can tell you about elongation of the polygon and the direction of the axis of elongation, which in my case might tell me something about wind direction.

I wrote this function which returns the vertex point coordinates, width, length, and long-axis angle of the minimum bounding rectangle around a set of polygons:

#' Return the minimum bounding rectangle around a set of polygons
#'
#' @param x sf object containing only geometries of type \code{POLYGON} or 
#'     \code{MULTIPOLYGON}
#'
#' @return list containing vertex points (\code{ptx}), 
#'     width (\code{width}) length (\code{length}) and angle of bounding 
#'     rectangle of each polygon.
#' 
#' @examples
#' dat <- sf::st_read(system.file("shape/nc.shp", package="sf"))
#' min_box_list <- minBox(dat)
#' min_box_list[[1]]
#'
#' min_box_sf <- do.call(rbind, lapply(min_box_list, function(x) {
#'   pts_sf <- sf::st_as_sf(as.data.frame(x$pts), coords = c("X", "Y"))
#'   sf::st_sf(geometry = sf::st_convex_hull(sf::st_union(pts_sf)), crs = sf::st_crs(dat))
#'   }))
#' plot(sf::st_geometry(min_box_sf), col = NA, border = "red")
#' plot(sf::st_geometry(dat), col = NA, border = "blue", add = TRUE)
#' 
#' @importFrom sf st_is st_coordinates
#' 
#' @export
#' 
minBox <- function(x) {
  stopifnot(all(sf::st_is(x, c("POLYGON", "MULTIPOLYGON"))))

  lapply(sf::st_geometry(x), function(y) {
    x_mat <- sf::st_coordinates(y)[,1:2]

    # Extract convex hull of polygon
    H <- chull(x_mat)
    n <- length(H)
    hull <- x_mat[H, ]

    # Get direction vector
    hDir <- diff(rbind(hull, hull[1,]))
    hLens <- sqrt(rowSums(hDir^2))
    huDir <- diag(1/hLens) %*% hDir
    ouDir <- cbind(-huDir[,2], huDir[,1])

    # Project hull vertices 
    projMat <- rbind(huDir, ouDir) %*% t(hull)

    # Get width and length of bounding rectangle
    rangeH <- matrix(numeric(n*2), ncol = 2)
    rangeO <- matrix(numeric(n*2), ncol = 2)
    widths <- numeric(n)
    lengths <- numeric(n)

    for(i in seq(along=numeric(n))) {
      rangeH[i,] <- range(projMat[i,])
      rangeO[i,] <- range(projMat[n+i,])
      widths[i] <- abs(diff(rangeH[i,]))
      lengths[i] <- abs(diff(rangeO[i,]))
    }

    eMin <- which.min(widths*lengths)
    hProj <- rbind(rangeH[eMin,], 0)
    oProj <- rbind(0, rangeO[eMin,])

    # Move projections to rectangle corners
    hPts <- sweep(hProj, 1, oProj[,1], "+")
    oPts <- sweep(hProj, 1, oProj[,2], "+")

    # Get corner coordinates
    basis <- cbind(huDir[eMin,], ouDir[eMin,])
    hCorn <- basis %*% hPts
    oCorn <- basis %*% oPts
    pts <- t(cbind(hCorn, oCorn[,c(2,1)]))

    # Angle
    dPts <- diff(pts)
    e <- dPts[which.max(rowSums(dPts^2)), ] 
    eUp <- e * sign(e[2])
    deg <- atan2(eUp[2], eUp[1])*180/pi

    return(list(pts = pts, length = lengths[eMin], 
        width = widths[eMin], angle = deg))
  })
}

Here is the output of the example code in the function definition, showing the minimum bounding rectangles of each county in North Carolina, which is a dataset included in the sf package.

Counties in North Carolina and their minimum bounding rectangles.