sp::over()。该点是否属于使用 OGRGeoJSON 文件标识的多边形之一？

Question

我正在尝试获取一个布尔向量，例如，v[i] =1 告诉我第 i 个点（纬度经度对，存在于 train 数据帧内）是否落在在 OGRGeoJSON 文件标识的地理区域之一内。

OGR file 的结构大致如下：

• 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 是一个较旧的软件包，正在逐步淘汰，取而代之的是较新的“简单功能”sf 软件包。让我知道您是否愿意使用 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"