3

I am currently looking for a solution where i need to find if a lat long is on the route or not. I cannot use Google Maps Geometry IsLocationOnEdge() or ContainsLocation() api's as i am working on a server side application.

The approach i am using is to find the perpendicular distance between the point and the line segment and then find the intersection point and check the distance if its within the acceptable limit.

The issue is at getting intersection point on the line segment.

I also tried the answer mentioned here. But still it doesn't work correctly. Here is the Sample data which i am testing.

  • Line Points ([60.3757749,11.2075521],[60.3755244,11.2074115])
  • Point to check (60.3757,11.2077)

I am using the below code snippet to find the intersection point on the line(perpendicular distance) 

def get_perp( X1, Y1, X2, Y2, X3, Y3):
"""************************************************************************************************ 
Purpose - X1,Y1,X2,Y2 = Two points representing the ends of the line segment
          X3,Y3 = The offset point 
'Returns - X4,Y4 = Returns the Point on the line perpendicular to the offset or None if no such
                    point exists
'************************************************************************************************ """
XX = X2 - X1 
YY = Y2 - Y1 
ShortestLength = ((XX * (X3 - X1)) + (YY * (Y3 - Y1))) / ((XX * XX) + (YY * YY)) 
X4 = X1 + XX * ShortestLength 
Y4 = Y1 + YY * ShortestLength
if X4 < X2 and X4 > X1 and Y4 < Y2 and Y4 > Y1:
    return X4,Y4
return None

The value of X4,Y4 i am getting is (60.3758949829255,11.2076194998376), which is not on the line segment.

In the image you can see that (X4,Y4) Is never gonna be perpendicular to it.

enter image description here

So question is what's wrong in the Pseudo code and how to make it correct.

[Updated] as per Kogexo New Projection 

import math
def get_perp( X1, Y1, X2, Y2, X3, Y3):
XX = X2 - X1 
YY = Y2 - Y1 
ShortestLength = ((XX * (X3 - X1)) + (YY * (Y3 - Y1))) / ((XX * XX) + (YY * YY)) 
X4 = X1 + XX * ShortestLength 
Y4 = Y1 + YY * ShortestLength
return X4,Y4
def compute_distance(lat1, lon1, lat2, lon2): 
R = 6378.137 # Radius of earth in KM
dLat = lat2 * math.pi / 180 - lat1 * math.pi / 180
dLon = lon2 * math.pi / 180 - lon1 * math.pi / 180
a = math.sin(dLat/2) * math.sin(dLat/2) + math.cos(lat1 * math.pi / 180) * math.cos(lat2 * math.pi / 180) * math.sin(dLon/2) * math.sin(dLon/2)
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
d = R * c
return d * 1000 # meters

def compute_distance_radians(lat1, lon1, lat2, lon2): 
# R = 6378.137 # Radius of earth in KM
dLat = lat2  - lat1 
dLon = lon2 - lon1
a = math.sin(dLat/2) * math.sin(dLat/2) + math.cos(lat1) * math.cos(lat2 ) * math.sin(dLon/2) * math.sin(dLon/2)
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
# d = R * c
# return d * 1000 # meters
return c

def compute_bearing(lat1, lon1, lat2, lon2):
lat1 = lat1*math.pi/180.0
lat2 = lat2*math.pi/180.0
lon1 = lon1*math.pi/180.0
lon2 = lon2*math.pi/180.0
y = math.sin(lon2-lon1) * math.cos(lat2)
x = math.cos(lat1)*math.sin(lat2) - math.sin(lat1)*math.cos(lat2)*math.cos(lon2-lon1)
brng = math.atan2(y, x)
return brng
X1=60.3757749
Y1=11.2075521
X2=60.3755244
Y2=11.2074115
X3=60.3757
Y3=11.2077
distance = compute_distance(X1, Y1, X2, Y2)
print("distance: "+str(distance))
bearing = compute_bearing(X1, Y1, X2, Y2)
print("bearing: "+str(bearing))
distance3 = compute_distance(X1, Y1, X3, Y3)
print("distance3: "+str(distance3))
bearing3 = compute_bearing(X1, Y1, X3, Y3)
print("bearing3: "+str(bearing3))
Y2p = distance * -1.0 * math.sin(bearing)
X2p = distance * math.cos(bearing)
X3p = distance3 * math.cos(bearing3)
Y3p = distance3 * -1.0 * math.sin(bearing3)
X4,Y4=get_perp( 0, 0, X2p, Y2p, X3p, Y3p)
print('X4,Y4: '+str(X4)+','+str(Y4))
distance4=compute_distance_radians(X2p, Y2p, X4, Y4)
bearing4 = bearing
R = 6378137
X4 = math.asin( math.sin(X1*math.pi/180)*math.cos(distance4/R) +                    math.cos(X1*math.pi/180)*math.sin(distance4/R)*math.cos(bearing4) )
Y4 = Y1 + math.atan2(math.sin(bearing4)*math.sin(distance4/R)*math.cos(X1*math.pi/180),
                     math.cos(distance4/R)-    math.sin(X1*math.pi/180)*math.sin(X4))*180.0/math.pi
X4 = X4*180.0/math.pi
print('X1,Y1: '+str(X1)+','+str(Y1))
print('X2,Y2: '+str(X4)+','+str(Y2))
print('X3,Y3: '+str(X4)+','+str(Y3))
print('X2p,Y2p: '+str(X2p)+','+str(Y2p))
print('X3p,Y3p: '+str(X3p)+','+str(Y3p))
print('X4,Y4: '+str(X4)+','+str(Y4))

Here are the values

 distance: 28.93889806426741
 bearing: -2.870952948395761
 distance3: 11.651304235535246
 bearing3: 2.368297140336325
 X4,Y4: -5.6452886445122985,1.5662681805137304
 X1,Y1: 60.3757749,11.2075521
 X2,Y2: 60.375750021792385,11.2074115
 X3,Y3: 60.375750021792385,11.2077
 X2p,Y2p: -27.885524191583833,7.7367539533484475
 X3p,Y3p: -8.337820728617235,-8.138404996462288
 X4,Y4: 60.375750021792385,11.20753813632277

It never clears the Condition if X4 < X2 and X4 > X1 and Y4 < Y2 and Y4 > Y1: X1 and Y1 is passed 0,0 in the snippet above

Improve this question

1 Answer 1

Reset to default
2

I think you have been fooled by CRS. Indeed, the lat/lng system kind of twists the forms. If you mathematically look at what you are trying to do, you can plot the points you give on a graph with equaled axis.

That gives you the following plot which highlights why the projected point (in yellow) of the reference point (in blue) is not on the segment (between red and green points). Plot highlighting the problem

If you want to fix this, you could maybe try to use another CRS for projection in meters or something more precise for the geographical zone you are using.

Projection in new coordinates systems

Using this website, you can have:

def compute_distance(lat1, lon1, lat2, lon2): 
    R = 6378.137 # Radius of earth in KM
    dLat = lat2 * math.pi / 180 - lat1 * math.pi / 180
    dLon = lon2 * math.pi / 180 - lon1 * math.pi / 180
    a = math.sin(dLat/2) * math.sin(dLat/2) + math.cos(lat1 * math.pi / 180) * math.cos(lat2 * math.pi / 180) * math.sin(dLon/2) * math.sin(dLon/2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
    d = R * c
    return d * 1000 # meters

def compute_bearing(lat1, lon1, lat2, lon2):
    lat1 = lat1*math.pi/180.0
    lat2 = lat2*math.pi/180.0
    lon1 = lon1*math.pi/180.0
    lon2 = lon2*math.pi/180.0
    y = math.sin(lon2-lon1) * math.cos(lat2)
    x = math.cos(lat1)*math.sin(lat2) - math.sin(lat1)*math.cos(lat2)*math.cos(lon2-lon1)
    brng = math.atan2(y, x)
    return brng

Computing new coordinates

That would give you, using point (X1,Y1) as a reference: 

distance = compute_distance(X1, Y1, X2, Y2)
bearing = compute_bearing(X1, Y1, X2, Y2)
distance3 = compute_distance(X1, Y1, X3, Y3)
bearing3 = compute_bearing(X1, Y1, X3, Y3)

Y2p = distance * -1.0 * math.sin(bearing)
X2p = distance * math.cos(bearing)

X3p = distance3 * math.cos(bearing3)
Y3p = distance3 * -1.0 * math.sin(bearing3)

New projection

Then, you can use the previous function get_perp() with the new coordinates. That will give you something like this:

results after projection

Back into latitude/longitude coordinates

Then, you will have to reproject into latitude and longitude system (always with this website):

bearing4 = bearing
R = 6378137
X4 = math.asin( math.sin(X1*math.pi/180)*math.cos(distance4/R) +
                    math.cos(X1*math.pi/180)*math.sin(distance4/R)*math.cos(bearing4) )
Y4 = Y1 + math.atan2(math.sin(bearing4)*math.sin(distance4/R)*math.cos(X1*math.pi/180),
                         math.cos(distance4/R)-math.sin(X1*math.pi/180)*math.sin(X4))*180.0/math.pi
X4 = X4*180.0/math.pi

That will give you:

Final result

Improve this answer
18
  • Yes. I read at lot many places that you might need to convert the lat longs in cartesian coordinates earth centric, but didn't got a concrete article which mentions about how to do conversions. I understand it need to be done through vectors and cross product but not sure what should be the steps for it.
    – Chinmay
    Commented May 28, 2018 at 19:52
  • @Chinmay I updated my answer with a proposition. The function to use is Haversine formula.
    – kogexo
    Commented May 29, 2018 at 8:33
  • Why multiplying it by -1 for Y coordinates Y2p = distance * -1.0 * math.sin(bearing). Will it work in all scenarios of Lat Longs or is it just to make this test case to work?
    – Chinmay
    Commented May 29, 2018 at 10:21
  • Because of the bearing definition (counted positive when counter clockwise) for the -1. See movable-type.co.uk/scripts/latlong.html#ellipsoid for accuracy and compatibility. It should work in all cases but your precision might not be enough everytime.
    – kogexo
    Commented May 29, 2018 at 10:26
  • I tried the above mentioned approach but seems to me that i am doing something wrong. I calculated the distance and bearing for X1Y1,X2Y2 pair and X1Y1,X3Y3 Pair. And then updated my X2Y2 and X3Y3 and passed it to same old function but X1Y1 is (0,0). But his doesn't seems to be working. Can you tell me what wrong i am doing or if you could share some pseudo code in fiddle.
    – Chinmay
    Commented May 29, 2018 at 13:09

Your Answer

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

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