# Minimum bounding rectangles around sf

## Found at: republic.circumlunar.space:70/~johngodlee/posts/2021-11-14-minimum_rectangle.txt

```TITLE: 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.](https://johngodlee.xyz/img_full/minimum_rectangle/exampl
```
```e.png)
```

```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.](https://johngodlee.xyz/img_full/minimum_rectangle/nc.pn
```
```g)
```