检查点是否落在多边形 Shapefile 中

Question

我有一个关于纽约市黄色出租车服务区的 shapefile：taxi_zones.shp。可以在这里下载：https://s3.amazonaws.com/nyc-tlc/misc/taxi_zones.zip

我想检查某些位置是否属于任何区域。这是我使用的 R 代码：

library(sf)
tt <- read_sf('taxi_zones.shp')

pnts <- data.frame(
  "x" = c(-73.97817,-74.00668,0,500),
  "y" = c(40.75798, 40.73178,0,400))

pnts_sf <- do.call("st_sfc",c(lapply(1:nrow(pnts), 
                                     function(i) {st_point(as.numeric(pnts[i, ]))}), list("crs" = 4326))) 

pnts_trans <- st_transform(pnts_sf, 2163) 
tt_trans <- st_transform(tt, 2163)   
   
zones <- apply(st_intersects(tt_trans, pnts_trans, sparse = FALSE), 2, 
                     function(col) { 
                       tt_trans[which(col), ]$LocationID 
                     })

前两个点在 shapefile 定义的区域内。然而，第三点不是。第四点的坐标不正确。我应该如何修改代码，以便对于区域外的点和坐标不正确的点，它 returns 'NA'?

Answer 1

我有自己的方法。那能满足你的要求吗？我不能告诉你你的代码具体有什么问题，但是这个代码也更干净一些：

library(sf)
tt <- read_sf('./Downloads/taxi_zones/taxi_zones.shp')


pnts <- data.frame(
  "x" = c(-73.97817, -74.00668, 0, 500),
  "y" = c(40.75798, 40.73178, 0, 400)
)

pnts_sf <- st_as_sf(pnts, coords = c('x', 'y'), crs = st_crs(4326))

pnts_trans <- st_transform(pnts_sf, 2163)
tt_trans <- st_transform(tt, 2163)


pnts_trans <- pnts_sf %>% mutate(
  intersection = as.integer(st_intersects( pnts_trans,tt_trans)))

结果会是

                    geometry intersection
1 POINT (-73.97817 40.75798)          161
2 POINT (-74.00668 40.73178)          158
3                POINT (0 0)           NA
4            POINT (500 400)           NA

Answer 2

我建议您考虑通过 sf::st_join() 连接您的空间对象，如下所示；它的作用是结合多边形对象和点对象的属性。

默认行为是“左”连接 = 缺少多边形的点将得到 NA。可以通过在连接参数中设置 left = FALSE 来调整它，从而导致“内部”连接行为 = 不包含在多边形中的点将从结果中省略。

library(sf)
tt <- read_sf('taxi_zones.shp')

pnts <- data.frame(
  "x" = c(-73.97817,-74.00668,0,500),
  "y" = c(40.75798, 40.73178,0,400))

pnts_sf <- sf::st_as_sf(pnts, coords = c("x", "y"), crs = 4326)


pnts_trans <- st_transform(pnts_sf, 2163) 
tt_trans <- st_transform(tt, 2163)  

res <- sf::st_join(pnts_trans, tt_trans)

print(res)

Simple feature collection with 4 features and 6 fields (with 1 geometry empty)
geometry type:  POINT
dimension:      XY
bbox:           xmin: 2152087 ymin: -130624.1 xmax: 9480615 ymax: 1178046
projected CRS:  NAD27 / US National Atlas Equal Area
  OBJECTID Shape_Leng   Shape_Area                          zone LocationID   borough                  geometry
1      161 0.03580391 7.191307e-05                Midtown Center        161 Manhattan POINT (2153474 -127064.5)
2      158 0.05480999 1.855683e-04 Meatpacking/West Village West        158 Manhattan POINT (2152087 -130624.1)
3       NA         NA           NA                          <NA>         NA      <NA>   POINT (9480615 1178046)
4       NA         NA           NA                          <NA>         NA      <NA>               POINT EMPTY

检查点是否落在多边形 Shapefile 中

Checking if a point falls within polygon Shapefile

r

geolocation

shapefile

sf