How to use Vectorization with NumPy arrays to calculate geodesic distance using Geopy library for a large dataset?

Question

I am trying to calculate geodesic distance from a dataframe which consists of four columns of latitude and longitude data with around 3 million rows. I used the apply lambda method to do it but it took 18 minutes to finish the task. Is there a way to use Vectorization with NumPy arrays to speed up the calculation? Thank you for answering.

My code using apply and lambda method:

from geopy import distance

df['geo_dist'] = df.apply(lambda x: distance.distance(
                              (x['start_latitude'], x['start_longitude']),
                              (x['end_latitude'], x['end_longitude'])).miles, axis=1)

Updates:

I am trying this code but it gives me the error: ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all(). Appreciate if anyone can help.

df['geo_dist'] = distance.distance(
                          (df['start_latitude'].values, df['start_longitude'].values),
                          (df['end_latitude'].values, df['end_longitude'].values)).miles

Answer 1

I think you might consider using geopandas for this, it's an extension of pandas (and therefore numpy) designed to do these types of calculations very quickly.

Specifically, it has a method for calculating the distance between sets of points in a GeoSeries, which can be a column of a GeoDataFrame. I’m fairly certain that this method leverages numexpr for vectorization.

It should look something like this, where you convert your data frame to a GeoDataFrame with (at least) two GeoSeries columns that you can use for the origin and point destinations. This should return a GeoSeries object:

import pandas as pd
import geopandas as gpd
from shapely.geometry import Point

geometry = [Point(xy) for xy in zip(df.longitude, df.latitude)]
gdf = gpd.GeoDataFrame(df, crs={'init': 'epsg:4326'}, geometry=geometry)

distances = gdf.geometry.distance(gdf.destination_geometry)

Answer 2

The answer to your question: You cannot do what you want to do with geopy. I am not familiar with this package but the error traceback shows that this function and possibly all other functions in this package were not written/designed with vectorized computations in mind.

Now, if you can do with great-circle distances, then I would suggest that you experiment with astropy.coordinates package that my be able to compute separations between points in a vectorial way.

Here is an example based on my answer to a different question: Finding closest point:

from astropy.units import Quantity
from astropy.coordinates import SkyCoord, EarthLocation
from astropy.constants import R_earth
import numpy as np

lon1 = Quantity([-71.312796, -87.645307, -87.640426, -87.635513,
                 -87.630629, -87.625793 ], unit='deg')
lat1 = Quantity([41.49008, 41.894577, 41.894647, 41.894713,
                 41.894768, 41.894830], unit='deg')
lon2 = Quantity([-81.695391, -87.645307 + 0.5, -87.640426, -87.635513 - 0.5,
                 -87.630629 + 1.0, -87.625793 - 1.0], unit='deg')
lat2 = Quantity([41.499498, 41.894577 - 0.5, 41.894647, 41.894713 - 0.5,
                 41.894768 - 1.0, 41.894830 + 1.0], unit='deg')

pts1 = SkyCoord(EarthLocation.from_geodetic(lon1, lat1, height=R_earth).itrs, frame='itrs')
pts2 = SkyCoord(EarthLocation.from_geodetic(lon2, lat2, height=R_earth).itrs, frame='itrs')

Then, distances between the two sets of points can be computed as:

>>> dist = pts2.separation(pts1)
>>> print(dist)
<Angle [ 7.78350849, 0.62435354, 0., 0.62435308, 1.25039805, 1.24353876] deg>

Approximate conversion to distance:

>>> np.deg2rad(pts2.separation(pts1)) * R_earth / u.rad
<Quantity [ 866451.17527216,  69502.31527953,      0.        ,
             69502.26348614, 139192.86680148, 138429.29874024] m>

Compare the first value with what you would get from the geopy's example:

>>> distance.distance((41.49008, -71.312796), (41.499498, -81.695391)).meters
866455.4329098687

EDIT: Actually, quite possibly this may actually give you the geodesic distance that you are after but make sure to check the description of EarthLocation.