等价于 `poly.counts` 来计算 lat/long 对落在带有 sf 包的多边形内

Question

sf 包提供了一种很好的处理地理特征的方法，但我想不出一个简单的等价于 GISTools 包中的 poly.counts 函数的方法，它需要 sp 个对象。

poly.counts 计算 SpatialPointsDataFrame 落在 SpatialPolygonsDataFrame 的多边形内的点数，可以按如下方式使用：

数据

## Libraries
library("GISTools")
library("tidyverse")
library("sf")
library("sp")
library("rgdal")
## Obtain shapefiles
download.file(url = "https://www2.census.gov/geo/tiger/TIGER2016/STATE/tl_2016_us_state.zip", destfile = "data-raw/states.zip")
unzip(zipfile = "data-raw/states.zip", exdir = "data-raw/states")
sf_us_states <- read_sf("data-raw/states")

## Our observations:
observations_tibble <- tribble(
  ~lat, ~long,
31.968599,  -99.901813,
35.263266,  -80.854385,
35.149534,  -90.04898,
41.897547,  -84.037166,
34.596759,  -86.965563,
42.652579,  -73.756232,
43.670406,  -93.575858
)

计算每个多边形的点数

我生成了我的 sp 个对象：

sp_us_states <- as(sf_us_states, "Spatial")

observations_spdf <- observations_tibble %>%
  select(long, lat) %>% # SPDF want long, lat pairs
  SpatialPointsDataFrame(coords = .,
                         data = .,
                         proj4string = sp_us_states@proj4string)

现在我可以使用了poly.counts

points_in_states <-
  poly.counts(pts = observations_spdf, polys = sp_us_states)

将其添加到 sp 对象中：

sp_us_states$points.in.state <- points_in_states

现在我已经完成转换回 sf 对象并且可以可视化如下：

library("leaflet")
updated_sf <- st_as_sf(sp_us_states)
updated_sf %>%
  filter(points.in.state > 0) %>%
  leaflet() %>%
  addPolygons() %>%
  addCircleMarkers(
    data = observations_tibble
  )

问题

我可以执行此操作而无需在 sf 和 sp 对象之间进行繁琐的转换吗？

Answer 1

尝试以下操作：

sf_obs = st_as_sf(observations_tibble, coords = c("long", "lat"), 
    crs = st_crs(sf_us_states))
lengths(st_covers(sf_us_states, sf_obs))
# check:
summary(points_in_states - lengths(st_covers(sf_us_states, sf_obs)))

st_covers returns 一个列表，其中包含每个州覆盖的点的索引； lengths returns 这些向量的长度向量，或点数。您将看到的警告表明，尽管您有地理坐标，但底层软件假定它们是笛卡尔坐标（对于这种情况，这很可能没有问题；如果您想以正确的方式摆脱它，请移至投影坐标)

等价于 `poly.counts` 来计算 lat/long 对落在带有 sf 包的多边形内

Equivalent of `poly.counts` to count lat/long pairs falling inside of polygons with the sf package

gis

r

sp

sf

数据

计算每个多边形的点数

问题