1

I'm trying to identify the common borders of two different polygons using the sf_intersection() function from the sf package.

I tried this simple approach in my data, which comes from a shapefile, but it's not working exactly as I expected.

My data is the shapefile names "zones" from this repository, and this is what I've tried:

library(sf)
library(ggplot2)

zones <- st_read('./Data/zones.shp')
zones$id <- seq(nrow(zones))
borders <- st_intersection(zones, zones)
borders <- borders[borders$id != borders$id.1, ]

ggplot() +
  geom_sf(data = zones, color='red', fill=NA) +
  geom_sf(data = borders, color = 'navy')

The final plot yields this result:

enter image description here

If you look carefully, you'll note that there are some portions of the inner line of the polygons that are not part of the line in borders (they are red and not blue). I don't know why this is happening. Any hint or advice will be much appreciated. Thanks!

Luise
  • 482
  • 1
  • 6
  • 21
  • As @David_O mentions, the boundaries are not actually shared and therefore don't intersect. Try, for example: `ggplot(data = zones) +geom_sf() + xlim(c(768824.7, 770092.5)) + ylim(c(2944198.6, 2944940.3))` you will see a section where the two linestrings diverge. – ander2ed Dec 18 '19 at 22:49
  • 1
    You might try to `st_buffer` your zones object before intersecting. – ander2ed Dec 18 '19 at 22:53
  • The best thing I've ever used for repairing polygon lattices to fix up non-overlapping borders or overlapping slivers is: https://github.com/tudelft3d/pprepair - some assembly required. – Spacedman Dec 19 '19 at 08:09

2 Answers2

3

It is local imprecision in the borders. With most vector data formats, shared POLYGON borders are duplicated in each of the neighbours. It doesn't take very much for slight differences in the coordinates to make the intersection of the two borders incomplete.

That's not a solution, I'm afraid.

This section shows the kind of problem. The view window is:

> par('usr')
[1]  764968.2  765650.8 2945266.2 2945890.9

example gap

That sliver is only about 3 metres wide.

EDIT: Just to add my attempt at a solution using st_snap. This seems to do the trick in some places but not consistently. It doesn't feel like it is working as intended. Also, just to note the projection uses US feet as units, which confused me.

z1 <- st_geometry(zones[1,])
z2 <- st_geometry(zones[2,])

z1 <- st_cast(z1, 'LINESTRING')
z2 <- st_cast(z2, 'LINESTRING')

z1s <- st_snap(z1, z2, 1000)

border <- st_intersection(z1s, z2)

That snap tolerance is way over the top - the gaps between the zone seem to be < 10 feet - but even with this huge tolerance the actual border has missing sections. More oddly, the result has a totally unexpected extension that heads off >6500 feet from the actual intersection.

enter image description here

David_O
  • 1,143
  • 7
  • 16
  • Thanks, that was the first idea that came into my mind when I saw the result. I'll try looking into a possible solution to this issue and will post it if I can find it. – Luise Dec 18 '19 at 22:49
  • I've played around with `st_snap` which my intuition says should snap the vertices of one polygon to another with a given tolerance, but it doesn't quite seem to work as I'm expecting – David_O Dec 18 '19 at 22:54
2

@David_O identifies the issue - the POLYGON borders don't actually touch throughout the shared boundary so st_intersection won't identify them as such.

One workaround may be to st_buffer your zones object before intersecting, although this is admittedly a crude workaround:

borders <- st_intersection(st_buffer(zones, 5), st_buffer(zones, 5))

borders <- borders[borders$id != borders$id.1, ]

ggplot() +
  geom_sf(data = zones, color='red', fill="transparent") +
  geom_sf(data = borders, color = 'navy')
ander2ed
  • 1,318
  • 1
  • 11
  • 19
  • That is a neat fix. A problem is that it will give you back a 10m wide POLYGON surrounding the border rather than a LINESTRING. Admittedly, I suspect you're already getting a Polygon or GeometryCollection. If this is just for display that may well work for you, but if you're after the border as a geometry for analytical purposes this may be an issue. – David_O Dec 18 '19 at 23:01
  • Great! indeed, using `st_buffer()` is a workaround to my issue – Luise Dec 18 '19 at 23:02
  • Hmm, yeah - I think it works okay for simple plotting but isn't very accurate. I think there is an `st_snap` solution. Let me investigate.. – ander2ed Dec 18 '19 at 23:03
  • My attempts with `st_snap` have been returning GEOMETRYCOLLECTION with an odd polygon on one end. Unfortunately I'm not good enough with `sf` objects and the spatial predicates to understand why this is happening or how to solve it. – ander2ed Dec 18 '19 at 23:40