R，sf：与多边形边界相交的线，提取这些相交点的坐标

Question

我是SF和stack的新手，希望我的问题足够清晰。 我已经能够创建一组连接1个点和遍布美国的一组点的线。 然后我可以将美国各县读入多边形中。 我的目标是找到并定位我创建的线穿过县界的所有点。

到目前为止，我已经能够从这些点创建线：

points_to_lines <- dt %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% 
  group_by(lineid) %>%
  summarize(do_union = FALSE) %>% lineid
  st_cast("LINESTRING")

这是行的标题

Simple feature collection with 1628 features and 1 field
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: 30.1127 ymin: -91.32484 xmax: 37.23671 ymax: -82.31262
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
# A tibble: 1,628 x 2
   lineid                                 geometry
    <int>                         <LINESTRING [°]>
 1      1  (33.51859 -86.81036, 36.16266 -86.7816)
 2      2 (33.51859 -86.81036, 34.61845 -82.47791)

这是县级数据集的标题。

Reading layer `US_county_1930_conflated' from data source `~/county_gis/1930' using driver `ESRI Shapefile'
Simple feature collection with 3110 features and 18 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -7115608 ymin: -1337505 xmax: 2258244 ymax: 4591848
epsg (SRID):    NA
proj4string:    +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs

我很天真地尝试给它们相同的坐标，然后使用st_intersects进行相交。非稀疏矩阵似乎表明所有线只与一个县相交。

gis1930_p <- st_set_crs(gis1930, 4326) %>% st_transform(4326)
st_intersects(points, gis1930_p, sparse=FALSE)

我也在县域上绘制了线条，但只绘制了美国县域的地图。

plot(gis1930_p[0], reset = FALSE)
plot(points[0], add = TRUE)

非常感谢您的帮助，如果需要我提供更多细节，请告诉我。

- Gianluca

Answer 1

您没有提供数据，因此我将使用以下链接中提供的数据集：https://gis.stackexchange.com/a/236922/32531

您需要的主要函数是st_intersection：

最初的回答：

library(sf)

line_1 <- st_as_sfc("LINESTRING(458.1 768.23, 455.3 413.29, 522.3 325.77, 664.8 282.01, 726.3 121.56)")
poly_1 <- st_as_sfc("MULTIPOLYGON(((402.2 893.03, 664.8 800.65, 611.7 666.13, 368.7 623.99, 215.1 692.06, 402.2 893.03)), ((703.9 366.29, 821.2 244.73, 796.1  25.93, 500.0 137.76, 703.9 366.29)))")
pnts   <- st_intersection(line_1, 
                          st_cast(poly_1, "MULTILINESTRING", group_or_split = FALSE))

plot(poly_1)
plot(line_1, add = TRUE)
plot(pnts, add = TRUE, col = "red", pch = 21)

- jsta

---

非常感谢。这似乎确实是我正在寻找的东西。不幸的是，我无法使其正常工作。当我在县和一条线之间进行st_interesection时，我得到了一个空数据集。 我注意到不同的是，在我的设置中，我有两个地理参考不同的对象，也许？我必须运行st_set_crs使它们相同（例如4326）。奇怪的是，如果我在县上绘制线，线不会出现。 - Gianluca

是的，这听起来像是一个投影问题。你应该使用 st_transform 来改变坐标参考系，而不是 st_set_crs。我使用 mapview 包来调试投影错误。 - jsta

嘿，非常感谢，我在代码中做了两个更改后它起作用了：

1. 我使用了st_transform然后st_set_src(4326)。
2. 我切换到st_as_sf(coords = c("lat","lon"))而不是st_as_sf(coords = c("lon","lat"))！ 此外，关于mapview的提示非常好，帮助我找出第二个错误！！

- Gianluca