如何使用 st_join() 与 sf 包进行空间连接

Question

这是我一直在努力的一个玩具示例

# Make points
point1 <- c(.5, .5)
point2 <- c(.6, .6)
point3 <- c(3, 3)
mpt <- st_multipoint(rbind(point1, point2, point3))  # create multipoint

# Make polygons
square1 <- rbind(c(0, 0), c(1, 0), c(1,1), c(0, 1), c(0, 0))
square2 <- rbind(c(0, 0), c(2, 0), c(2,2), c(0, 2), c(0, 0))
square3 <- rbind(c(0, 0), c(-1, 0), c(-1,-1), c(0, -1), c(0, 0))
mpol <- st_multipolygon(list(list(square1), list(square2), list(square2)))  # create multipolygon

# Convert to class 'sf'
pts <- st_sf(st_sfc(mpt))
polys <- st_sf(st_sfc(mpol))

# Determine which points fall inside which polygons
st_join(pts, polys, join = st_contains)

最后一行产生

Error in as.data.frame.default(x[[i]], optional = TRUE, stringsAsFactors = stringsAsFactors) : 
  cannot coerce class "c("sfc_MULTIPOINT", "sfc")" to a data.frame

如何进行空间连接以确定哪些点位于哪些多边形内？

Answer 1

我也在努力解决sf 包的功能，所以如果这不正确，我们深表歉意或有更好的方法。我认为这里的一个问题是，如果像在您的示例中那样构建几何图形，您将无法获得您的想法：

> pts
Simple feature collection with 1 feature and 0 fields
geometry type:  MULTIPOINT
dimension:      XY
bbox:           xmin: 0.5 ymin: 0.5 xmax: 3 ymax: 3
epsg (SRID):    NA
proj4string:    NA
                     st_sfc.mpt.
1 MULTIPOINT(0.5 0.5, 0.6 0.6...

> polys
Simple feature collection with 1 feature and 0 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 0 ymin: 0 xmax: 2 ymax: 2
epsg (SRID):    NA
proj4string:    NA
                    st_sfc.mpol.
1 MULTIPOLYGON(((0 0, 1 0, 1 ...

您可以看到pts 和polys 中都只有一个“功能”。这意味着您正在构建一个“多面体”特征（即由 3 个部分构成的面），而不是三个不同的面。积分也是一样。

经过一番挖掘，我发现使用 WKT 表示法构建几何图形的这种不同（并且在我看来更容易）方式：

polys <- st_as_sfc(c("POLYGON((0 0 , 0 1 , 1 1 , 1 0, 0 0))",
                     "POLYGON((0 0 , 0 2 , 2 2 , 2 0, 0 0 ))", 
                     "POLYGON((0 0 , 0 -1 , -1 -1 , -1 0, 0 0))")) %>% 
  st_sf(ID = paste0("poly", 1:3))    

pts <- st_as_sfc(c("POINT(0.5 0.5)",
                   "POINT(0.6 0.6)",
                   "POINT(3 3)")) %>%
  st_sf(ID = paste0("point", 1:3))

> polys
Simple feature collection with 3 features and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -1 ymin: -1 xmax: 2 ymax: 2
epsg (SRID):    NA
proj4string:    NA
     ID                              .
1 poly1 POLYGON((0 0, 0 1, 1 1, 1 0...
2 poly2 POLYGON((0 0, 0 2, 2 2, 2 0...
3 poly3 POLYGON((0 0, 0 -1, -1 -1, ...

> pts
Simple feature collection with 3 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 0.5 ymin: 0.5 xmax: 3 ymax: 3
epsg (SRID):    NA
proj4string:    NA
      ID              .
1 point1 POINT(0.5 0.5)
2 point2 POINT(0.6 0.6)
3 point3     POINT(3 3)

您可以看到现在polys 和pts 都具有三个特征。

我们现在可以使用以下方法找到“交集矩阵”：

# Determine which points fall inside which polygons
pi <- st_contains(polys,pts, sparse = F) %>% 
  as.data.frame() %>% 
  mutate(polys = polys$ID) %>% 
  select(dim(pi)[2],1:dim(pi)[1])
colnames(pi)[2:dim(pi)[2]] = levels(pts$ID)

> pi
  polys point1 point2 point3
1 poly1   TRUE   TRUE  FALSE
2 poly2   TRUE   TRUE  FALSE
3 poly3  FALSE  FALSE  FALSE

意思是（正如 cmets 中的 @symbolixau 所指出的那样）多边形 1 和 2 包含点 1 和 2，而多边形 3 不包含任何点。相反，点 3 不包含在任何多边形中。

HTH。

Answer 2

我看到了不同的输出：

> # Determine which points fall inside which polygons
> st_join(pts, polys, join = st_contains)
Simple feature collection with 1 feature and 0 fields
geometry type:  MULTIPOINT
dimension:      XY
bbox:           xmin: 0.5 ymin: 0.5 xmax: 3 ymax: 3
epsg (SRID):    NA
proj4string:    NA
                        geometry
1 MULTIPOINT(0.5 0.5, 0.6 0.6...

这是最新的 CRAN 版本的 sf 吗？

Answer 3

请注意，原始的多点和多多边形集可以“投射”为点和多边形，而无需创建新对象：

st_contains(polys %>% st_cast("POLYGON"), pts %>% st_cast("POINT"), sparse = F)
#      [,1]  [,2]  [,3]
#[1,]  TRUE  TRUE FALSE
#[2,]  TRUE  TRUE FALSE
#[3,] FALSE FALSE FALSE