CODEWITHSUNDEEP
rgeolocationshapefiler-sf
Allen
Allen

Reputation: 447

Checking if a point falls within polygon Shapefile

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'?

Upvotes: 1

Views: 3355

Answers (2)

Jindra Lacko
Jindra Lacko

Reputation: 8749

I suggest you consider joining your spatial objects via sf::st_join(), as shown bellow; what it does is that it combines the attributes of your polygon objects and points objects.

The default behaviour is "left" join = points lacking polygons will get NA. It can be tweaked by setting left = FALSE in join parameters, resulting in "inner" join behaviour = points not contained in polygons will be omitted from result.

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

Upvotes: 2

sequoia
sequoia

Reputation: 480

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

Upvotes: 4

Related Questions