在 R 中使用两个几何对 sf 进行几何运算

Question

我有一个包含两个几何列的 sf 数据集。这是它的样子：

> trip_geo

   dstid sourceid                    dest_geom                    source_geom
1    1        1   MULTIPOLYGON (((-2.607985 5... MULTIPOLYGON (((-2.607985 5...
2    1        2   MULTIPOLYGON (((-2.607985 5... MULTIPOLYGON (((-2.57022 51...
3    1        3   MULTIPOLYGON (((-2.607985 5... MULTIPOLYGON (((-2.593213 5...
4    1        4   MULTIPOLYGON (((-2.607985 5... MULTIPOLYGON (((-2.608686 5...
5    1        5   MULTIPOLYGON (((-2.607985 5... MULTIPOLYGON (((-2.512852 5...

活动几何是dest_geom。

每一行对应于社区之间的一次旅行。对于每次旅行，我都想知道哪些街区会阻碍旅行。也就是说，如果要在每一行的 source_geom 和 dest_geom 之间画一条直线，哪些几何图形会接触这条直线？我想获取该行的所有接触几何图形，然后将它们合并。

我有另一个数据集，其几何对应于每个 id：

> id_geo

   dstid                       geometry
1      1 MULTIPOLYGON (((-2.607985 5...
2      2 MULTIPOLYGON (((-2.57022 51...
3      3 MULTIPOLYGON (((-2.593213 5...
4      4 MULTIPOLYGON (((-2.608686 5...
5      5 MULTIPOLYGON (((-2.512852 5...

我想第一步是在每次行程的 source_geom 和 dest_geom 的质心之间定义一条线。然后，创建一个 sf，其中几何列包含接触线的多边形列表（我不知道是否可以在 sf 的一列中有多个几何图形）。然后，合并同一行/列表中包含的几何图形。

我不认为我解决问题的方式是正确的，因为据我所知，无法对 sf 的两个几何图形执行操作，例如定义一条线。此外，我不知道如何将列表集成到数据框/sf。

如果您能提出一种更现实的解决问题的方法，我将不胜感激。

Answer 1

在一个 sf 数据框中有两个几何图形并没有错，但是当它与采用隐式几何图形的函数一起使用时，其中只有一个可以是几何图形，例如 st_centroid(foo) 它获取集合几何图形的质心。

对于其他几何，您可以处理该命名列，例如st_centroid(foo$source_geom)。

对于具有两个多边形几何形状的数据框nc，您可以通过首先合并这些点以形成 MULTIPOINT 然后将它们转换为 LINESTRING 来计算质心之间的线。例如，对于第一行：

> st_cast(st_union(c(st_centroid(nc$source_geom[1]), st_centroid(nc$dest_geom[1]))),"LINESTRING")
Geometry set for 1 feature 
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: -81.49826 ymin: 36.42228 xmax: -77.41056 ymax: 36.4314
epsg (SRID):    NA
proj4string:    NA
LINESTRING (-81.49826 36.4314, -77.41056 36.42228)

您必须逐行执行此操作，否则您最终会在整个几何向量上进行操作。

完整示例。执行library(spdep) 和example(poly2nb) 以获得nc.sids。

首先将其缩减为两列和 5 个随机行：

> nc = nc.sids[,c("NAME","FIPS")]
> nc = nc[sample(1:nrow(nc.sids),5),]
> nc
Simple feature collection with 5 features and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -83.73952 ymin: 34.36442 xmax: -78.16968 ymax: 36.54607
epsg (SRID):    NA
proj4string:    NA
         NAME  FIPS                       geometry
82 Cumberland 37051 MULTIPOLYGON (((-78.49929 3...
96     Bladen 37017 MULTIPOLYGON (((-78.2615 34...
13  Granville 37077 MULTIPOLYGON (((-78.74912 3...
78      Macon 37113 MULTIPOLYGON (((-83.10629 3...
14     Person 37145 MULTIPOLYGON (((-78.8068 36...

假设特征 1 转到特征 2、2 到 3 等。创建新几何列：

> nc$dest_geom = nc$geometry[c(2,3,4,5,1)]

现在制作连接质心的线：

> nc$join_geom = st_sfc(sapply(1:nrow(nc),function(i)st_cast(st_union(c(st_centroid(nc$geometry[i]), st_centroid(nc$dest_geom[i]))),"LINESTRING")))

剧情：

> plot(nc$geom)
> plot(nc$join_geom,add=TRUE,lty=2)

Answer 2

@Spacedman 的回复非常适合在感兴趣的多边形之间画线。下面是 (1) 绘制线和 (2) 选择与每条线相交的多边形并合并选择的代码。

#Draw a line between each source and destination
library(sf)
trip_geo$join_geom = st_sfc(sapply(1:nrow(trip_geo),function(j)
+ st_cast(st_union(c(st_centroid(trip_geo$dest_geom[j]),
+ st_centroid(trip_geo$source_geom[j]))),"LINESTRING")), crs=4326)

#For each trip, select the polygons that intersect the line, then union them together
df=data.frame()
for (i in 1:nrow(trip_geo)){
    id_geo$touch=apply(st_intersects(id_geo, trip_geo$join_geom[i]), 1, any)
    touch=subset(id_geo, touch==T)
    touch=st_union(touch)
    df[i,1]=touch #append the unioned polygons to an empty dataframe
  }

#Add the unioned polygons to the original dataset
sf=st_sf(df, crs=4326) 
bound=cbind(trip_geo, sf)

计算需要一些时间（画线大约需要 1 分钟，循环大约需要 1 分钟，有 3,500 个多边形），我确信我的代码可以改进。