R/despecle_multipolygons.R

In davelovellCARU/carumapr: Handling Anglican map data

## Clean parish map DL 191205

# This function cleans complicated geometries, returning the largest polygon for a give vector of multipolygons.
## Robust against mixed multipolygon / polygon vectors.

despecle_multipolygon =
  function(geomCol,threshold = 1000) {
    ### Save the attributes of the geometry, which would otherwise be lost
    geomCol %>%
      attributes  ->
      oldGeomAttribs

    ### Pluck the BIGGEST element from each list-within-a-list  [`purrr::keep()`]
    ### Put it inside a one element list                      [`list()`]
    ### Turn the list into a polygon                          [`sf::st_polygon()`]
    ### Return a list of such polygons                        [`purrr::map()`]
    ### turn that list into a geometry column                 [`sf::st_sfc()`]

    ##### Also if it's  a MULTIPOLYGON we need to pluck twice

    (purrr::map)(geomCol,
        ~ {
          if ("MULTIPOLYGON" %in% class(.)) {
            # Unlist multipolygons::::::::::::::::::::::::::::::
            unlist(., recursive = FALSE) %>%
              # Manually recast each multipolygon into a list of polygons (counterintuitive because that's literally what they are)
              (purrr::map)( ~ st_polygon(list(.))) ->
              thePolygons

            # Find the area of the biggest polygon in each multipolygon ::::::::::::::::::::::::::
            thePolygons %>%
              (purrr::map_dbl)(sf::st_area) %>%
              max ->
              bigArea

            # Keep the biggest polygon in each MP. (Annoyingly class is lost here, thought I could avoid it)
            goodPolygons =
              thePolygons %>%
              (purrr::keep)( ~ (sf::st_area)(.) > bigArea*(1/threshold))

            #Unlist those 1 element lists and turn the lists inside of them into polygons (wat?)
            biggestPolygon =
              (sf::st_multipolygon)(goodPolygons)

            return(biggestPolygon) # return those polygons

          } else if ("POLYGON" %in% class(.)) {
            return(.)   # Think I should have done this all along, we don't need to change any polygons
          } else
            stop("Element class is not in c(\"MULTIPOLYGON\",\"POLYGON\")")
        }) ->
      newGeom

    ### Change the attributes of new column to match those of the old:
    attributes(newGeom) = oldGeomAttribs

    #### Return the new column in place of the old
    return(newGeom)
  }