luukvdmeer / sfnetworks

Tidy Geospatial Networks in R
https://luukvdmeer.github.io/sfnetworks/
Other
347 stars 20 forks source link

Blending fails when point is very close to its nearest edge #98

Closed luukvdmeer closed 3 years ago

luukvdmeer commented 3 years ago

Describe the bug When trying to blend a point that is extremely close to its nearest edge, st_network_blend fails.

Reproducible example

library(sf)
#> Linking to GEOS 3.7.1, GDAL 2.2.2, PROJ 4.9.2
library(sfnetworks)

p1 = st_point(c(0.532442, 1.95422))
p2 = st_point(c(0.53236, 1.95377))
p3 = st_point(c(0.53209, 1.95328))
l1 = st_sfc(st_linestring(c(p1, p2, p3)))

p4 = st_point(c(0.53209, 1.95345))
p5 = st_point(c(0.53245, 1.95345))
l2 = st_sfc(st_linestring(c(p4, p5)))

# The two lines share an intersection point.
st_intersection(l1, l2)
#> Geometry set for 1 feature 
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 0.5321837 ymin: 1.95345 xmax: 0.5321837 ymax: 1.95345
#> CRS:            NA
#> POINT (0.5321837 1.95345)

# But this intersection point does not intersects the lines themselves!
st_intersects(l1, st_intersection(l1, l2), sparse = FALSE)
#>       [,1]
#> [1,] FALSE
st_intersects(l2, st_intersection(l1, l2), sparse = FALSE)
#>      [,1]
#> [1,] TRUE

# The intersection point is instead located a tiny bit next to the line.
st_distance(l1, st_intersection(l1, l2))
#>              [,1]
#> [1,] 4.310191e-17

# When we try to blend this point it fails.
net = as_sfnetwork(l1)
p = st_intersection(l1, l2)
st_network_blend(net, p)
#> Warning: st_network_blend assumes attributes are constant over geometries
#> Error in valid_numeric_matrix(x): !anyNA(x) is not TRUE

# The reason:
# Because it is located so close to the line.
# The 'st_nearest_point' line therefore has the same start and endpoint.
# This is an invalid geometry and extending the line fails.
st_nearest_points(p, l1)
#> Geometry set for 1 feature 
#> geometry type:  LINESTRING
#> dimension:      XY
#> bbox:           xmin: 0.5321837 ymin: 1.95345 xmax: 0.5321837 ymax: 1.95345
#> CRS:            NA
#> LINESTRING (0.5321837 1.95345, 0.5321837 1.95345)

Created on 2020-12-14 by the reprex package (v0.3.0)

Expected behavior Blend the point no matter the distance to its nearest edge.

luukvdmeer commented 3 years ago

This is solved now with the new st_network_blend implementation:

library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.0, PROJ 7.2.0
library(sfnetworks)

p1 = st_point(c(0.532442, 1.95422))
p2 = st_point(c(0.53236, 1.95377))
p3 = st_point(c(0.53209, 1.95328))
l1 = st_sfc(st_linestring(c(p1, p2, p3)))

p4 = st_point(c(0.53209, 1.95345))
p5 = st_point(c(0.53245, 1.95345))
l2 = st_sfc(st_linestring(c(p4, p5)))

net = as_sfnetwork(l1)
p = st_intersection(l1, l2)
st_network_blend(net, p)
#> # A sfnetwork with 3 nodes and 2 edges
#> #
#> # CRS:  NA 
#> #
#> # A rooted tree with spatially explicit edges
#> #
#> # Node Data:     3 x 1 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: 0.53209 ymin: 1.95328 xmax: 0.532442 ymax: 1.95422
#>                     x
#>               <POINT>
#> 1  (0.532442 1.95422)
#> 2   (0.53209 1.95328)
#> 3 (0.5321837 1.95345)
#> #
#> # Edge Data:     2 x 3
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: 0.53209 ymin: 1.95328 xmax: 0.532442 ymax: 1.95422
#>    from    to                                                      x
#>   <int> <int>                                           <LINESTRING>
#> 1     1     3 (0.532442 1.95422, 0.53236 1.95377, 0.5321837 1.95345)
#> 2     3     2                   (0.5321837 1.95345, 0.53209 1.95328)

Created on 2021-03-03 by the reprex package (v0.3.0)