зр :: над (). Принадлежит ли точка одному из полигонов, обозначенных файлом OGRGeo JSON?

Question

Я пытаюсь получить булевский вектор, где, например, v [i] = 1 сообщает мне, попадает ли i-я точка (пара широты и долготы, присутствующая в train dataframe) в одна из географических областей, указанных в файле OGRGeo JSON.

Файл OGR имеет примерно такую ​​структуру:

• Район 1: 24 многоугольных
• Район 2: 4 многоугольных
• Район 3: 27 многоугольных
• Район 4: 18 многоугольных
• Район 5: 34 многоугольных

Это то, что я пытался сделать.

Однако полученные результаты не верны , потому что сгенерированный многоугольник представляет собой смесь всех различных областей, присутствующих в OGR файл.

library(rgdal)
library(httr)
library(sp)

r <- GET('https://data.cityofnewyork.us/api/geospatial/tqmj-j8zm?method=export&format=GeoJSON')
nyc_neighborhoods <- readOGR(content(r,'text'), 'OGRGeoJSON', verbose = F)

#New York City polygonal
pol_lat <- c(nyc_neighborhoods_df$lat)
pol_long <- c(nyc_neighborhoods_df$long)
xy <- cbind(pol_lat, pol_long)
p = Polygon(xy)
ps = Polygons(list(p),1)
pol = SpatialPolygons(list(ps))

#Points to analyse (pair of coordinates)
ny_lat <- c(train$pickup_latitude, train$dropoff_latitude)
ny_long <- c(train$pickup_longitude, train$dropoff_longitude)
ny_coord <- cbind(ny_lat, ny_long)
pts <- SpatialPoints(ny_coord)

#Query: Does the point to analyze fall in or out NYC?
over(pts, pol, returnList = TRUE)

Как это исправить, чтобы получить правильный результат?

Answer 1

sp - это более старый пакет, который постепенно сокращается в пользу более нового пакета "Simple Features" sf. Дайте мне знать, если вы открыты для использования оператора pipe %>% из пакета magrittr, поскольку он прекрасно работает с пакетом sf (как и dplyr и purrr).

Используя sf, вы могли бы сделать:

library(sf)

# Replace this with the path to the geojson file
geojson_path <- "path/to/file.geojson"

boroughs <- sf::st_read(dsn = geojson_path, stringsAsFactors = FALSE)

Теперь создаем очень простой пространственный точечный объект для замены данных «поездов».

# Make test data.frame

test_df <- 
  data.frame(
  # Random test point I chose, a couple of blocks from Central Park
      a = "manhattan_point", 
      y = 40.771959, 
      x = -73.964128, 
      stringsAsFactors = FALSE)

# Turn the test_df into a spatial object
test_point <-
  sf::st_as_sf(
    test_df,
    # The coords argument tells the st_as_sf function
    # what columns store the longitude and latitude data
    # which it uses to associate a spatial point to each
    # row in the data.frame
    coords = c("x", "y"), 
    crs = 4326 # WGS84
    )

Теперь мы готовы к определить, в какой полигон (ы) попадает наша точка:

# Get the sparse binary predicate. This will give a list with as 
# many elements as there are spatial objects in the first argument, 
# in this case, test_point, which has 1 element.
# It also has attributes which detail what the relationship is
# (intersection, in our case) 
sparse_bin_pred <- sf::st_intersects(test_point, boroughs)

# Output the boro_name that matched. I think the package purrr
# offers some more intuitive ways to do this, but
lapply(
  sparse_bin_pred, 
  function(x) boroughs$boro_name[x]
  )

Эта последняя часть выводит:

[[1]]
[1] "Manhattan"