8wDlpd.png
8wDFp9.png
8wDEOx.png
8wDMfH.png
8wDKte.png

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

Madhura D 1月前

20 0

我正在使用 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)
  • 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

  • 我不相信您可以 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

  • 怎么回事?

    这可能是浮点运算的一个缺点,如下所述: 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')
    

    enter image description here

    st_touches(line_n, point_n)
    #Sparse geometry binary predicate list of length 1, where the predicate was #`touches'
    # 1: (empty)
    
返回
作者最近主题: