Finding if point exists in shapefile using shapely performance?

Question

I have a data frame that consists of 2 million rows (lat, lon) and I am trying to figure out if each of these points are in a boundary defined by a shapefile or not (intact I am looking at multiple shapefiles). Here is the code, I have used to run this, but it is taking me around 30+ hours to run it. I am positive that there might be a better/faster way to do this and I am doing something incorrect. I have also tried the approaches suggested here but to no improvement in performance (R-Tree) R-Tree Spatial Index, Aggregating Points

Here is the code that I am currently using: To upload the shapefiles, I'm using geopandas, shapley, and pandas. I have downloaded the shapefiles from: www.zillow.com/howto/api/neighborhood-boundaries.htm, geo.nyu.edu/catalog/nyu_2451_34154, and geo.nyu.edu/catalog/nyu_2451_34496

import fiona as fn
import shapely as sh
import geopandas as gpd
from shapely.geometry import Point
import pandas as pd

nyc_boroughs_df = gpd.read_file('./data/shapefiles/NYshapefiles/nyc_boroughs/nybb-geo_json.json')
nyc_neighborhood_df = gpd.read_file('./data/shapefiles/NYshapefiles/ZillowNeighborhoods-NY/ZillowNeighborhoods-NY.shp')
nj_neighborhood_df = gpd.read_file('./data/shapefiles/NJshapefiles/ZillowNeighborhoods-NJ/ZillowNeighborhoods-NJ.shp')
ct_neighborhood_df = gpd.read_file('./data/shapefiles/CTshapefiles/ZillowNeighborhoods-CT/ZillowNeighborhoods-CT.shp')
nyc_counties_df = gpd.read_file('./data/shapefiles/NYshapefiles/nyc_metro_counties/nyu_2451_34496-geojson.json')

def get_boro_city_county_name(row):
lat = row.lat
lon = row.lon
bn_series = nyc_boroughs_df[nyc_boroughs_df.geometry.contains(Point(lon, lat))].BoroName
if len(bn_series) == 0:
    bn_series = nyc_neighborhood_df[nyc_neighborhood_df.geometry.contains(Point(lon, lat))].CITY
    if len(bn_series) == 0:
        bn_series = nj_neighborhood_df[nj_neighborhood_df.geometry.contains(Point(lon, lat))].CITY
        if len(bn_series) == 0:
            bn_series = ct_neighborhood_df[ct_neighborhood_df.geometry.contains(Point(lon, lat))].CITY
            if len(bn_series) == 0:
                bn_series = nyc_counties_df[nyc_counties_df.geometry.contains(Point(lon, lat))].name

if len(bn_series) > 0:
    return str(bn_series.iat[0])
else: 
    return ''

df is a dataframe of 2 mil rows (columns: lat, lon) and when I apply the function on each of the rows, it takes around 30 hours to run.

df['name'] = df.apply(lambda x: get_boro_city_county_name, axis=1)

Answer 1

They are many problems in your script

1) why import fiona, shapely and pandas if you do not use them. In addition, GeoPandas use Fiona to import and save the shapefiles, Shapely for the geometry and Pandas for the datas (the module import them)

You need only:

import geopandas as gpd
from shapely.geometry import Point

2) you need to use here spatial joins (spatial_joins.ipynb)

Your solution

nyc_boroughs_df.geometry.contains(Point(1000482.953749179,190236.2289404646))
0    False
1    False
2    False
3     True
4    False
dtype: bool

With a spatial join

point2.head()
             geometry                                id
0   POINT (1000482.953749179 190236.2289404646)     None
pointInPoly = gpd.sjoin(point2, nyc_boroughs_df, how='left',op='within')
print pointinPoly
               geometry                       id    index_right BOROCODE BORONAME  SHAPE_AREA    SHAPE_LEN    SHAPE_LENG
0  POINT (1000482.953749179 190236.2289404646) None     3         3       Brooklyn 1.991519e+09 592422.077454 592422.078518

if not pointInPoly.empty:
     print(pointInPoly.BORONAME)
0    Brooklyn
if pointInPoly.empty:
   ...