I have a shapefile about NYC Yellow cab service zones: taxi_zones.shp. It can be download here: https://s3.amazonaws.com/nyc-tlc/misc/taxi_zones.zip
I want to check whether certain locations fall into any of the zones. Here is the R code I use:
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
})
The first two points are within the zones defined by the shapefile. However, the third point is not. And the fourth point has incorrect coordinates. How should I modify the code so that for points outside the zones and points with incorrect coordinates, it returns 'NA'?
I have my own approach. Would that fulfill your requirements? I can't tell you what specifically is wrong with your code, but this one is also a bit cleaner:
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)))
The result would be
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