2

I am trying to calculate some kind of similarity between points with coordinates in geographical space. I will use an example to make things a bit more clear:

import pandas as pd
import geopandas as gpd
from geopy import distance
from shapely import Point

df = pd.DataFrame({
    'Name':['a','b','c','d'],
    'Value':[1,2,3,4],
    'geometry':[Point(1,0), Point(1,2), Point(1,0), Point(3,3)]
})
gdf = gpd.GeoDataFrame(df, geometry=df.geometry)
print(gdf)
  Name  Value                 geometry
0    a      1  POINT (1.00000 0.00000)
1    b      2  POINT (1.00000 2.00000)
2    c      3  POINT (1.00000 0.00000)
3    d      4  POINT (3.00000 3.00000)

I need a new dataframe containing the distance between each pair of points, their similarity (Manhattan distance in this case) and also their other possible variables (in this case there is only name as additional variable).

My solution is the following:

def calc_values_for_row(row, sourcepoint):  ## sourcepoint is a row of tdf
    sourcename = sourcepoint['Name']
    targetname = row['Name']
    manhattan = abs(sourcepoint['Value']-row['Value'])
    sourcecoord = sourcepoint['geometry']
    targetcoord = row['geometry']
    dist_meters = distance.distance(np.array(sourcecoord.coords), np.array(targetcoord.coords)).meters

    new_row = [sourcename, targetname, manhattan, sourcecoord, targetcoord, dist_meters]
    new_row = pd.Series(new_row)
    new_row.index = ['SourceName','TargetName','Manhattan','SourceCoord','TargetCoord','Distance (m)']
    return new_row

def calc_dist_df(df):
    full_df = pd.DataFrame()
    for i in df.index:
        tdf = df.loc[df.index>i]
        if tdf.empty == False:
            sliced_df = tdf.apply(lambda x: calc_values_for_row(x, df.loc[i]), axis=1)
            full_df = pd.concat([full_df, sliced_df])
    return full_df.reset_index(drop=True)

calc_dist_df(gdf)


### EXPECTED RESULT
  SourceName TargetName  Manhattan  SourceCoord  TargetCoord   Distance (m)
0          a          b          1  POINT (1 0)  POINT (1 2)  222605.296097
1          a          c          2  POINT (1 0)  POINT (1 0)       0.000000
2          a          d          3  POINT (1 0)  POINT (3 3)  400362.335920
3          b          c          1  POINT (1 2)  POINT (1 0)  222605.296097
4          b          d          2  POINT (1 2)  POINT (3 3)  247555.571681
5          c          d          1  POINT (1 0)  POINT (3 3)  400362.335920

It works good as expected, BUT it is extremely slow for big datasets. I am iterating once over each row of the dataframe to slice the gdf once and then I use .apply() on the sliced gdf, but I was wondering if there is a way to avoid the first for loop or maybe another solution to make this algorithm much faster.

NOTE
combination from itertools might not be the solution because the geometry column can contain repeated values
EDIT
This is the distribution of repeated values for the 'geometry' column. As you can see most of the points are repeated and only a few are unique. Distribution of PointID counts

4
  • How many unique POINTs do you actually have in your big datasets? If the cardinality is not very high, I think just lru_cache()ing the distance computation results could help.
    – AKX
    May 29 at 12:25
  • Actually there are many points that are exactly the same, furthermore in this case I'm much more interested in those pairs (source=target) than the distant ones. I never used lru_cache(), but if you think it still might help I can give it a try.
    – Mirko
    May 29 at 12:38
  • 1
    Where distance.distance comes from?
    – Corralien
    May 29 at 12:44
  • Thanks for noticing i forgot to add the import. distance.distance is the best way to calculate distances in WGS84 crs.
    – Mirko
    May 29 at 12:57

3 Answers 3

2

You can use scipy.spatial.distance_matrix. Use .x and .y properties to extract coordinates from shapely Point:

from scipy.spatial import distance_matrix

RADIUS = 6371.009 * 1e3  # meters

cx = gdf.add_prefix('Source').merge(gdf.add_prefix('Target'), how='cross')
coords = np.radians(np.stack([gdf['geometry'].x, gdf['geometry'].y], axis=1))
cx['Distance'] = distance_matrix(coords, coords, p=2).ravel() * RADIUS

r, c = np.triu_indices(len(gdf), k=1)
cx = cx.loc[r * len(df) + c]

Output:

>>> cx
   SourceName  SourceValue           Sourcegeometry TargetName  TargetValue           Targetgeometry       Distance
1           a            1  POINT (1.00000 0.00000)          b            2  POINT (1.00000 2.00000)  222390.167448
2           a            1  POINT (1.00000 0.00000)          c            3  POINT (1.00000 0.00000)       0.000000
3           a            1  POINT (1.00000 0.00000)          d            4  POINT (3.00000 3.00000)  400919.575947
6           b            2  POINT (1.00000 2.00000)          c            3  POINT (1.00000 0.00000)  222390.167448
7           b            2  POINT (1.00000 2.00000)          d            4  POINT (3.00000 3.00000)  248639.765971
11          c            3  POINT (1.00000 0.00000)          d            4  POINT (3.00000 3.00000)  400919.575947
1
  • yes thank you, i was trying to drop the repeated values too but you beat me to it. Thanks a lot it is much faster then my code.
    – Mirko
    May 29 at 15:31
1

How about something like this? If the cardinality of the points is low enough, you can precalculate the distances between all unique pairs (since x.distance(y) == y.distance(x)), and just apply that to the df.

For my RNG example, there's 39,800 rows and 19,900 pairs :)

import math
import random
from itertools import product

import pandas as pd
from shapely import Point

# Generate example dataframe
rng = random.Random(1309510)
N = 200
df = pd.DataFrame({
    'Name': [f'n{x}' for x in range(N)],
    'Value': [math.sin(x) for x in range(N)],
    'geometry': [Point(round(rng.uniform(0, 10), 2), round(rng.uniform(0, 10), 2)) for x in range(N)],
})

# Get all unique point pairs within the dataframe
point_pairs = {frozenset((x, y)): (x, y) for (x, y) in product(df.geometry, df.geometry) if x != y}

# Calculate distances between all point pairs
point_distances = {pair_key: pair[0].distance(pair[1]) for pair_key, pair in point_pairs.items()}

# Generate dataframe with all point pairs and their associated data
df = df.merge(df, how='cross', suffixes=('_1', '_2'))
df = df[df.Name_1 != df.Name_2]

# Read distances from precalculated dictionary
df["distance"] = df.apply(lambda x: point_distances[frozenset((x.geometry_1, x.geometry_2))], axis=1)

print(df)
-1

In pandas you have the method diff() which calculates the difference between one value and it's preivous index in the same column. In a case like this you need to create a row with a continuous value and apply the .diff(). Remember that you need to insert a NaN value as the first index since the diff() method will create a list of values with len-1. This is:

import pandas as pd


df = pd.DataFrame({
    'Name':['a','b','c','d'],
    'Value':[1,2,3,4],
    'geometry':[(1,0), (1,2), (1,0), (3,3)]
})

df['first_val'] = df.geometry.str[0]
df['second_val'] = df.geometry.str[1]

df['first_diff'] = df.first_val.diff()
df['second_diff'] = df.second_val.diff()

row_list = []
for idx, rows in df.iterrows():
    my_list = [rows.first_diff, rows.second_diff]
    row_list.append(my_list)

df['geometrical_distance'] = row_list
print(df)

       Name Value geometry   first_val second_val first_diff geometrical_distance
 0      a     1    (1, 0)       1           0       NaN NaN     [nan, nan]
 1      b     2    (1, 2)       1           2       0.0 2.0     [0.0, 2.0]
 2      c     3    (1, 0)       1           0       0.0 -2.0    [0.0, -2.0]
 3      d     4    (3, 3)       3           3       2.0 3.0     [2.0, 3.0]

In case you are tyring to calculate the distance between geographical points you can use haversine. It has a method to calculate the distance between geographical points passing their geometrical coordinate. This method is:

import haversine
from haversine import Unit

loc1=(35.526954, 44.659832)
loc2=(36.215489, 45.625896)
haversine.haversine(loc1, loc2, unit=Unit.METERS)

distances = []
for row_index in range(len(df)):
distances.append(
    haversine.haversine(
        df['geometrical_point_1'].iloc[row_index],
        df['geometrical_point_2'].iloc[row_index], unit=Unit.METERS)
    )
5
  • Sorry this is really not what i was looking for. The result does not even contain the right number of rows. In general for a set of n elements there are n(n-1)/2 pairs, so in this case 6 elements. Plus in this way you cannot apply a function but only calculate differences between rows (not what i asked). Please edit your answer so that is related to the question.
    – Mirko
    May 29 at 12:07
  • Sorry if it was not the solution you were looking for. Seems like it's not very clear to me what you are trying to acomplish in this case. If it's just the difference between geometrical points you can use harversine as I said in the edit. If it's none of both solutions then I'll delete my answer since it's incorrect. May 29 at 12:11
  • You are calculating distances only between consecutive rows, I need distances between all possible rows. The problem is not how to calculate the distance (in my code you can see the line dist_meters = distance.distance(np.array(sourcecoord.coords), np.array(targetcoord.coords)).meters is doing exactly that). The problem is that for each source point with index i there are n-i possible distances, and this makes the code very slow for large n.
    – Mirko
    May 29 at 12:16
  • Check itertools.combinations where you can extract a set of all possible combinations in a list. So, from there, you can extract the distances. May 29 at 12:20
  • @JoseM.González OP already mentions combinations may not be suitable in the original post...
    – AKX
    May 29 at 12:24

Your Answer

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

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