sfc_multipoint对象中的点数

Question

我想根据sfc_multipoint的大小过滤sfc_multipoint数据，即每一行中包含的点/坐标对的数目。假设我希望只保留以下数据集中不超过3个点的行(基于mapview::breweries )

library(mapview)
library(sf)
library(tidyverse)

pts <- breweries %>% 
  group_by(zipcode) %>% 
  summarise()

为了做到这一点，我需要知道每个sfc_multipoint中的点数。理论上，我可以通过as("spatial")和计数行导出坐标，但这不仅很难看，而且对于任何实际大小的数据集也是不可行的。有什么建议吗？

Answer 1

您可以尝试dplyr::count()或summarize(n = n())来读取给定zipcode中的行数，但是breweries数据集似乎有一些重复项，因此这可能会产生误导。

breweries %>% 
  count(zipcode)

#----------
  zipcode n                       geometry
1   90562 1      POINT (11.15795 49.53495)
2   90614 1       POINT (10.85194 49.4208)
3   90763 1      POINT (10.99625 49.44171)
4   91054 3 MULTIPOINT ((11.00901 49.59...
5   91097 2 MULTIPOINT ((10.76099 49.59...

或者只有唯一的点(注意91054的变化)

breweries %>%
  distinct(zipcode, geometry) %>%
  count(zipcode)

#-----
  zipcode n                       geometry
1   90562 1      POINT (11.15795 49.53495)
2   90614 1       POINT (10.85194 49.4208)
3   90763 1      POINT (10.99625 49.44171)
4   91054 2 MULTIPOINT ((11.00901 49.59...
5   91097 2 MULTIPOINT ((10.76099 49.59...

你也可以试试mapview::npts()

breweries %>%
  group_by(zipcode) %>%
  summarize() %>%
  rowwise() %>%
  mutate(n = npts(geometry))

#----
   zipcode                                              geometry     n
 * <chr>                                          <GEOMETRY [°]> <dbl>
 1 90562                               POINT (11.15795 49.53495)     1
 2 90614                                POINT (10.85194 49.4208)     1
 3 90763                               POINT (10.99625 49.44171)     1
 4 91054   MULTIPOINT ((11.00901 49.59511), (11.00505 49.60255))     2
 5 91097   MULTIPOINT ((10.76099 49.59044), (10.76954 49.59015))     2

Answer 2

您可以在sf::st_coordinates()基础上构建--对于多点，产生1)坐标矩阵(duh.)有一个有趣的副作用。( 2) L1列，识别坐标所指的原始多点特征。

因此，考虑这段代码，建立在原始pts示例的基础上(从原始几何学集合到多点转换)。

library(mapview)
library(sf)
library(tidyverse)

pts <- breweries %>% 
  group_by(zipcode) %>% 
  summarise() %>% 
  st_cast("MULTIPOINT")

# get count of coordinates per zipcode
pts$breweries <- st_coordinates(pts) %>% 
  as.data.frame() %>% 
  group_by(L1) %>% # id of the row from original multipoint object
  tally() %>% 
  pull(n)

# zipcodes with less than 3 breweries
pts %>% 
  filter(breweries < 3)

#Simple feature collection with 65 features and 2 fields
#Geometry type: MULTIPOINT
#Dimension:     XY
#Bounding box:  xmin: 9.462785 ymin: 48.90074 xmax: 11.93539 ymax: 50.44162
#Geodetic CRS:  WGS 84
# A tibble: 65 × 3
#   zipcode                                   geometry breweries
# * <chr>                             <MULTIPOINT [°]>     <int>
# 1 90562                        ((11.15795 49.53495))         1
# 2 90614                         ((10.85194 49.4208))         1
# 3 90763                        ((10.99625 49.44171))         1
# 4 91054   ((11.00901 49.59511), (11.00505 49.60255))         2
# 5 91097   ((10.76099 49.59044), (10.76954 49.59015))         2
# 6 91126                         ((10.9281 49.27472))         1
# 7 91207                        ((11.22997 49.55471))         1
# 8 91217                        ((11.42834 49.50683))         1
# 9 91220                         ((11.36851 49.5618))         1
#10 91227                        ((11.30872 49.45008))         1
## … with 55 more rows