как объединить шейп-файл с фреймом данных с данными широты / долготы

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 соответствует гранулярности налоговой партии), и подготовить соответствующую карту этого. Мне нужно сделать это, используя sf.

По сути что-то похожее на

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

но откуда рассчитывается каждая налоговая партия из подсчета того, сколько точек в моем кадре данных широты / долготы попадают в них.

Конечно, в моем небольшом примере у меня просто есть 2 точки в одной и той же налоговой партии, так что это всего лишь одна единая налоговая партия во всей области. Мои реальные данные содержат гораздо больше очков.

Я думаю, что есть простой способ сделать это, но я не смог его найти. Спасибо!

Answer 1

Я пробовал это на ваших данных, но пересечение пусто для обоих наборов точек, которые вы указали. Тем не менее, код должен работать.

РЕДАКТИРОВАТЬ: Упрощенный 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 подход от этот ответ .

Answer 2

Так я бы сделал это с произвольными данными многоугольника и точки. Я бы не стал объединять их, а просто использовал бы предикат геометрии, чтобы получить желаемое количество. Здесь мы:

1. Используйте встроенный набор данных nc и преобразуйте его в 3857 crs, который проецируется, а не широтно-длинный (избегает предупреждения в st_contains)
2. Создайте 1000 случайных точек в пределах рамки nc, используя st_bbox и runif. Обратите внимание, что st_as_sf может превратить фрейм data.frame с длинными столбцами в sf точек.
3. Используйте lengths(st_contains(polygons, points), чтобы получить количество точек на многоугольник. sgbp объекты, созданные с помощью предиката геометрии, в основном «для каждой геометрии в sf x, какие индексы геометрии в sf y удовлетворяют предикату». Таким образом, lengths1 фактически дает количество точек, удовлетворяющих предикату для каждой геометрии, в этом случае количество точек, содержащихся в каждом многоугольнике.
4. Как только счетчики в объекте 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()

Создано в 2018-05-02 пакетом представ (v0.2.0).