Accelerating GeoPandas for selecting points inside polygon

Question

I have one dataset with about 10 million points (lat, long). I would like to select the points that fall within a map.

I have this map as a shapely object (.shp). To do this, I transformed the localization of each point into a Point object using points_from_xy. So I transformed my DataFrame into a GeoDataFrame.

My map is in the form of a Multipolygon object. So I used geometry.unary_union to convert it into a unified polygon (I'm not sure if it is correct). Then I used the within method in Geopandas to selects the Points inside the map.

import pandas as pd
import geopandas as gpd
from tqdm import tqdm

map = gpd.read_file('foo.shp')

df = pd.read_csv('foo1.csv')
points = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.longitude, df.latitude))

map = map.geometry.unary_union

within_points = [points.geometry[i].within(map) for i in tqdm(range(points.geometry.count()))]
within_points = points[within_points]

My problem is that this process takes TOO long (about 5 days on my core i5 laptop). I would like to know if there is any way to accelerate it?

Answer 1

With the help of @Taras I could solve my efficiency problem.

As the first try, I thought maybe the tqdm library is causing this delay. So I modified my code to

within_points = [points.geometry[i].within(map) for i in range(points.geometry.count())]

As I wasn't using tqdm to know the remaining time, I waited about one day to see if the process terminates but as it took too long I aborted it. Using tqdm the rate of calculation was about 9.3 it/s and took about 114 hour!

The successful try was to rewrite the code and using geopandas.sjoin method. Before using this method I had to check if both GeoDataFrames have the same coordinate reference system (CRS).

In my the map.crs returned {'init': 'epsg:4674'} but points.crs returned nothing. To fix this, I assigned the result of map.crs to points.crs. After that, I used geopandas.sjoin with `op='within``.

As a comment, there two other options for op in sjoin, 'contains' e 'intersects'.

The modified code was surprisingly rapid and finished in less than 5min.

import pandas as pd
import geopandas as gpd

map = gpd.read_file('foo.shp')

df = pd.read_csv('foo1.csv')
points = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.longitude, df.latitude))
points.crs = {'init': 'epsg:4674'}
within_points = gpd.sjoin(points, map, op = 'within')

I made a time testing with a portion of my dataset. the process with my previous code took about 44.2s while with new code it reduced to 1.05s.