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)