25

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().

Phil
  • 1,129
  • 2
  • 11
  • 15
  • 1
    Do you have any access to a Douglas–Peucker algorithm? It is widely know for features simplification in GIS world. https://stackoverflow.com/questions/17217413/why-is-rs-implementation-of-the-douglas-peucker-algorithm-so-slow & https://www.r-bloggers.com/simplifying-spatial-polygons-in-r/ – swiss_knight 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 05:43
  • @LoBu Yep, that's it! I can't believe I missed that! Care to answer and I'll accept it, as that's probably the canonical sf answer. Thanks. – Phil Jun 12 '17 at 09:57
  • @s.k See my edited question and thanks for the suggestion – Phil Jun 12 '17 at 09:59
  • 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
  • @obrl_soil preserveTopology = TRUE seemed to work for tol up to about 1000 (of course YMMV), after that things got a bit psychedelic – Phil Jun 12 '17 at 12:33
  • @LoBu I'm honestly not sure if its a true bug or a semantic thing, 'preserve topology' seems to mean different things to different people in this space. It does deserve disussion, though. – obrl_soil Jun 12 '17 at 12:36
  • 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
  • 1
    Yeah, I managed to break it at 25m on a more local-scale dataset. Issue opened here - https://github.com/edzer/sfr/issues/381 – obrl_soil Jun 12 '17 at 13:39

1 Answers1

26

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()
obrl_soil
  • 3,752
  • 15
  • 34
  • 2
    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 09: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
  • 11
    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
  • 10
    As of rmapshaper version 0.3.0, the call to as( , "Spatial") is no longer required. – Luke1018 May 07 '18 at 01:16