如何使用sf方法（如st_intersects()）过滤R语言的简单要素集合？

Question

SF是一个针对tidy语法（例如dyplr和pipes）设计的R-Spatial软件包。

我想在简单特征集合对象上执行简单的空间过滤器。给定一个简单的特征集合，我希望返回符合某些几何条件的所有要素。特别是，我想找到与另一个对象相交的要素。

SF提供了函数st_intersects(x, y, ...)来实现这一点，但我无法将其与dplyr配合使用。

我正在使用R 3.5.2和从github安装的最新sf。

library(tidyverse)
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

# many multipolygons:
nc <- st_read(system.file("shape/nc.shp", package="sf"))

#> Simple feature collection with 100 features and 14 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> epsg (SRID):    4267
#> proj4string:    +proj=longlat +datum=NAD27 +no_defs

# A point in Ashe County:
ash_point <- nc %>% 
  filter(NAME == "Ashe") %>% 
  st_point_on_surface()

# how many counties intersect ash_point? 
nc %>% 
  st_intersects(ash_point, sparse = FALSE) %>% 
  sum()
#> [1] 1

# return the features which intersect ash_point:
nc %>% 
  filter(st_intersects(ash_point, sparse = FALSE)) 
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
#> Simple feature collection with 100 features and 14 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> epsg (SRID):    4267
#> proj4string:    +proj=longlat +datum=NAD27 +no_defs
#> First 10 features:
#>     AREA PERIMETER CNTY_ CNTY_ID        NAME  FIPS FIPSNO CRESS_ID BIR74
#> 1  0.114     1.442  1825    1825        Ashe 37009  37009        5  1091
#> 2  0.061     1.231  1827    1827   Alleghany 37005  37005        3   487
#> 3  0.143     1.630  1828    1828       Surry 37171  37171       86  3188
#> 4  0.070     2.968  1831    1831   Currituck 37053  37053       27   508
#> 5  0.153     2.206  1832    1832 Northampton 37131  37131       66  1421
#> 6  0.097     1.670  1833    1833    Hertford 37091  37091       46  1452
#> 7  0.062     1.547  1834    1834      Camden 37029  37029       15   286
#> 8  0.091     1.284  1835    1835       Gates 37073  37073       37   420
#> 9  0.118     1.421  1836    1836      Warren 37185  37185       93   968
#> 10 0.124     1.428  1837    1837      Stokes 37169  37169       85  1612
#>    SID74 NWBIR74 BIR79 SID79 NWBIR79                       geometry
#> 1      1      10  1364     0      19 MULTIPOLYGON (((-81.47276 3...
#> 2      0      10   542     3      12 MULTIPOLYGON (((-81.23989 3...
#> 3      5     208  3616     6     260 MULTIPOLYGON (((-80.45634 3...
#> 4      1     123   830     2     145 MULTIPOLYGON (((-76.00897 3...
#> 5      9    1066  1606     3    1197 MULTIPOLYGON (((-77.21767 3...
#> 6      7     954  1838     5    1237 MULTIPOLYGON (((-76.74506 3...
#> 7      0     115   350     2     139 MULTIPOLYGON (((-76.00897 3...
#> 8      0     254   594     2     371 MULTIPOLYGON (((-76.56251 3...
#> 9      4     748  1190     2     844 MULTIPOLYGON (((-78.30876 3...
#> 10     1     160  2038     5     176 MULTIPOLYGON (((-80.02567 3...

^{此内容由reprex包（v0.3.0.9000）于2019-07-12创建}

单独使用st_intersects（）可返回正确的逻辑矩阵，但在过滤器中使用时，即使逻辑矩阵为“FALSE”的要素也会返回所有结果。

- jmw

Answer 1

请注意，st_intersection(, sparse = TRUE) 返回一个逻辑型的 matrix，而 filter 需要一个向量。我们可以通过对矩阵进行子集选择来获取选择向量：

nc %>%
  filter(st_intersects(., ash_point, sparse = FALSE)[1,])

需要使用点号.，以便nc也成为st_intersects的参数，而不仅仅是filter的参数。

如果filter.sf方法能够直接敏感地处理st_intersects的输出，而不需要sparse=FALSE和[1,]就好了。我会把它放在某个TODO列表上。

- Edzer Pebesma

---

我明白了，所以我猜想，提供nc的geometry列而不是.作为st_intersects的输入，同样可以将输出更改为向量，就像使用[1,]方法一样。 我同意，对于dplyr动词，有一种预期的工作流程，额外的样板代码来使sf几何验证函数在这个流程中发挥作用并不理想。 - jmw

当我使用[1,]时，得到的结果是不正确的，但是当我使用[,1]时，得到的结果是正确的。也许应该使用[,1]来使用列而不是行？ - Bjørn Kallerud

这个答案对我也不起作用，我得到了“Error in FALSE[1, ] : incorrect number of dimensions”错误。 - Nebulloyd

Answer 2

另一种实现这个的方法是使用st_filter()函数和.predicates = st_intersects参数。

在这种情况下，代码如下：

nc %>% st_filter(ash_point, .predicates = st_intersects)

您还可以将.predicates参数替换为任何类似的st_methods，例如.predicates = st_within，以仅获取完全位于内部的观测值。

- ttalVlatt

Answer 3

显然，为了让dplyr动词与sf函数配合使用，您需要指定列名“geometry”。

更正后的版本：

nc %>% 
  filter(st_intersects(geometry, ash_point, sparse = FALSE))

- jmw

---

另外，请注意st_point_on_surface()会将您的多边形强制转换为点，以便它仅与其边界多边形相交。如果您将其保留为多边形st_polygonize()，它将与其他4个多边形相交。 - Rich Pauloo

是的。我认为这是一个更简单的例子，只使用一个相交多边形的单个点。 - jmw

Answer 4

提醒一下：我使用建议的方法没有成功完成这个用例。但是，直接操作 st_intersects 的输出以创建指示变量对我有用：r - 将 sf::st_within 的输出转换为向量

- Sebastian Kraus

---

欢迎@Sebastian Krays，看起来并不是很有帮助。你可以放置一些代码片段来解决问题。这是我获取更多声誉的建议。 - Willian