How to check how many points are in a specific polygon?

Viewed 58

I am quite new to R and working on a project where I want to draw a Heatmap of New York plotting Pickup Locations from Taxi/Uber rides in different Taxi Zones. I have a dataset of 98k rows of Coordinates defining 263 different Taxi Zones of NYC. In addition, I have a df containing all Pickup Locations from July, 2010 which are around 14.5 million locations. I want to check which Location lies in which Taxizone and sum the number of trips started in each Taxizone to plot it into a heat map. This is the code I am currently working on but the Column "TripsCountPerId" doesn't increase and I think the loop will take really long for 15 million rows.The point.in.polygon function is from the sp package. There are probably some other functions doing similar things and I am open to use these.

for (row in 1:nrow(tlc1007)) {    
  point.x = tlc1007[row,"latPu"]
  point.y = tlc1007[row, "longPu"]
  for (i in 1:nrow(ZoneList)){
    if (point.in.polygon(point.x,point.y,as.numeric(unlist(ZoneList[i,"lat"][[1]])),as.numeric(unlist(ZoneList[i,"long"][[1]])))){
      ZoneList$TripsCountPerId[i]  <- ZoneList$TripsCountPerId[i] +1
    }
    if (i%%100000==0){
      print(paste(i," rows checked"))
    }
  }
} 

tlc1007 is the data frame containing the trips in NYC in July,2010 and ZoneList contains coordinates for each Taxizone. I created the TripsCountPerId column to have a sum of trips for each Taxizone because from this point I can plot the map easier.

For your understanding I tried to build a simple example which differs a little:

library(sp)
library(dplyr)

lat <- c(0,2)
long <- c(0,2)
square <- as.data.frame(cbind(lat,long))
point1 <- c(1,1)
point2 <- c(1,3)
points <- as.data.frame(rbind(point1,point2)) %>% 
  rename("lat"="V1","long"="V2")


for (row in 1:2) {     
  point.x = points[row,"lat"]
  point.y = points[row, "long"]
  for (i in 1:2){
    if (point.in.polygon(point.x,point.y, square$lat[i],square$long[i])){
      print("TRUE") 
    } else {
      print("FALSE")
  }
} 
1 Answers

I think all you need is this:

ZoneList <- list()
for (i in 1:2){
  ZoneList[[i]] <- sum(point.in.polygon(points$lat, points$long, square$lat[i],square$long[i]))
} 

point.in.polygon takes can take vectors of points, and this will likely be much faster than looping over the points yourself.

Related