如何将 shapefile 与具有 latitude/longitude 数据的数据框合并

Question

我正在努力解决以下问题

我已经从这里 https://www1.nyc.gov/site/planning/data-maps/open-data/dwn-pluto-mappluto.page

下载 NYC tax lots PLUTO NYC Manhattan Shapefile

我可以在 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 数据框，这样我就可以计算出每个征税地段中我有多少 latitude/longitude 个观察值（请记住，mydf 是在tax lot粒度），并绘制相应的地图。我需要使用 sf.

本质上类似于

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

但是每个征税地块的计数来自计算我的 latitude/longitude 数据框中有多少点落入其中。

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

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

Answer 1

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

使用内置的 nc 数据集并转换为 3857 crs，它是投影的而不是经纬度的（避免 st_contains 中的警告）
使用 st_bbox 和 runif 在 nc 的边界框内创建 1000 个随机点。请注意，st_as_sf 可以将具有经纬度列的 data.frame 转换为 sf 点。
使用lengths(st_contains(polygons, points) 获取每个多边形的点数。 sgbp 几何谓词创建的对象基本上是 "for each geometry in sf x, what indices of geometries in sf y satisfy the predicate"。因此 lengths1 有效地给出了满足每个几何体谓词的点数，在本例中是每个多边形中包含的点数。
一旦计数在 sf 对象中作为一列，我们就可以 select 并使用 plot.sf 方法绘制它们。

对于您的数据，只需将 nc 替换为 mydf 并省略对 tibble 的调用，而是将 data.frame 与正确的经纬度对一起使用。

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()

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

Answer 2

我在你的数据上试过了，但是你提供的两组点的交集都是空的。但是，代码应该可以工作。

编辑： 简化 group_by + mutate 与 add_count:

mydf = st_read("MN_Dcp_Mappinglot.shp")
xydf = data.frame(lat=c(40.758896,40.758896), lon=c(-73.985130, -73.985130))
xysf = st_as_sf(xydf, coords=c('lon', 'lat'), crs=st_crs(mydf))
## NB: make sure to st_transform both to common CRS, as Calum You suggests
xysf %>% 
    sf::st_intersection(mydf) %>% 
    dplyr::add_count(LOT)

可重现的例子：

nc = sf::st_read(system.file("shape/nc.shp", package="sf"))
ncxy = sf::st_as_sf(data.frame(lon=c(-80, -80.1, -82), lat=c(35.5, 35.5, 35.5)), 
           coords=c('lon', 'lat'), crs=st_crs(nc))
ncxy = ncxy %>% 
           sf::st_intersection(nc) %>%
           dplyr::add_count(FIPS)

## a better approach
ncxy = ncxy %>%
           sf::st_join(nc, join=st_intersects) %>%
           dplyr::add_count(FIPS)

新列 n 包括每个 FIPS 代码的总点数。

ncxy %>% dplyr::group_by(FIPS) %>% dplyr::distinct(n)
> although coordinates are longitude/latitude, st_intersects assumes 
  that they are planar
  # A tibble: 2 x 2
  # Groups:   FIPS [2]
    FIPS     n
   <fctr> <int>
  1  37123     2
  2  37161     1

我不确定为什么您的数据会导致空交集，但由于代码在上面的示例中有效，所以一定有一个单独的问题。

HT：st_join 来自 this answer 的方法。

如何将 shapefile 与具有 latitude/longitude 数据的数据框合并

how to merge a shapefile with a dataframe with latitude/longitude data

r

sf