R/aggregate_polys.R

In mapscanner: Print Maps, Draw on Them, Scan Them Back in

#' Aggregate disparate polygons
#'
#' Planar partition from disparate polygon inputs. Overlaps aggregate to `n`.
#'
#' Input is a single simple features polygon data frame. No attribute data is
#' considered.
#' @param p input (multi-)polygons (assumed to be overlapping)
#' @return Set of \pkg{sf}-format polygons with additional column, `n`, denoting
#' number of overlaps contributing to each of the resultant polygons.
#' @export
#' @examples
#' g <- sf::st_sfc(list(sf::st_point(cbind(0, 0)),
#'                      sf::st_point(cbind(0, 1)),
#'                      sf::st_point(cbind(1, 0))))
#' pts <- sf::st_sf(a = 1:3,  geometry = g)
#' overlapping_polys <- sf::st_buffer(pts, 0.75)
#'
#' ## decompose and count space-filling from overlapping polygons
#' x <- ms_aggregate_polys(overlapping_polys)
#' plot(x)
#' \dontrun{
#' library(ggplot2)
#' ggplot(x) + geom_sf() + facet_wrap(~n)
#' }
#'
#' library(sf)
#' set.seed(6)
#' pts <- expand.grid (x = 1:8, y = 1:10) %>% st_as_sf (coords = c("x", "y"))
#' xsf <- sf::st_buffer (pts, runif (nrow (pts), 0.2, 1.5))
#' \dontrun{
#' out <- ms_aggregate_polys (xsf)
#' }
ms_aggregate_polys <- function (p) {

    p <- sf::st_intersection (p)

    n <- seq (max (p$n.overlaps))
    g <- lapply (n, function (i) {
                     gi <- p$geometry [p$n.overlaps >= i]
                     out <- sf::st_union (gi)
                     if (sf::st_geometry_type (out) == "GEOMETRYCOLLECTION") {
                         out <- sf::st_collection_extract (out) %>%
                             sf::st_union ()
                     }
                     return (out)
    })

    sf::st_sf (n = n, geometry = do.call (c, g))
}