Open DOSull opened 4 days ago
Please add set.seed
to any example usong the RNG. What precision settings do you prefer to apply?
Point taken about set.seed
- although in this case, part of the point is the variability. When I set precision I usually go with 1e8
as that's a number that I think you've said in the past has some previous currency! However, in this example I have to set it much more crudely than that st_set_precision(1)
! for the outcome to be consistently what I would hope/expect.
I know that line intersection is notoriously challenging, but that seems rather extreme!
1e8
in this context is the inverse, so:
> formatC(1/1e8, format="f", digits=8)
[1] "0.00000001"
and so on. Precision means converting the coordinates from double precision to an integer grid by multiplying by say 1e8
and checking for integer equality on that grid. sf
by default uses floating point double precision rather than the integer grid (which was the default in defunct rgeos
). So any jitter greater than the chosen precision will lead to the point being seen as not equal.
> set.seed(1)
> sd <- 0.01
> lines <- st_multilinestring(list(matrix(c(-1, 1, -1, 1) +
+ rnorm(4, 0, sd), 2, 2),
+ matrix(c(-1, 1, 1, -1) +
+ rnorm(4, 0, sd), 2, 2))) |>
+ st_sfc() |> st_cast("LINESTRING")
> pt <- st_intersection(lines[1], lines[2]) |> st_sfc()
>
> st_intersects(pt, lines)
Sparse geometry binary predicate list of length 1, where the predicate
was `intersects'
1: (empty)
> st_precision(lines) <- 1e8
> pt <- st_intersection(lines[1], lines[2]) |> st_sfc()
> st_intersects(pt, lines)
Sparse geometry binary predicate list of length 1, where the predicate
was `intersects'
1: (empty)
> st_precision(lines) <- 1e2
> pt <- st_intersection(lines[1], lines[2]) |> st_sfc()
> st_intersects(pt, lines)
Sparse geometry binary predicate list of length 1, where the predicate
was `intersects'
1: 1, 2
I know that precision is set in the opposite sense to how one might expect i.e., using a 'scale' multiplier, per the JTS documentation at https://locationtech.github.io/jts/javadoc/org/locationtech/jts/geom/PrecisionModel.html linked from the st_as_binary
help. I made that mistake a while back in a different issue #1855!
But the above is not what I get...
> set.seed(1)
> sd <- 0.1
> lines <- st_multilinestring(list(matrix(c(-1, 1, -1, 1) +
+ rnorm(4, 0, sd), 2, 2),
+ matrix(c(-1, 1, 1, -1) +
+ rnorm(4, 0, sd), 2, 2))) |>
+ st_sfc() |> st_cast("LINESTRING")
>
> st_precision(lines) <- 1e8
> pt <- st_intersection(lines[1], lines[2]) |> st_sfc()
> st_intersects(pt, lines)
Sparse geometry binary predicate list of length 1, where the predicate was `intersects'
1: (empty)
>
> st_precision(lines) <- 1e2
> pt <- st_intersection(lines[1], lines[2])
> st_intersects(pt, lines)
Sparse geometry binary predicate list of length 1, where the predicate was `intersects'
1: (empty)
In any case... I'm not jittering the point, I'm changing the ends of the lines, and they are intersecting to yield the point, which by definition ought therefore to intersect the lines that generated it? It also seems wrong that to get this result we need a grid as crude as 0.01 resolution (and on my machine I have to go to 0.1, i.e. st_precision(...) <- 10
, and even then I frequently get 'misses'!
> set.seed(2)
> sd <- 0.1
> lines <- st_multilinestring(list(matrix(c(-1, 1, -1, 1) +
+ rnorm(4, 0, sd), 2, 2),
+ matrix(c(-1, 1, 1, -1) +
+ rnorm(4, 0, sd), 2, 2))) |>
+ st_sfc() |> st_cast("LINESTRING")
>
> st_precision(lines) <- 1e8
> pt <- st_intersection(lines[1], lines[2]) |> st_sfc()
> st_intersects(pt, lines)
Sparse geometry binary predicate list of length 1, where the predicate was `intersects'
1: (empty)
>
> st_precision(lines) <- 1e1
> pt <- st_intersection(lines[1], lines[2])
> st_intersects(pt, lines)
Sparse geometry binary predicate list of length 1, where the predicate was `intersects'
1: (empty)
Also worth noting how apparent the poor resolution is when these objects are plotted. The plot below is for the objects generated with set.seed(2)
and st_precision(lines) <- 1e1
on my machine:
> plot(lines)
> plot(pt, add = TRUE)
For what it's worth geos::geos_intersects
yields the same incorrect result on my machine.
> geos::geos_intersects(pt, lines)
[1] FALSE FALSE
Is there some way in which this could be machine-dependent? MacOS ARM architecture here.
How was sf
installed, CRAN binary or otherwise?
Trying on an aarch64 with R 4.4.1 and sf
1.0-16, I get your results, so a double-precision 64-bit/extended precision 80-bit difference most likely. My GEOS on Intel (Fedora 40) is newer too, but not much. What you are "jittering" are the endpoints of the lines, so indeed something is going on here. @paleolimbot any thoughts?
I would have installed from CRAN. I'm glad to hear that I am probably not missing something obvious!
Also... this:
> st_distance(lines, pt)
[,1]
[1,] 8.21824e-17
[2,] 0.00000e+00
where I get some variation on 0
running many times with different random line end offsets, but the intersections still come up empty.
It looks as though the lines would need to be noded at the intersection:
> set.seed(1)
> lines <- st_multilinestring(list(matrix(c(-1, 1, -1, 1) +
+ rnorm(4, 0, sd), 2, 2),
+ matrix(c(-1, 1, 1, -1) +
+ rnorm(4, 0, sd), 2, 2))) |>
+ st_sfc() |> st_cast("LINESTRING")
> pt <- st_intersection(lines[1], lines[2]) |> st_sfc()
> st_intersects(pt, lines)
Sparse geometry binary predicate list of length 1, where the predicate
was `intersects'
1: (empty)
> st_intersects(pt, st_union(lines))
Sparse geometry binary predicate list of length 1, where the predicate
was `intersects'
1: 1
> st_coordinates(lines)
X Y L1
[1,] -1.0062645 -1.0083563 1
[2,] 1.0018364 1.0159528 1
[3,] -0.9967049 1.0048743 2
[4,] 0.9917953 -0.9926168 2
> st_coordinates(st_union(lines))
X Y L1 L2
[1,] -1.006264538 -1.008356286 1 1
[2,] -0.001176252 0.004844437 1 1
[3,] -0.001176252 0.004844437 2 1
[4,] 1.001836433 1.015952808 2 1
[5,] -0.996704922 1.004874291 3 1
[6,] -0.001176252 0.004844437 3 1
[7,] -0.001176252 0.004844437 4 1
[8,] 0.991795316 -0.992616753 4 1
> st_coordinates(pt)
X Y
[1,] -0.001176252 0.004844437
> st_length(st_nearest_points(pt, lines))
[1] 1.612109e-16 6.930893e-17
> st_length(st_nearest_points(pt, st_union(lines)))
[1] 0
See also https://github.com/libgeos/geos/issues/585:
> st_is_within_distance(pt, lines, .Machine$double.eps)
Sparse geometry binary predicate list of length 1, where the predicate
was `is_within_distance'
1: 1, 2
So this is obviously the 'feature-not-a-bug' that prevents me using lwgeom::st_split
to successfully break linestrings at all their points of intersection. I've written code that works around it by inserting all the intersection points into the linestrings but it's not practical because it's slow! Frustrating!
Describe the bug A point I generate by intersecting two linestrings fails expected
st_intersects
and does not yield expectedst_relate
pattern.To Reproduce
So the point resulting from intersecting two lines does not intersect them, and is not related to them as expected (
"0FFFFF102"
). Inspection of the objects shows that the point is at the intersection of the two lines close to (0, 0).Additional context If I make lines running from (-1 -1) to (1 1) and (-1 1) to (1 -1) exactly (by setting
sd <- 0
) then they intersect at (0 0) and that point satisfies expected spatial predicates against the generating intersecting lines, i.e. thest_relate
pattern is"0FFFFF102"
.I have experimented with setting precision on the geometries but have only obtained expected results reliably when precision is set very crudely. If the problem is deemed 'upstream' in GEOS, as I suspect it might be (
geos::geos_relate
produces the same results on these geometries assf
) what might be a viable workaround insf
?I initially encountered this problem working with
lwgeom::st_split()
trying to split lines by their intersection points, and getting many fewer output linestrings than expected. On looking more closely this was what I ran into.