st_union and floating point coordinates

[DOSull]

Apologies in advance for how long this post became, when it is really just a non-urgent request for clarification

Not so much an issue as a question about how floating point affects st_union (and I guess other spatial operations) in sf. Normally for 'dissolve' or union, I would use (and have used) the group_by() %>% summarise() pipeline, but I'm working on a project where I am making a lot of sf geometries from scratch, translating, rotating, intersecting, unioning them and so on.

Today I came across a situation where st_union fails to dissolve two simple polygons that (at least mathematically!) are touching. Here's some code that reproduces the problem.

# make a square from 4 points at 45, 135 ... etc around 0, 0
a <- seq(1, 7, 2) * pi / 4
a <- c(a, a[1]) # append a copy of the first corner to close the polygon

# make a polygon
p1 <- st_polygon(list(cbind(cos(a), sin(a))))
# make a copy shifted by the width of the first one
bb <- st_bbox(p1)
p2 <- p1 + c(bb$xmax - bb$xmin, 0)

Now when I union these polygons I get a MULTIPOLYGON which when plotted is clearly two polygons

> st_union(p1, p2)
MULTIPOLYGON (((0.7071068 0.7071068, -0.7071068 0.7071068, -0.7071068 -0.7071068, 0.7071068 -0.7071068, 0.7071068 0.7071068)), ((2.12132 0.7071068, 0.7071068 0.7071068, 0.7071068 -0.7071068, 2.12132 -0.7071068, 2.12132 0.7071068)))

For what it's worth integer coordinate squares dissolve as expected.

> p3 <- st_polygon(list(matrix(c(-1, -1, 1, -1, 1, 1, -1, 1, -1, -1), 5, 2, byrow = 2)))
> p4 <- p3 + c(2, 0)
> st_union(p3, p4)
POLYGON ((1 -1, -1 -1, -1 1, 1 1, 3 1, 3 -1, 1 -1))

So is this purely a floating point numbers issue? Inspection of the coordinates of corresponding points in p1 and p2 reveals they do indeed not overlap:

> as.matrix(p1)[2:3, ] == as.matrix(p3)[c(1, 4),]
      [,1]  [,2]
[1,] FALSE FALSE
[2,] FALSE FALSE

but... well... the differences in the coordinates are ~1e-16, And for what it's worth the same operation in rgeos::gUnion dissolves the two squares as I would hope.

> rgeos::gUnion(as(p1, "Spatial"), as(p2, "Spatial"))
class       : SpatialPolygons 
features    : 1 
extent      : -0.7071068, 2.12132, -0.7071068, 0.7071068  (xmin, xmax, ymin, ymax)
crs         : NA 

I can't find a way to change the tolerance that st_union applies. I poked around with st_set_precision() a bit but this didn't seem to get me any closer to fixing my problem. I've experienced presumably related issues when for example I close the square above by generating coordinates at 45, 135... 405 degrees where the error that the polygon is not closed is raised. Presumably because:

> cos(pi/4) - cos(9*pi/4)
[1] -1.110223e-16

I get it that floating point is unsafe for high-precision operations like these, but is there no tolerance setting I can use to avoid having to round coordinates to arbitrary precision as I go along? Or why does sf not apply the machine epsilon?

> .Machine$double.eps
[1] 2.220446e-16

I know I can use st_snap

> st_union(st_snap(p1, p2, 1e-10), p2)
POLYGON ((0.7071068 -0.7071068, -0.7071068 -0.7071068, -0.7071068 0.7071068, 0.7071068 0.7071068, 2.12132 0.7071068, 2.12132 -0.7071068, 0.7071068 -0.7071068))

but it's not obvious to me how to use this to union a set of polygons (keep in mind I am making and manipulating geometries not working with datasets much of the time, so the group_by() %>% summarise() option is not in play). I've had fewer problems when dealing with projected coordinates in sf datasets, so I wonder if there is some under-the-hood magic happening in such cases? I don't usually work with sf for handling 'pure geometry` in this way, so perhaps my inexperience is showing and there is a switch I can flip to fix things?

[edzer]

rgeos sets a default precision value of 1e8, sf doesn't do that. All geometry operations in sf are carried out by GEOS.

[rsbivand]

p1s <- st_as_sfc(list(p1))
p2s <- st_as_sfc(list(p2))
st_precision(p1s)
st_precision(p2s)
st_precision(p1s) <- 1e14
st_precision(p2s) <- 1e14
st_union(p1s, p2s)

seems to work; using 1e8 would ensure rgeos-compatible output, I think.

[DOSull]

This is super helpful!

RTFM but properly! I was setting the precision with a negative exponent... having somehow or other misread or misinterpreted the details. My bad.

It would be nice if there was a package setting like sf_use_s2(), that set precision for subsequent operations, rather than needing to package sfgs up as sfcs and apply it case by case.

I have a hard time imagining situations where (say) 1e8 would not be a sensible default precision (see: https://xkcd.com/2170/), but as you note this is a GEOS thing not sf as such. But really, unless someone is using a CRS with 10km intervals as the units or something (and then you'd have other problems) something in the 1e8 to 1e10 range or even 1e14 seems like a much better choice than floating point with all its weirdnesses.

But I can completely understand if a package option/switch doesn't make sense for general architectural reasons of some kind. My current use-case is rather off the wall.

[edzer]

I like the idea of a globally settable precision value different from zero; feel free to raise it as an issue, or propose a PR.