I am trying to calculate/plot a point on a line, closest to another line. I can calculate the point on that line, both in cartesian coordinates and with great circle distances etc. Where I am having trouble is accurately calculating the distance in miles and plotting it.
I have a few lines defined by lats/lons. I am trying to take the end point of one of the lines, as the point, p3. P1 and p2 are two points along the adjacent line.
When I calculate the cross track distance it doesn't appear to be correct. (image below)
As a side note: I can accurately calculate, in Cartesian coordinates, using arbitrary numbers. Using my code below:
import numpy as np
import matplotlib.pyplot as plt
x1=1; y1=5
x2=6; y2=7
x0=0; y0=0
p1 = x1,y1
p2 = x2,y2
p3 = x0,y0
d = abs((x2-x1)*(y1-y0) - (x1-x0)*(y2-y1)) / np.sqrt(np.square(x2-x1) + np.square(y2-y1))
print('d = '+str(d))
A = y1-y2
B = x2-x1
C = (x1*y2)-(x2*y1)
d2 = abs( A*x0 + B*y0 + C ) / np.sqrt(np.square(B) + np.square(A))
print('d2 = '+str(d2))
# point along the line at distance d
# https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#:~:text=In%20Euclidean%20geometry%2C%20the%20distance,nearest%20point%20on%20the%20line.
x4 = ( B*(B*x0-A*y0) - A*C ) / (A*A + B*B)
y4 = ( A*(A*y0-B*x0) - B*C ) / (A*A + B*B)
P4 = x4;y4
%matplotlib notebook
fig,ax = plt.subplots(figsize=(5,5))
plt.scatter(x1,y1,color="red",marker=".",label="p1")
plt.scatter(x2,y2,color="blue",marker=".",label="p2")
plt.scatter(x0,y0,color="black",marker=".",label="p3")
plt.scatter(x4,y4,color="green",marker=".",label="p4")
# plt.axline((x1, y1), (x2, y2), linestyle='--', color='black', zorder=0)
plt.plot((x0, x4), (y0, y4), linestyle=':', color='black', zorder=0)
plt.axis('equal')
plt.legend()
Here is the plot:
So, I can make the equations work, I'm just having difficulty translating them to/from cartesian/spheical/eliptical and getting accurate distances.
I know I'm plotting these lat/lons as a cartesian plot but it still doesn't add up in my head.
Here is the code I used to calculate/plot:
Found at this website: https://gis.stackexchange.com/questions/209540/projecting-cross-track-distance-on-great-circle
from math import radians, cos, sin, asin, sqrt, degrees, atan2
def validate_point(p):
lat, lon = p
assert -90 <= lat <= 90, "bad latitude"
assert -180 <= lon <= 180, "bad longitude"
# original formula from http://www.movable-type.co.uk/scripts/latlong.html
def distance_haversine(p1, p2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
Haversine
formula:
a = sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ sin²(Δλ/2)
_ ____
c = 2 ⋅ atan2( √a, √(1−a) )
d = R ⋅ c
where φ is latitude, λ is longitude, R is earth’s radius (mean radius = 6,371km);
note that angles need to be in radians to pass to trig functions!
"""
lat1, lon1 = p1
lat2, lon2 = p2
for p in [p1, p2]:
validate_point(p)
R = 6371 # km - earths's radius
# convert decimal degrees to radians
lat1, lon1, lat2, lon2 = map(radians, [lat1, lon1, lat2, lon2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * asin(sqrt(a)) # 2 * atan2(sqrt(a), sqrt(1-a))
d = R * c # convert to km
d = d * 0.621371 # convert km to miles
return d
def spherical2Cart(lat,lon):
clat=(90-lat)*np.pi/180.
lon=lon*np.pi/180.
x=np.cos(lon)*np.sin(clat)
y=np.sin(lon)*np.sin(clat)
z=np.cos(clat)
return np.array([x,y,z])
def cart2Spherical(x,y,z):
r=np.sqrt(x**2+y**2+z**2)
clat=np.arccos(z/r)/np.pi*180
lat=90.-clat
lon=np.arctan2(y,x)/np.pi*180
lon=(lon+360)%360
return np.array([lat,lon,np.ones(lat.shape)])
def greatCircle(lat1,lon1,lat2,lon2,r=None,verbose=True):
'''Compute the great circle distance on a sphere
<lat1>, <lat2>: scalar float or nd-array, latitudes in degree for
location 1 and 2.
<lon1>, <lon2>: scalar float or nd-array, longitudes in degree for
location 1 and 2.
<r>: scalar float, spherical radius.
Return <arc>: great circle distance on sphere.
'''
if r is None:
r=6371 # km
d2r=lambda x:x*np.pi/180
lat1,lon1,lat2,lon2=map(d2r,[lat1,lon1,lat2,lon2])
dlon=abs(lon1-lon2)
numerator=(cos(lat2)*sin(dlon))**2 + \
(cos(lat1)*sin(lat2) - sin(lat1)*cos(lat2)*cos(dlon))**2
numerator=np.sqrt(numerator)
denominator=sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(dlon)
dsigma=np.arctan2(numerator,denominator)
arc=r*dsigma
return arc
def getCrossTrackPoint(lat1,lon1,lat2,lon2,lat3,lon3):
'''Get the closest point on great circle path to the 3rd point
<lat1>, <lon1>: scalar float or nd-array, latitudes and longitudes in
degree, start point of the great circle.
<lat2>, <lon2>: scalar float or nd-array, latitudes and longitudes in
degree, end point of the great circle.
<lat3>, <lon3>: scalar float or nd-array, latitudes and longitudes in
degree, a point away from the great circle.
Return <latp>, <lonp>: latitude and longitude of point P on the great
circle that connects P1, P2, and is closest
to point P3.
'''
x1,y1,z1=spherical2Cart(lat1,lon1)
x2,y2,z2=spherical2Cart(lat2,lon2)
x3,y3,z3=spherical2Cart(lat3,lon3)
D,E,F=np.cross([x1,y1,z1],[x2,y2,z2])
a=E*z3-F*y3
b=F*x3-D*z3
c=D*y3-E*x3
f=c*E-b*F
g=a*F-c*D
h=b*D-a*E
tt=np.sqrt(f**2+g**2+h**2)
xp=f/tt
yp=g/tt
zp=h/tt
result1=cart2Spherical(xp,yp,zp)
result2=cart2Spherical(-xp,-yp,-zp)
d1=greatCircle(result1[0],result1[1],lat3,lon3,r=1)
d2=greatCircle(result2[0],result2[1],lat3,lon3,r=1)
if d1>d2:
return result2[0],result2[1]
else:
return result1[0],result1[1]
def getCrossTrackDistance(lat1,lon1,lat2,lon2,lat3,lon3,r=None):
'''Compute cross-track distance
<lat1>, <lon1>: scalar float or nd-array, latitudes and longitudes in
degree, start point of the great circle.
<lat2>, <lon2>: scalar float or nd-array, latitudes and longitudes in
degree, end point of the great circle.
<lat3>, <lon3>: scalar float or nd-array, latitudes and longitudes in
degree, a point away from the great circle.
Return <dxt>: great cicle distance between point P3 to the closest point
on great circle that connects P1 and P2.
NOTE that the sign of dxt tells which side of the 3rd point
P3 is on.
'''
if r is None:
r=CONS.EARTH_RADIUS
# get angular distance between P1 and P3
delta13=greatCircle(lat1,lon1,lat3,lon3,r=1.)
# bearing between P1, P3
theta13=getBearing(lat1,lon1,lat3,lon3)*np.pi/180
# bearing between P1, P2
theta12=getBearing(lat1,lon1,lat2,lon2)*np.pi/180
dtheta=np.arcsin(sin(delta13)*sin(theta13-theta12))
dxt=r*dtheta
return dxt
def getAlongTrackDistance(lat1,lon1,lat2,lon2,lat3,lon3,r=None):
'''Compute the distance from the start point to closest point to 3rd point
<lat1>, <lon1>: scalar float or nd-array, latitudes and longitudes in
degree, start point of the great circle.
<lat2>, <lon2>: scalar float or nd-array, latitudes and longitudes in
degree, end point of the great circle.
<lat3>, <lon3>: scalar float or nd-array, latitudes and longitudes in
degree, a point away from the great circle.
Return <dat>: distance from the start point to the closest point on
the great circle connecting P1 and P2.
See also getCrossTrackDistance(), getCrossTrackPoint().
'''
if r is None:
r=CONS.EARTH_RADIUS
# angular distance from P1 to P3
delta13=greatCircle(lat1,lon1,lat3,lon3,r=1.)
# angular distance from Pcloset to P3
dxt=getCrossTrackDistance(lat1,lon1,lat2,lon2,lat3,lon3,r=1)
dat=r*np.arccos(cos(delta13)/cos(dxt/r))
return dat
def getBearing(lat1, long1, lat2, long2):
import math
dLon = (long2 - long1)
y = math.sin(dLon) * math.cos(lat2)
x = math.cos(lat1) * math.sin(lat2) - math.sin(lat1) * math.cos(lat2) * math.cos(dLon)
brng = math.atan2(y, x)
brng = np.rad2deg(brng)
return brng
x1 = -75.9671285; y1 = 41.67304238
x2 = -75.96262168; y2 = 41.66425696
x3 = -75.96017288; y3 = 41.67049662
p1=[x1,y1]
p2=[x2,y2]
p3=[x3,y3]
pp=getCrossTrackPoint(p1[0],p1[1],p2[0],p2[1],p3[0],p3[1],)
print('pp',pp)
dxt=getCrossTrackDistance(p1[0],p1[1],p2[0],p2[1],p3[0],p3[1],r=1)
print('dxt',dxt)
dat=getAlongTrackDistance(p1[0],p1[1],p2[0],p2[1],p3[0],p3[1],r=1)
print('dat',dat)
dxt2=greatCircle(pp[0],pp[1],p3[0],p3[1],r=1)
print('dxt2',dxt2)
dat2=greatCircle(pp[0],pp[1],p1[0],p1[1],r=1)
print('dat2',dat2)
%matplotlib notebook
fig,ax = plt.subplots(figsize=(8,8))
plt.plot(p1[0],p1[1],color="blue",marker="o",label="Great Circle Point p1",markersize=12)
plt.plot(p2[0],p2[1],color="green",marker="o",label="Great Circle Point p2",markersize=12)
plt.plot(p3[0],p3[1],color="red",marker="o",label="Point p3",markersize=12)
plt.plot(pp[0],pp[1],color="Black",marker="*",label="Crosstrack Point",markersize=12)
plt.plot([pp[0],p3[0]],[pp[1],p3[1]],linestyle='--',label='Cross track')
print('Cross track dist = ',distance_haversine(pp,p3))
for q in range(df.shape[0]):
toe_lat = df['WGS84_toe_lat'].iloc[q]
toe_lon = df['WGS84_toe_lon'].iloc[q]
heel_lat = df['WGS84_heel_lat'].iloc[q]
heel_lon = df['WGS84_heel_lon'].iloc[q]
plt.plot([toe_lon,heel_lon],[toe_lat,heel_lat],'green')
if q==0:
plt.plot([toe_lon,heel_lon],[toe_lat,heel_lat],'green',label="Well Lateral")
plt.ylabel('Latitude')
plt.xlabel('Longitude')
plt.grid(True)
plt.axis('equal')
plt.legend()
When I calculate the haversine distance from p1 to p3, it calculates 0.48 miles but the GIS software says 0.4 miles. Which is not nearly as accurate as I need.
Does this mean the lines/points I am evaluating are so close that cartesian coordinates will be more accurate? If so, How do I calculate/plot them accurately in miles using the lats/lons?
When I try to convert the lat/lons to radians
With a function simlar to:
def conv_latlon2xy(lat=None,lon=None):
lat, lon = np.deg2rad(lat), np.deg2rad(lon)
# R = 6371 # radius of the earth
x = R * np.cos(lat) * np.cos(lon)
y = R * np.cos(lat) * np.sin(lon)
z = R *np.sin(lat)
return x,y
I do not get accurate results.
In this situation, what is the best way to use lat/lons to calculate the closest distance from a point to a line and be able to accurately measure/plot that distance in miles, with distances less than 1 mile? I can calculate it both in cartesian and spherical/ellpitical but I can't make the jump between the two.
I found the best way to address this problem is using pyproj with code similar to below:
This transform the lat/lons into feet or meters.