3

In order to select/identify the border polygons of a shapefile, I would like to use a function able to select/identify polygon that share a line segment with a source polygon.

With figures:

I have this kind of shapefile:

enter image description here

Using gUnionCascaded from rgeos package, I have a second shapefile with the "contour polygon"

enter image description here

Now I am looking for a function that can select/identify border polygons (shaded on the fig) i.e. polygons of the first shapefile that share a line segment with the polygon of the second shapefile . :

enter image description here

DJack
  • 4,850
  • 3
  • 21
  • 45
  • 1
    Have a look at the very discriminating (but initially hard-to-grock) `rgeos::gRelate()` function. The [DE-9IM Wikipedia page](http://en.wikipedia.org/wiki/DE-9IM) is very helpful in understanding how you could use it, as might be a few SO posts that demonstrate its use. [Here's one example](http://stackoverflow.com/questions/12294185/how-to-create-new-polygons-by-simplifying-from-two-spatialpolygonsdataframe-obje/12327602#12327602) of `rgeos::gRelate()` in action. – Josh O'Brien Nov 05 '13 at 08:57
  • 1
    In particular, you are wanting the component polygons that share borders with the exterior of the large polygon, and whose interiors are overlapped by the interior of the large polygon. – Josh O'Brien Nov 05 '13 at 09:05
  • @JoshO'Brien Thx, I ve used your method and it work great. I did'nt know the DE-9IM code, looks like very interesting. – DJack Nov 05 '13 at 10:30
  • @JoshO'Brien "grock" ? eh? I don't think Heinlein himself could grok "grock" :-) – Carl Witthoft Nov 05 '13 at 12:25
  • @CarlWitthoft -- Ha! Does answering while +/- sleep-walking count as an excuse? Now off to repeat 50 times: "glock, grok; glock, grok; glock, grock; ..." ;-) – Josh O'Brien Nov 05 '13 at 14:01

2 Answers2

3

As proposed by Josh O'Brien, I have used the rgeos::gRelate() function. I get 3 DE-9IM cases:

x <- gRelate(shapefile.1, shapefile.2, byid = TRUE)
table(x)
2FF10F212 2FF11F212 2FF1FF212 
       63      2470    174495

The resulted DE-9IM string codes can be interpreted as follow:

1) 2FF1FF212: represent polygons from the first shapefile that don't intersect the border of the polygon of the second shapefile

2) 2FF11F212: represent polygons from the first shapefile that intersect the border of the polygon of the second shapefile with a line

3) 2FF10F212: represent polygons from the first shapefile that intersect the border of the polygon of the second shapefile with a point

The two last cases are my border polygons that I was looking for. I have got their ID with:

poly.border <- which(x %in% c("2FF10F212","2FF11F212"))
DJack
  • 4,850
  • 3
  • 21
  • 45
1

I want to second the above solution, which really helped me a lot. To add a bit more to the DE-9IM solution (see wikipedia page on DE-9IM), see below for application - in this case, how I used it to subset counties which don't border any county that has at least one CAFO/IHO (concentrated animal feeding operations / industrial hog operations) in it. That is, the blue counties below were compared against a super-shape (like the original question asker using gUnionCascaded), and that vector was used to subset the noIHO counties.

IHOblob = gUnionCascaded(hasIHOcounties)
plot(IHOblob)
touch.v = gRelate(noIHO.counties, IHOblob, byid = T)
counties.not.touching<-(touch.v %in% c("FF2FF1212")) #DE-9IM is super cool.

notouchIHO.counties = noIHO.counties[counties.not.touching,]
plot(notouchIHO.counties, co="light blue", add=T)
invisible(text(getSpPPolygonsLabptSlots(notouchIHO.counties),
  labels=as.character(notouchIHO.counties$NAME), cex=0.4))
#wish I had a better way to label polys than above...`

You can see it here and here:

IHO map

IHO blob

Brian Tompsett - 汤莱恩
  • 5,753
  • 72
  • 57
  • 129
Mike Dolan Fliss
  • 217
  • 2
  • 11