Finding a Point On a Line Using Lat/Lon with Distances Less Than a Mile Apart

341 views Asked by At

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)

enter image description here

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:

enter image description here

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.

enter image description here

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.

1

There are 1 answers

0
someds On

I found the best way to address this problem is using pyproj with code similar to below:

transformer = pyproj.Transformer.from_crs(4326, 2271) # WSG84 -> NAD83PAnorth ft 
df['xx_mid'],df['yy_mid'] =   transformer.transform(df['WGS84_mid_lat'].copy(),df['WGS84_mid_lon'].copy())

This transform the lat/lons into feet or meters.