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])
    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)

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

0    False
1    False
2    False
3     True
4    False
dtype: bool

With a spatial join

             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:
0    Brooklyn
if pointInPoly.empty:
  • Thanks Gene! This speeded things up considerably using spatial joins. Now, I can run 15 million records in roughly 1 hour. Commented Sep 25, 2016 at 7:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.