# Distance of a Lat long from a Polyline

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. 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

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

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). 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: # 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: • 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 May 28 '18 at 19:52
• @Chinmay I updated my answer with a proposition. The function to use is Haversine formula. – kogexo May 29 '18 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 May 29 '18 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 May 29 '18 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 May 29 '18 at 13:09