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?

  • How big is your foo1.csv', 10 Gb? Have you seen these articles Geopandas Intersects Speed and GeoPandas performance: optimizing vectorized operations? – Taras Jan 5 '20 at 19:14
  • Hi @Taras! It's too big, is about 1GB. Are you saying that I can use "sjoin" instead of "within"? I think the problem is in foo.shp. The geometry (multipolygon) has many details. It is a map o a biome in Brazil. – Vahid Jan 5 '20 at 19:19
  • 1
    Exactly, I am pointing to spatial-join as was mentioned by @gene. Can you apply that option? – Taras Jan 5 '20 at 19:24
  • @Taras, well suggestion! I will try to use this future. I will inform you here, if I succeed. – Vahid Jan 5 '20 at 19:54
  • 1
    Thank you for endorsing my comment, but I would highly suggest you rewrite your code, expand it with some details/information and put it as a valid answer. Hence it could be much helpful first of all for you, and secondly, let others comprehend explicitly the difference between your initial code and final code. IMHO, you may also extend it with a sort of a processing time calculated for instance with time package. – Taras Jan 6 '20 at 12:27

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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