14

I am running some geoprocessing tasks in R, in which I am trying to create some polygons for clipping rasters of environmental information. I am buffering somewhat complex polygons, and this leaves small subgeometries that I would like to get rid of. In ArcGIS, I think this would involve converting my polygon from multipart to singlepart (or something along those lines) and then dissolving, but I don't know how to do this in R.

Here's an example that illustrates the problem:

require(maptools)
require(rgeos)

data(wrld_simpl)
wrld_simpl[which(wrld_simpl@data$NAME=='Greece'),]->greece
proj4string(greece)<-CRS('+proj=lonlat +datum=WGS84')
gBuffer(greece,width=0.5)->buf
plot(buf)

What I really want is the outer boundary of the polygon, with nothing else inside. Any ideas?

Pascal
  • 1,590
  • 2
  • 16
  • 35
  • I get this error from your code: `A new CRS was assigned to an object with an existing CRS: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 without reprojecting. For reprojection, use function spTransform in package rgdal` – nico Sep 30 '12 at 18:13
  • That's just a warning because I re-defined the projection on the object so that it would play nice with the gBuffer function. You can ignore that. – Pascal Sep 30 '12 at 18:15
  • Correct me if I'm wrong, but is your question going to boil down to calculating the minimum convex hull of the vertexes of your polygons? – BenBarnes Sep 30 '12 at 19:34
  • The convex hull approach had crossed my mind, but that really generalizes the polygon alot. In this case, it works pretty well, but there are polygon shapes where the MCP will really change the end result. So, you're right! But I was hoping for another approach. – Pascal Sep 30 '12 at 20:10

2 Answers2

13

If you just want to get the one ring that forms the boundary of your buffer, then this:

plot(SpatialPolygons(list(Polygons(list(buf@polygons[[1]]@Polygons[[1]]),ID=1))),lwd=2)

is a very ad-hoc way of doing it (and plotting it) for your case.

What you really really want is to get all the rings with ringDir=1, since the rest will be holes. You need all the rings because your buffer might still be two disconnected islands.

outerRings = Filter(function(f){f@ringDir==1},buf@polygons[[1]]@Polygons)
outerBounds = SpatialPolygons(list(Polygons(outerRings,ID=1)))
plot(outerBounds)

might do the trick... Try it with width=0.1 and you'll see it work with multiple islands, but still removing a hole.

MichaelChirico
  • 33,841
  • 14
  • 113
  • 198
Spacedman
  • 92,590
  • 12
  • 140
  • 224
  • This seems to be exactly what I was hoping for. Thank you! I really need to learn more about the structure of these polygon objects... – Pascal Sep 30 '12 at 23:59
  • 1
    There should be some functions that make it easier to work with the inner guts of sp-objects, but I can never find them... str(buf) is your friend. – Spacedman Oct 01 '12 at 06:46
  • I couldn't find any documentation on `ringDir`, any pointers? – MichaelChirico Nov 16 '16 at 01:09
  • the second thing works! Awesome. Isn't there a standard way to do this? Using some rgdal function? – Tomas Aug 10 '18 at 14:57
3

If you want the convex hull that will fit Greece, you can use the gConvexHull function in the rgeos package. Note that this is not necessarily the approach to take if you are dealing with polygons with holes in them, as I thought was the case from the question's title. However, from your example, it looks like the below approach will get you where you want.

myCH <- gConvexHull(greece)

plot(myCH)

which will produce something like

enter image description here

And to check that everything fits,

plot(myCH)
plot(greece,add=TRUE)

enter image description here

BenBarnes
  • 19,114
  • 6
  • 56
  • 74