1

Given a vector of two points with direction from point 1 to point 2:

latitude_1 longitud_1 latitude_2 longitud_2
41.378118   2.169557   41.389474    2.156421

Find points in the dataset that are in the same "direction":

Latitude Longitud 
41.383093   2.181116
41.373258   2.159358
41.385252   2.168779
41.390692   2.148911

In the image below, the 3 red points in front of the vector would be in the same direction. For direction we can divide 360 degrees in 6 equal parts. A 4 part split is also good.

Expected output would be:

latitude_1 longitud_1 latitude_2 longitud_2   Number of points
41.378118   2.169557   41.389474    2.156421  1

Things I have tried with a 4 quadrant split:

def relative_location(row):
    if row["lat_dif"] > 0 and row["long_dif"] > 0:
        return "North East"
    elif row["lat_dif"] < 0 and row["long_dif"] > 0:
        return "South East"
    elif row["lat_dif"] > 0 and row["long_dif"] < 0:
        return "North West"
    elif row["lat_dif"] < 0 and row["long_dif"] < 0:
        return "South West"
    elif row["lat_dif"] == 0 and row["long_dif"] < 0:
        return "West"
    elif row["lat_dif"] == 0 and row["long_dif"] > 0:
        return "East"
    elif row["lat_dif"] > 0 and row["long_dif"] == 0:
        return "North"
    elif row["lat_dif"] < 0 and row["long_dif"] == 0:
        return "South"
    else: return "Same point"

df["relative_location"] = df.apply(relative_location, axis=1)

enter image description here

2
  • In your graphic, the 0 or 360° position is not correct (mathematical and not geographic) – gene Nov 9 '20 at 15:37
  • You can use a clustering algorithm like kmeans or DBscan (I thinks is more adapt for this case) for regroup the points according their angle with the reference point. – Atm Nov 9 '20 at 15:46
5

Try using pyproj:

import geopandas as gpd
import pyproj
import numpy as np

geodesic = pyproj.Geod(ellps='WGS84') #See: https://stackoverflow.com/questions/54873868/python-calculate-bearing-between-two-lat-long

df = gpd.read_file('/home/bera/Desktop/tempgis/random_points.shp')

#Find azimuth of the two points by using their indexes
p1 = df.iloc[2595]
p2 = df.iloc[3110]
fwd_azimuth_goal,back_azimuth,distance = geodesic.inv(p1.geometry.x, p1.geometry.y, 
                                                 p2.geometry.x, p2.geometry.y)
#fwd_azimuth_goal
#Out[36]: 50.920292718041104

#Find azimuth from the first point to all other points and add as column
az = []
for ix in df.index.tolist():
    fwd_azimuth,back_azimuth,distance = geodesic.inv(df.iloc[2595].geometry.x, df.iloc[2595].geometry.y,
                                                     df.iloc[ix].geometry.x, df.iloc[ix].geometry.y)
    az.append(fwd_azimuth)
df['az'] = az

df['match'] = np.where((df['az']<fwd_azimuth_goal+5)&(df['az']>fwd_azimuth_goal-5),1,0)
df.to_file('/home/bera/Desktop/tempgis/random_points_result.shp')

Filter match=1: enter image description here

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.