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

  • Four part split shoud be relatively simple comparison of x and y coordinates. Given point 1, check if x and y coords of point 2 are smaller or larger. Any other point fulfilling the same condition is in the same quadrant. – martinfleis Nov 9 at 13:55
  • yeah. The case for 4 quadrants is easy. I was hoping I could get help on the 6 directions case. – italo Nov 9 at 14:00
  • Sorry @BERA the graph doesn't match the point given. The answer would be different. – italo Nov 9 at 14:34
  • @BERA I updated the number of points value to match the expected result. :) Thanks for the link – italo Nov 9 at 15:26
  • In your graphic, the 0 or 360° position is not correct (mathematical and not geographic) – gene Nov 9 at 15:37
4

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')

#You could bin using pandas cut. But I'm to tired right now to get it to work: 
#bins = list(range(-405, 450, 90))
#labels = ['N','E','S','W','N','E','S','W','N']
#df['compass']=pd.cut(df['az'], bins=bins, labels=labels)

Filter match=1: enter image description here

| improve this answer | |

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.