Android Question Get shortest distance from point to line on map

[RB Smissaert]

Have a point and a line on a map and need the shortest distance (in meters) between that point and that line.
The point is set by a latitude, longitude type and the lines is set by 2 of these types.

This seems the simplest way to do this:

For distance up To a few thousands meters I would simplify the issue from sphere To plane. Then, the issue Is pretty simply As a easy triangle calculation can be used:
We have points A And B And look For a distance X To line AB. Then:

Location a;
Location b;
Location x;

double ax = a.distanceTo(x)
double alfa = ((Math.abs(a.bearingTo(b) - a.bearingTo(x))) / 180) * Math.PI
double distance = Math.sin(alfa) * ax

Having some trouble though to translate this to B4A and sofar been unable to get the right result.

For the distance I have this Sub:

Public Sub DistVincenty(tLL1 As TMapLatLng, tLL2 As TMapLatLng) As ResumableSub
    'INPUTS: Latitude and Longitude of initial and destination points in decimal format.
    'OUTPUT: Distance between the two points in Meters.
    '
    '=================================================================================
    ' Calculate geodesic distance (in m) between two points specified by latitude/longitude (in numeric [decimal] degrees)
    ' using Vincenty inverse formula for ellipsoids
    '=================================================================================
    ' Code has been ported by lost_species from www.aliencoffee.co.uk to VBA from javascript published at:
    ' http://www.movable-type.co.uk/scripts/latlong-vincenty.html
    ' * from: Vincenty inverse formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the
    ' *       Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975
    ' *       http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
    'Additional Reference: http://en.wikipedia.org/wiki/Vincenty%27s_formulae
    '=================================================================================
    ' Copyright lost_species 2008 LGPL http://www.fsf.org/licensing/licenses/lgpl.html
    '=================================================================================
    ' Code modifications to prevent "Formula Too Complex" errors in Excel (2010) VBA implementation
    ' provided by Jerry Latham, Microsoft MVP Excel Group, 2005-2011
    ' July 23 2011
    '=================================================================================

    Dim low_a As Double
    Dim low_b As Double
    Dim f As Double
    Dim L As Double
    Dim U1 As Double
    Dim U2 As Double
    Dim sinU1 As Double
    Dim sinU2 As Double
    Dim cosU1 As Double
    Dim cosU2 As Double
    Dim lambda As Double
    Dim lambdaP As Double
    Dim iterLimit As Int
    Dim sinLambda As Double
    Dim cosLambda As Double
    Dim sinSigma As Double
    Dim cosSigma As Double
    Dim sigma As Double
    Dim sinAlpha As Double
    Dim cosSqAlpha As Double
    Dim cos2SigmaM As Double
    Dim C As Double
    Dim uSq As Double
    Dim upper_A As Double
    Dim upper_B As Double
    Dim deltaSigma As Double
    Dim s As Double    ' final result, will be returned rounded to 3 decimals (mm).
    'added by JLatham to break up "too Complex" formulas
    'into pieces to properly calculate those formulas as noted below
    'and to prevent overflow errors when using Excel 2010 x64 on Windows 7 x64 systems
    Dim P1 As Double    ' used to calculate a portion of a complex formula
    Dim P2 As Double    ' used to calculate a portion of a complex formula
    Dim P3 As Double    ' used to calculate a portion of a complex formula

    'See http://en.wikipedia.org/wiki/World_Geodetic_System
    'for information on various Ellipsoid parameters for other standards.
    'low_a and low_b in meters
    ' === GRS-80 ===
    ' low_a = 6378137
    ' low_b = 6356752.314245
    ' f = 1 / 298.257223563
    '
    ' === Airy 1830 ===  Reported best accuracy for England and Northern Europe.
    low_a = 6377563.396
    low_b = 6356256.910
    f = 1 / 299.3249646
    '
    ' === International 1924 ===
    ' low_a = 6378388
    ' low_b = 6356911.946
    ' f = 1 / 297
    '
    ' === Clarke Model 1880 ===
    ' low_a = 6378249.145
    ' low_b = 6356514.86955
    ' f = 1 / 293.465
    '
    ' === GRS-67 ===
    ' low_a = 6378160
    ' low_b = 6356774.719
    ' f = 1 / 298.247167

    '=== WGS-84 Ellipsoid Parameters ===
'    low_a = 6378137       ' +/- 2m
'    low_b = 6356752.3142
'    f = 1 / 298.257223563
    '====================================
    L = ToRad(tLL2.fLng - tLL1.fLng)
    U1 = ATan((1 - f) * Tan(ToRad(tLL1.fLat)))
    U2 = ATan((1 - f) * Tan(ToRad(tLL2.fLat)))
    sinU1 = Sin(U1)
    cosU1 = Cos(U1)
    sinU2 = Sin(U2)
    cosU2 = Cos(U2)
    
    lambda = L
    lambdaP = 2 * cPI
    iterLimit = 100    ' can be set as low as 20 if desired.

    Do While (Abs(lambda - lambdaP) > dEPSILON) And (iterLimit > 0)
        iterLimit = iterLimit - 1

        sinLambda = Sin(lambda)
        cosLambda = Cos(lambda)
        sinSigma = Sqrt(Power(cosU2 * sinLambda, 2) + Power(cosU1 * sinU2 - sinU1 * cosU2 * cosLambda, 2))
        If sinSigma = 0 Then
            Return 0  'co-incident points
         End If
        cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda
        
        '----------------------------------------------------------------------------------
        sigma = ATan2(sinSigma, cosSigma) 'arguments need to be reversed
        'sigma = ATan2(cosSigma, sinSigma) '<<<<<<<<<<<<<<<<<<<<<<<<<<<< original VBA code!
        '----------------------------------------------------------------------------------
        
        sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
        cosSqAlpha = 1 - sinAlpha * sinAlpha

        If cosSqAlpha = 0 Then    'check for a divide by zero
            cos2SigmaM = 0    '2 points on the equator
        Else
            cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha
        End If
        
        C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
        lambdaP = lambda

        'the original calculation is "Too Complex" for Excel VBA to deal with
        'so it is broken into segments to calculate without that issue
        'the original implementation to calculate lambda
        'lambda = L + (1 - C) * f * sinAlpha * _
         '(sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * (cos2SigmaM ^ 2))))
        'calculate portions
        P1 = -1 + 2 * Power(cos2SigmaM, 2)
        P2 = (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * P1))
        'complete the calculation
        lambda = L + (1 - C) * f * sinAlpha * P2
        
    Loop

    If iterLimit < 1 Then
        Log("iteration limit has been reached, something didn't work.")
        Return 0
    End If
    
    uSq = cosSqAlpha * (Power(low_a, 2) - Power(low_b, 2)) / Power(low_b, 2)

    'the original calculation is "Too Complex" for Excel VBA to deal with
    'so it is broken into segments to calculate without that issue
    'the original implementation to calculate upper_A
    'upper_A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
    'calculate one piece of the equation
    P1 = (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
    'complete the calculation
    upper_A = 1 + uSq / 16384 * P1

    'oddly enough, upper_B calculates without any issues - JLatham
    upper_B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))

    'the original calculation is "Too Complex" for Excel VBA to deal with
    'so it is broken into segments to calculate without that issue
    'the original implementation to calculate deltaSigma
    'deltaSigma = upper_B * sinSigma * (cos2SigmaM + upper_B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM ^ 2) _
    ' - upper_B / 6 * cos2SigmaM * (-3 + 4 * sinSigma ^ 2) * (-3 + 4 * cos2SigmaM ^ 2)))
    'calculate pieces of the deltaSigma formula
    'broken into 3 pieces to prevent overflow error that may occur in
    'Excel 2010 64-bit version.
    P1 = (-3 + 4 * Power(sinSigma, 2)) * (-3 + 4 * Power(cos2SigmaM, 2))
    P2 = upper_B * sinSigma
    P3 = (cos2SigmaM + upper_B / 4 * (cosSigma * (-1 + 2 * Power(cos2SigmaM, 2)) - upper_B / 6 * cos2SigmaM * P1))
    'complete the deltaSigma calculation
    deltaSigma = P2 * P3

    'calculate the distance
    s = low_b * upper_A * (sigma - deltaSigma)
    
    'round distance to millimeters
    Return s

End Sub

And for the bearing (line direction) I have this Sub (thanks to Klaus):

Sub GetDirectionOfTwoLatLngs(tLL1 As TMapLatLng, tLL2 As TMapLatLng) As Double
    
    Dim x1 As Double
    Dim y1 As Double
    Dim x2 As Double
    Dim y2 As Double
    Dim dAngle As Double
    
    x1 = tLL1.fLng * CosD(tLL1.fLat)
    y1 = tLL1.flat
    x2 = tLL2.fLng * CosD(tLL2.fLat)
    y2 = tLL2.flat
    
    dAngle = ATan2D((y1 - y2), (x2 - x1)) 'trigonometric angle 3 o'clock, positive clockwise
    dAngle = dAngle + 90 'angle 12 o'clock, positve clockwise
    dAngle = dAngle + 360 'only positive values
    dAngle = dAngle Mod 360 'values between 0 and 360
    
    Return dAngle
End Sub

Tried various coding, for example:

'Supplied map point and both ends of map line
Sub GetDistanceFromRoute(tLLX As TMapLatLng, tLLA As TMapLatLng, tLLB As TMapLatLng) As ResumableSub
    
'   For distance up to a few thousands meters I would simplify the issue from sphere to plane.
'   Then, the issue is pretty simply as an easy triangle calculation can be used:
'    We have points A and B And look for a distance X to line AB. Then:

'    Location a;
'    Location b;
'    Location x;

'    double ax = a.distanceTo(x)
'    double alfa = ((Math.abs(a.bearingTo(b) - a.bearingTo(x))) / 180)    * Math.PI
'    double distance = Math.sin(alfa) * ax

    Dim dDistanceFromRoute As Double
    Dim dBearingDiff As Double
    
    Dim rs As ResumableSub = MapUtilities.DistVincenty(tLLA, tLLX)
    Wait For (rs) Complete (dDistanceLX As Double)
    
    dBearingDiff = (Abs(GetDirectionOfTwoLatLngs(tLLA, tLLB) - GetDirectionOfTwoLatLngs(tLLA, tLLX)) / 180) * dPi
    
    dDistanceFromRoute = Abs(Sin(dBearingDiff)) * dDistanceLX
    
    Return dDistanceFromRoute
    
End Sub

But as said, not got the right results yet.
Any suggestion how this should be coded?

RBS

 

[RB Smissaert]

and, not that it matters much, but is there a reason for retrieving data from the resultset a second time here:



rather than just reusing what was already retrieved eg:

    End If
    tLLA = tLLB    'start of next segment = end of this segment
Loop

Yes, I was aware of this and did exactly the same solution:
tLLA = tLLB

There were some problems (likely not related to this change) and I reverted back as I had already enough trouble.

RBS

 

[emexes]

Need to find the bug now caused by supplying 2 identical lat/lng types to that Sub.

Check whether they're equal (or ridiculously close) and if so, then just return the distance between the line-that-is-a-point, and the target point?

 

[RB Smissaert]

Check whether they're equal (or ridiculously close) and if so, then just return the distance between the line-that-is-a-point, and the target point?

If the 2 types have the same lat/lng values then it should just return 0 (zero), given as a Double.
I added this to the Sub.

RBS

 

[peacemaker]

https://chat.deepseek.com/share/4xim651fjk2ghk15tt

 

[emexes]

If the 2 types have the same lat/lng values then it should just return 0 (zero), given as a Double.

Ah, I thought "that sub" was the sub for finding the closest point on a line segment to a target point, and that the 2 types (variables) were the 2 ends of the line.

 

[RB Smissaert]

Any orthogonal equal-scale coordinate system should work. If the end use is as screen coordinates, then may as well use those.

But if end use is as longitude-latitude coordinates, interfacing to a gps library, then better to convert degrees of longitude to be same size as degrees of latitude by multiplying by cos(latitude), find the closest line point to the target point, then rescale the longitude (x) coordinate back to degrees of latitude by dividing it by the same cos(latitude) factor.

I think I may need to do this.
When I put the X on the AB line then the closest point (as determined in the below Sub) end up below that line.
Will give your Cos(latitude) factor a try.

Sub NearestPointOnline2(iIDX As Int, tLLX As TMapLatLng, tLLA As TMapLatLng, tLLB As TMapLatLng) As ResumableSub 'ignore
    
    Dim tLL As TMapLatLng
    Dim ABx As Double = tLLB.fLng - tLLA.fLng
    Dim ABy As Double = tLLB.fLat - tLLA.fLat
    Dim APx As Double = tLLX.fLng - tLLA.fLng
    Dim APy As Double = tLLX.fLat - tLLA.fLat
    
    Dim ab2 As Double = ABx * ABx + ABy * ABy
    Dim ap_ab As Double = APx * ABx + APy * ABy
    
    Dim t As Double = ap_ab / ab2
    
'    If iIDX = 35 Or iIDX = 34 Then
'        Log("iIDX: " & iIDX & ", t: " & t)
'    End If
    
    ' Limit t to correspond to the segment, if you want the segment, not the infinite line
    If t < 0 Then t = 0
    If t > 1 Then t = 1
    
    Dim cx As Double = tLLX.fLng + t * ABx
    Dim cy As Double = tLLX.fLat + t * ABy
    
    tLL.Initialize
    tLL.fLng = cx
    tLL.fLat = cy
    
    'SetMapLocation(tLL)
    'Sleep(3000)
    
    Return tLL
     
End Sub

RBS

 

[RB Smissaert]

I think I may need to do this.
When I put the X on the AB line then the closest point (as determined in the below Sub) end up below that line.
Will give your Cos(latitude) factor a try.

Sub NearestPointOnline2(iIDX As Int, tLLX As TMapLatLng, tLLA As TMapLatLng, tLLB As TMapLatLng) As ResumableSub 'ignore
   
    Dim tLL As TMapLatLng
    Dim ABx As Double = tLLB.fLng - tLLA.fLng
    Dim ABy As Double = tLLB.fLat - tLLA.fLat
    Dim APx As Double = tLLX.fLng - tLLA.fLng
    Dim APy As Double = tLLX.fLat - tLLA.fLat
   
    Dim ab2 As Double = ABx * ABx + ABy * ABy
    Dim ap_ab As Double = APx * ABx + APy * ABy
   
    Dim t As Double = ap_ab / ab2
   
'    If iIDX = 35 Or iIDX = 34 Then
'        Log("iIDX: " & iIDX & ", t: " & t)
'    End If
   
    ' Limit t to correspond to the segment, if you want the segment, not the infinite line
    If t < 0 Then t = 0
    If t > 1 Then t = 1
   
    Dim cx As Double = tLLX.fLng + t * ABx
    Dim cy As Double = tLLX.fLat + t * ABy
   
    tLL.Initialize
    tLL.fLng = cx
    tLL.fLat = cy
   
    'SetMapLocation(tLL)
    'Sleep(3000)
   
    Return tLL
    
End Sub

RBS

Got this all working fine now with the code given by TILogistics (Sub NearestPointOnline).
Have it slightly altered to work with my Lat/Lng types.
Note there was no need to alter the supplied data, eg do multiply by Cos(Latitude) or convert to screen XY values.

Type TMapLatLng(fLat As Double, fLng As Double)

Sub NearestPointOnLineABFromX(tLLA As TMapLatLng, tLLB As TMapLatLng, tLLX As TMapLatLng) As TMapLatLng
    
    Dim tLL As TMapLatLng
    Dim ABx As Double = tLLB.fLng - tLLA.fLng
    Dim ABy As Double = tLLB.fLat - tLLA.fLat
    Dim APx As Double = tLLX.fLng - tLLA.fLng
    Dim APy As Double = tLLX.fLat - tLLA.fLat
    
    Dim ab2 As Double = ABx * ABx + ABy * ABy
    Dim ap_ab As Double = APx * ABx + APy * ABy
    
    Dim t As Double = ap_ab / ab2
    
    ' Limit t to correspond to the segment, if you want the segment, not the infinite line
    If t < 0 Then t = 0
    If t > 1 Then t = 1
    
    Dim cx As Double = tLLA.fLng + t * ABx
    Dim cy As Double = tLLA.fLat + t * ABy
    
    tLL.fLat = cy
    tLL.fLng = cx
    
    Return tLL
    
End Sub

Having the right point on the line and having fixed the type mismatch bug in Sub DistVincenty I get the distance fine
and all is OK now.

RBS

 

[William Lancee]

For completeness...

I wasn't aware of the fact that returning an 0 from a ResumableSub could cause problems if the return value was meant to be a Double.

I wasn't aware either so I checked and it is probably related to a casting bug:
https://www.b4x.com/android/forum/t...inline-casting-not-working.143870/post-912110

I tested it and the error messages look the same.

 

[emexes]

there was no need to alter the supplied data, eg do multiply by Cos(Latitude)

for your application, probably true - i am guessing that the route segments are short enough that the error doesn't matter

and I suspect that for the most part, the closeness bias cancels out anyway,

plus the closer you are to the equator, the less the difference in size between degrees of longitude and degrees of latitude

one other thing to watch out for is: what happens at longitude wraparound? eg 359 to 0 or 180 to -179

 

[TILogistic]

another solution without using vicenty algorithms, If the UTM zones are different.

Class.

Sub Class_Globals
    Private SemiMajorAxis As Double = 6378137 ' Semi-major axis (WGS84)
    Private Flattening As Double = 1 / 298.257223563 ' flattening
End Sub

'Initializes the object. You can add parameters to this method if needed.
Public Sub Initialize
    
End Sub

Private Sub DegToRad(Deg As Double) As Double
    Return Deg * cPI / 180
End Sub

Private Sub RadToDeg(Rad As Double) As Double
    Return Rad * 180 / cPI
End Sub

Private Sub LatLonToECEF(Lat As Double, Lon As Double, H As Double) As Double()
    Dim e2 As Double = 2 * Flattening - Flattening * Flattening
    Dim phi As Double = DegToRad(Lat)
    Dim lambda As Double = DegToRad(Lon)
    Dim N As Double = SemiMajorAxis / Sqrt(1 - e2 * Sin(phi) * Sin(phi))

    Dim X As Double = (N + h) * Cos(phi) * Cos(lambda)
    Dim Y As Double = (N + h) * Cos(phi) * Sin(lambda)
    Dim Z As Double = (N * (1 - e2) + h) * Sin(phi)

    Return Array As Double(X, Y, Z)
End Sub

Private Sub ECEFToLatLon(X As Double, Y As Double, Z As Double) As Double()
    Dim e2 As Double = 2 * Flattening - Flattening * Flattening
    Dim b As Double = SemiMajorAxis * (1 - Flattening)
    Dim ep As Double = Sqrt((SemiMajorAxis*SemiMajorAxis - b * b) / (b * b))
    Dim p As Double = Sqrt(X * X + Y * Y)
    Dim theta As Double = ATan2(Z * SemiMajorAxis, p * b)
    Dim lon As Double = ATan2(Y, X)
    Dim lat As Double = ATan2(Z + ep*ep * b * Sin(theta) * Sin(theta) * Sin(theta), p - e2 * SemiMajorAxis * Cos(theta) * Cos(theta) * Cos(theta))
    Dim N As Double = SemiMajorAxis / Sqrt(1 - e2 * Sin(lat) * Sin(lat))
    Dim h As Double = p / Cos(lat) - N

    Return Array As Double(RadToDeg(lat), RadToDeg(lon), h)
End Sub

Public Sub GeodeticClosestPoint(LatA As Double, LonA As Double, LatB As Double, LonB As Double, LatP As Double, LonP As Double) As Double()
    Dim h As Double = 0 ' We assume height 0 for all points
    
    Dim A() As Double = LatLonToECEF(LatA, LonA, h)
    Dim B() As Double = LatLonToECEF(LatB, LonB, h)
    Dim P() As Double = LatLonToECEF(LatP, LonP, h)

    Dim AB() As Double = Array As Double(B(0) - A(0), B(1) - A(1), B(2) - A(2))
    Dim AP() As Double = Array As Double(P(0) - A(0), P(1) - A(1), P(2) - A(2))

    Dim ab2 As Double = AB(0) * AB(0) + AB(1) * AB(1) + AB(2) * AB(2)
    
    Dim ap_ab As Double = AP(0) * AB(0) + AP(1) * AB(1) + AP(2) * AB(2)

    Dim t As Double = ap_ab / ab2
    If t < 0 Then t = 0
    If t > 1 Then t = 1

    Dim C() As Double = Array As Double(A(0) + t * AB(0), A(1) + t * AB(1), A(2) + t * AB(2))

    Dim Result() As Double = ECEFToLatLon(C(0), C(1), C(2))
    Return Result ' [latitude, longitude, height]
End Sub

use:

    Dim LatA As Double = -33.45   
    Dim LonA As Double = -70.65   
    Dim LatB As Double = -33.46 
    Dim LonB As Double = -70.64   
    Dim LatP As Double = -33.455 
    Dim LonP As Double = -70.645

    Dim Geo As B4XGeodeticClosestPoint
    Geo.Initialize
    Dim Result() As Double = Geo.GeodeticClosestPoint(LatA, LonA, LatB, LonB, LatP, LonP)
    Log($"Nearest geodetic point: Lat = ${Result(0)}, Lon = ${Result(1)}"$)

 

[TILogistic]

If you want, I have another class that uses Vincenty's algorithms, which are more accurate due to their precision-seeking iterations.

Vincenty's formulae - Wikipedia

en.wikipedia.org

What the class does:
Calculating the closest point on a geodesic line using Vincenty's formulas, you can combine the inverse method with Vincenty's direct method to project a point onto that line.

To find the closest point using Vincenty
Use Vincenty's inverse method to calculate the distance and azimuth between the endpoints A and B, and also between A and point P.

Perform an iterative search or analytical formula to find the point on line AB that minimizes the distance to point P:

Use Vincenty's direct method to calculate the position on the geodesic from A by advancing a distance d.

Adjust d so that the distance between this point and P is minimum.

The point that achieves the minimum distance is the closest point on the geodesic line.

Advantages and disadvantages:
This method uses Vincenty's precision and accuracy for distances and positions on the ellipsoid.

The iterative search may be slower but is accurate.

It is especially useful for large distances or different UTM zones.

 

[emexes]

The Vincenty distance calculation though I added myself from code from the net.

If performance becomes an issue, the first thing to do would be replace Vincenty with something simpler and faster.

For land travel size route segment lengths and distances, a spheroid model will be accurate enough. The mm accuracy of Vincenty will be swamped by the fact that land travel (mostly) occurs at differing radii (altitudes) above sea level.

Even for air travel (aviation flight plan) size route segment lengths and distances, a spheroid model will still be accurate enough for the purpose of measuring the distance, although I suspect the closest-point-on-line-to-target-point method might need some attention, ie change from flat-earth to spheroid-earth.