r - 如何将 shapefile 与具有纬度/经度数据的数据框合并

Question

我正在努力解决以下问题

我已经下载了PLUTO NYC Manhattan Shapefile对于纽约市税务地段 从这里 https://www1.nyc.gov/site/planning/data-maps/open-data/dwn-pluto-mappluto.page

我可以在 sf 中阅读它们用一个简单的st_read

> mydf
Simple feature collection with 42638 features and 90 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 971045.3 ymin: 188447.4 xmax: 1010027 ymax: 259571.5
epsg (SRID):    NA
proj4string:    +proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs
First 10 features:
   Borough Block  Lot  CD CT2010 CB2010 SchoolDist Council ZipCode FireComp PolicePrct HealthCent HealthArea
1       MN  1545   52 108    138   4000         02       5   10028     E022         19         13       3700

我的问题如下:我有一个数据框如下

> data_frame('lat' = c(40.785091,40.785091), 'lon' = c(-73.968285, -73.968285))
# A tibble: 2 x 2
        lat        lon
      <dbl>      <dbl>
1 40.785091 -73.968285
2 40.785091 -73.968285

我想将此数据合并到 mydf上面的数据框，这样我就可以计算每个税收批处理内有多少纬度/经度观测值(请记住， mydf 处于税收批处理粒度)，并绘制相应的 map 。我需要使用 sf .

本质上类似于

pol <- mydf %>% select(SchoolDist)
plot(pol)

但是每个税收批处理的计数来自计算我的纬度/经度数据框中有多少点落入其中。

当然，在我的小例子中，我在同一个税区只有 2 个点，所以这只会突出整个区域中的一个税区。我的真实数据包含更多的点。

我认为有一种简单的方法可以做到这一点，但我找不到它。
谢谢!

Answer 1

这就是我对任意多边形和点数据的处理方式。我不会将两者合并，而只是使用几何谓词来获取您想要的计数。在这里，我们:

使用内置 nc数据集并转换为 3857 crs，它是投影而不是经纬度(避免 st_contains 中的警告)

在 nc 的边界框内创建 1000 个随机点, 使用 st_bbox和 runif .请注意 st_as_sf可以将具有 lat long 列的 data.frame 转换为 sf点。

使用lengths(st_contains(polygons, points)获取每个多边形的点数。 sgbp几何谓词创建的对象基本上是“对于 sf x 中的每个几何，sf y 中的几何索引满足谓词”。所以lengths1有效地给出满足每个几何的谓词的点数，在这种情况下，每个多边形中包含的点数.

一旦计数在 sf对象作为列，我们可以 select并用 plot.sf 绘制它们方法。

对于您的数据，只需替换 nc与 mydf并忽略对 tibble 的调用，而是使用您的 data.frame与正确的 lat long 对。

library(tidyverse)
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3
nc <- system.file("shape/nc.shp", package="sf") %>%
  read_sf() %>%
  st_transform(3857)
set.seed(1000)
points <- tibble(
  x = runif(1000, min = st_bbox(nc)[1], max = st_bbox(nc)[3]),
  y = runif(1000, min = st_bbox(nc)[2], max = st_bbox(nc)[4])
) %>%
  st_as_sf(coords = c("x", "y"), crs = 3857)

plot(nc$geometry)
plot(points$geometry, add = TRUE)

nc %>%
  mutate(pt_count = lengths(st_contains(nc, points))) %>%
  select(pt_count) %>%
  plot()

创建于 2018-05-02 由 reprex package (v0.2.0)。