R st_nearest_point 返回不接触第二个几何体的线串-r-IT问答社区-解决你的IT疑问

R st_nearest_point 返回不接触第二个几何体的线串

Madhura D 1月前

我正在使用 st_nearest_point 函数来查找线上距离附近点最近的点。我不确定问题是什么，但我无法让函数返回带有任何点的线串...

我正在使用该 st_nearest_point 函数查找线上与附近点最近的点。我不确定问题是什么，但我无法让函数返回 linestring 任何点接触线的点。我需要使用 4326 crs，因为这是 leaflet 地图使用的，在我的实际应用中，我的起点是 map_click 捕获。 st_nearest_point 返回一个 linestring 接近但不接触线的点。怎么回事？请参见下面的示例。

library(sf)
library(dplyr)
library(mapview)

#Creating some dummy data to test with
# Create a line object (LINESTRING)
line <- st_sfc(st_linestring(matrix(c(-122.4, 37.7, -122.5, 37.8), ncol = 2, byrow = TRUE)), crs = 4326)
# Create a point object (POINT) that is close to but does not touch the line
point <- st_sfc(st_point(c(-122.45, 37.78)), crs = 4326)
# Create sf objects
line_sf <- st_sf(geometry = line)
point_sf <- st_sf(geometry = point)


#Plotting the objects so you can see they don't touch
plot(st_geometry(line_sf), col = 'blue', lwd = 2)
plot(st_geometry(point_sf), col = 'red', add = TRUE)

#Confirming with st_touches
st_touches(line_sf,point_sf)


#Lets find the closest point on the line and use that

geo <- st_nearest_points(point_sf, line_sf)
lngs = geo |> st_length() # Next few lines are redundant here but important if there are multiple lines
nearestp = min(lngs)
my_linestring <- geo[nearestp == lngs]
new_point <- st_cast(my_linestring, 'POINT')[c(FALSE,TRUE)] |> st_as_sf()


#Now lets plot this new point   
plot(st_geometry(line_sf), col = 'blue', lwd = 2)
plot(st_geometry(point_sf), col = 'red', add = TRUE)
plot(st_geometry(new_point), col = 'green', add = TRUE)

#Success? Not so fast lets confirm with st_touches

st_touches(line_sf,new_point)

#Lets use mapview to see if we can see by zooming in the farthest we can

mapview(line_sf) + mapview(new_point, col.region  = "green") + mapview(point_sf, col.region  = "red")

#Looks good but if you zoom in you will see the points indeed do not touch.

帖子版权声明 1、本帖标题：R st_nearest_point 返回不接触第二个几何体的线串
本站网址：http://xjnalaquan.com/
2、本网站的资源部分来源于网络，如有侵权，请联系站长进行删除处理。
3、会员发帖仅代表会员个人观点，并不代表本站赞同其观点和对其真实性负责。
4、本站一律禁止以任何方式发布或转载任何违法的相关信息，访客发现请向站长举报
5、站长邮箱：yeweds@126.com 除非注明，本帖由Madhura D在本站《r》版块原创发布，转载请注明出处！

最新回复 (0)

最新倒序只看楼主

啊鹿Dizzyi 1月前 0 只看Ta

引用 2楼
```
library(sf)
#> Linking to GEOS 3.12.0, GDAL 3.7.1, PROJ 9.2.1; sf_use_s2() is TRUE
```
确保 sf_use_s2() 设置为 TRUE
```
sf::sf_use_s2(TRUE)

line <- st_sfc(st_linestring(matrix(c(-122.4, 37.7, -122.5, 37.8), ncol = 2, byrow = TRUE)), crs = 4326)
point <- st_sfc(st_point(c(-122.45, 37.78)), crs = 4326)
line_sf <- st_sf(geometry = line)
point_sf <- st_sf(geometry = point)

geo <- st_nearest_points(point_sf, line_sf) |>
  sf::st_sf()
```
在这种情况下 st_nearest_points() 和 point_sf 之间的最短线串 line_sf ，无需对其进行排序：
```
geo
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -122.4684 ymin: 37.76848 xmax: -122.45 ymax: 37.78
#> Geodetic CRS:  WGS 84
#>   st_nearest_points.point_sf..line_sf.
#> 1       LINESTRING (-122.45 37.78, ...
```
两条线相交：
```
sf::st_intersects(line_sf, geo)
#> Sparse geometry binary predicate list of length 1, where the predicate
#> was `intersects'
#>  1: 1
```
并且请求 new_point 本身就是一个交叉点：
```
new_point <- st_intersection(line_sf, geo)

tmap::tmap_mode("view")
#> tmap mode set to 'view'
tmap::tm_shape(line_sf) +
  tmap::tm_lines(col = "blue") +
  tmap::tm_shape(point_sf) +
  tmap::tm_dots(size = 0.5) +
  tmap::tm_shape(geo) +
  tmap::tm_lines(col = "green") +
  tmap::tm_shape(new_point) +
  tmap::tm_dots(col = "green", size = 1, shape = 12)
```
但是，在将经度和纬度投影到 X、Y 坐标时，投影 line_sf 为起点和终点之间的直线，而不是曲线。因此 new_point 放置的位置有点偏离。

我会将坐标转换为任何投影的 CRS，进行计算，如果需要，再将其转换回纬度/经度。

^{Created on 2024-09-14 with reprex v2.1.0}
ababo 1月前 0 只看Ta

引用 3楼
我不相信您可以 st_touches() 在这种确切的情况下使用它，因为 \'touches\' 谓词的定义方式（ https://en.wikipedia.org/wiki/DE-9IM ）：
- A 接触 B：它们至少有一个共同点，但它们的内部不相交。
- 线串
  - interior: points that are left when the boundary points are removed.
  - boundary: two end points
换句话说，点只能接触线的端点。或许可以改为评估相交？
```
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
l  <- st_as_sfc("LINESTRING (0 0, 1 1, 2 0)")
p0 <- st_as_sfc("POINT (0 0)")
p1 <- st_as_sfc("POINT (1 1)")

              # 1, 2G, 3R 
st_touches(   c(l, p0, p1), sparse = FALSE)
#>       [,1]  [,2]  [,3]
#> [1,] FALSE  TRUE FALSE
#> [2,]  TRUE FALSE FALSE
#> [3,] FALSE FALSE FALSE
# only p0 touches
st_intersects(c(l, p0, p1), sparse = FALSE)
#>      [,1]  [,2]  [,3]
#> [1,] TRUE  TRUE  TRUE
#> [2,] TRUE  TRUE FALSE
#> [3,] TRUE FALSE  TRUE
# both points intersect
plot(c(l, p0, p1), col = c("darkblue", "darkgreen", "darkred"), pch = 16, cex = 2, axes = TRUE)
```
^{Created on 2024-09-14 with reprex v2.1.1}
Baptiste 1月前 0 只看Ta

引用 4楼
怎么回事？

这可能是浮点运算的一个缺点，如下所述： https://github.com/r-spatial/sf/issues/1503

有一个解决方法，使用 sfnetworks https://luukvdmeer.github.io/sfnetworks/articles/sfn03_join_filter.html#blending-points-into-a-network

这就是下面代码中采用的方法。虽然点最终位于线上，但 st_touches() 仍然返回 false 或空。
```
library(sfnetworks)

#blend the line & points into an sf_newtork
blended <- st_network_blend(as_sfnetwork(line_sf), new_point)
#select the edges & return a linestring
line_n <- blended %>% activate('edges') %>% st_as_sf() %>% st_union() %>% st_as_sf() %>% st_cast('LINESTRING')
# select the nodes & return as points
points_n <- blended %>% activate("nodes") %>% st_as_sf() %>% st_combine() %>% st_as_sf() %>% st_cast('POINT')
point_n <- points_n[3,] #just picking the right point (not the end nodes)

mapview(line_n) + mapview(points_n, col.region = 'green') + mapview(point_n, col.region = 'red')
```
```
st_touches(line_n, point_n)
#Sparse geometry binary predicate list of length 1, where the predicate was #`touches'
# 1: (empty)
```