14

How do I simplify an sf polygon without introducing gaps and slivers?

With a shapefile, for example, I would use rmapshaper::ms_simplify():

library("pryr")
library("rgdal")
library("rmapshaper")

download.file("https://borders.ukdataservice.ac.uk/ukborders/easy_download/prebuilt/shape/England_gor_2011.zip",
              destfile = "regions.zip")
unzip("regions.zip")
regions <- readOGR(".", "england_gor_2011")
object_size(regions)
# ~13MB

regions <- ms_simplify(regions)
object_size(regions)
# < 1MB

I've tried sf::st_cast() which, from the man pages, states:

Cast geometry to another type: either simplify, or cast explicitly

and:

to argument: character; target type, if missing, simplification is tried; when x is of type sfg (i.e., a single geometry) then to needs to be specified.

When I've left to as missing this hasn't worked as expected though (I knew it was too good to be true!):

library("sf")
regions <- sf::read_sf("england_gor_2011.shp")
object_size(regions)
# ~13MB

regions <- sf::st_cast(regions)
object_size(regions)
# Still 13MB

Currently I'm opening the file with rgdal::readOGR(), simplifying it, saving this, then loading this again with sf.

Is there a better way?


rgeos::gSimplify()

@s.k 's suggestion of rgeos::gSimplify() can do topologically-aware simplifications (i.e. simplifies without creating slivers) when specified with the following arguments:

library("rgeos")
regions_gSimplify <- gSimplify(regions, tol = 0.05, topologyPreserve = TRUE)

gSimplify doesn't preserve the @data frame, though, so we should re-create it:

regions_df <- regions@data
regions_gSimplify <- sp::SpatialPolygonsDataFrame(regions_gSimplify, regions_df)

And this does indeed result in a smaller file size (can tweak tol argument to get it smaller) and I confirmed this hadn't created any slivers by examining it in QGIS.

object_size(regions_gSimplify)
# ~8MB

So although this is a valid alternative to rmapshaper::ms_simplify() I still have the same problem, namely that it doesn't work with sf:

regions_sf <- sf::read_sf("england_gor_2011.shp")
object_size(regions_sf)

regions_gSimplify <- gSimplify(regions_sf, topologyPreserve = TRUE, tol = 0.05)
# Error in gSimplify(regions_sf, topologyPreserve = TRUE, tol = 0.05) : 
# no slot of name "proj4string" for this object of class "sf"

@obrl_soil 's answer can also be applied to gSimplify(), just use it in place of ms_simplify().

  • 1
    Do you have any access to a Douglas–Peucker algorithm? It is widely know for features simplification in GIS world. stackoverflow.com/questions/17217413/… & r-bloggers.com/simplifying-spatial-polygons-in-r – s.k Jun 11 '17 at 21:38
  • 1
    Isn't st_simplify supposed to do that? (didn't use it, yet) – lbusett Jun 12 '17 at 5:43
  • 2
    ooh, I hadn't noticed st_simplify, thanks for pointing it out. I do prefer the algorithm that rmapshaper::ms_simplify defaults to over all the others I've tried so far, but I'll have a play with the new option (update: whoa proceed with caution, preserveTopology = TRUE definitely doesn't work properly yet) – obrl_soil Jun 12 '17 at 10:13
  • 1
    Good to know. What about opening a bug report on this ? – lbusett Jun 12 '17 at 11:14
  • 1
    @obrl_soil It works for tolerances up to about 1000 on the polygons I've used in the question (regions) but beyond that it no longer preserves topology. As it breaks at a certain point I'd say that's not intended behaviour – Phil Jun 12 '17 at 12:44
16

You can cast an sf object to sp, for packages that don't yet support sf - I do this a fair bit for raster/polygon interactions. So you could do:

simplepolys <- rmapshaper::ms_simplify(input = as(sfobj, 'Spatial')) %>%
  st_as_sf()
  • 1
    This technique - casting as a spatial object, simplifying, then re-casting as an sf object - worked perfectly, and can be used with rmapshaper::ms_simplify() or rgeos::gSimplify(). Thanks for the suggestion! – Phil Jun 12 '17 at 9:54
  • Cool cool, just be aware that the topology of interlocking polygons is only truly preserved with rmapshaper's approach. If your input data is all isolated, non-interlocking polygons, you can safely use any of the available simplification algs. – obrl_soil Jun 12 '17 at 10:24
  • I'm accepting this as the answer because the more canonical sf::st_simplify() isn't robust at high tolerances at the time of writing, although obviously this may change. – Phil Jun 12 '17 at 12:34
  • 8
    I am currently working on support for sf objects in rmapshaper. ms_simplify is available for sf objects in the development version. I'd love early testers - if you want to try it out, you can install with devtools::install_github("ateucher/rmapshaper", ref = "sf") – andyteucher Jun 12 '17 at 17:15
  • 6
    As of rmapshaper version 0.3.0, the call to as( , "Spatial") is no longer required. – Luke1018 May 7 '18 at 1:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.