Android Question Get shortest distance from point to line on map

R

RB Smissaert

Well-Known Member
Licensed User
Longtime User
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:

Code: 
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):

Code: 
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:

Code: 
'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
 
Sort by date Sort by votes
RB Smissaert said:
For distance up To a few thousands meters I would simplify the issue from sphere To plane.
Click to expand...

+1

and depending on how accurate the answer needs to be, I'd be considering up to a few degrees ie hundreds of km.

Another sometimes-useful simplification of doing stuff like this is to convert the latitude-longitude-radius into x-y-z.

But using just x-y on a plane, I'm pretty sure it can be done with arithmetic, no trig required, possibly even no square-root required for pythagoras if you just need the closest point and don't need to know the actual distance.

First thing would be that the closest point on the line to the point is one of:
1/ a
2/ b
3/ somewhere in between

If it is somewhere-in-between then it will be where a line at right-angle to line a-b that passes through location x, intersects with the line a-b.

If the line equation is (re)stated as ux + vy = k then I'm pretty sure that swapping u and v will give a right-angle (perpendicular) line to that, and then it's a case of substituting-in the x and y of location x to find k, and then seeing if that specific perpendicular line passes through line a-b between points a and b.

Lol feels easy enough, which is a red flag. I'll try it in a spreadsheet to find out what I've missed.
 
Upvote 0
emexes said:
I'll try it in a spreadsheet
Click to expand...

Lol I ended up doing a worked example on paper with a diagram so that I could see what was happening, figured out that the line needed to be restated as two equations and a "sliding" variable that I called i, that varied between value 0 at one end of the line, to value 1 at the other end of the line. Then swap the x and y factors and negate one of them, to do the 90 degree perpendicularisation (is that a word?) and then I thought, heck, why don't I just ask Grok to finish it for me?


I haven't checked it manually myself, but it looks correct, eg Grok also had a "sliding" variable but called it t instead of i.

Also Grok tightened up some of my phrasing, eg it's a line segment, not a line. Go Grok!!!
 
Last edited:
Upvote 0
emexes said:
Lol I ended up doing a worked example on paper with a diagram so that I could see what was happening, figured out that the line needed to be restated as two equations and a "sliding" variable that I called i, that varied between value 0 at one end of the line, to value 1 at the other end of the line. Then swap the x and y factors and negate one of them, to do the 90 degree perpendicularisation (is that a word?) and then I thought, heck, why don't I just ask Grok to finish it for me?


I haven't checked it manually myself, but it looks correct, eg Grok also had a "sliding" variable but called it t instead of i.

Also Grok tightened up some of my phrasing, eg it's a line segment, not a line. Go Grok!!!
Click to expand...
This approach:

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

Looks a lot simpler than your Grok approach.
Shame it doesn't seem to give the right results!
I think it is the right approach though to think about where the distance line should join the A to B line rather than trying to get the distance directly.
Will start (as you did) on a piece of paper.
Will work it out.

RBS
 
Upvote 0
= Haversine(latitus,longitus, Lat2s, Lon2s)

Distancias2 return the distance

Where latitus,longitus are the the inicial spot
and Lat2s, Lon2s will be the gps reading

Use 00.000000 its very accurate


Try This one

Code: 
Public  Sub Haversine(lat1 As Double, lon1 As Double, lat2 As Double, lon2 As Double) As Double
   
    Dim lat1Rad As Double = DegreesToRadians(lat1)
    Dim lon1Rad As Double = DegreesToRadians(lon1)
    Dim lat2Rad As Double = DegreesToRadians(lat2)
    Dim lon2Rad As Double = DegreesToRadians(lon2)
 
    ' Haversine formula
    Dim dlon As Double = lon2Rad - lon1Rad
    Dim dlat As Double = lat2Rad - lat1Rad
    Dim a As Double = Power(Sin(dlat / 2), 2) + Cos(lat1Rad) * Cos(lat2Rad) * Power(Sin(dlon / 2), 2)
    Dim c As Double = 2 * ATan2(Sqrt(a), Sqrt(1 - a))
   Dim radiusOfEarth As Double] = 6717407                 '<------------- Meter/Miles/KMS
    Dim distance As Double = radiusOfEarth * c
'    Log(distance)
    Return distance
End Sub

Meter = 6717407
Miles 6378.137
KMS 6371

This works pretty good

Salud
 
Last edited:
Upvote 0
Situ LLC said:
= Haversine(latitus,longitus, Lat2s, Lon2s)

Distancias2 return the distance

Where latitus,longitus are the the inicial spot
and Lat2s, Lon2s will be the gps reading

Use 00.000000 its very accurate


Try This one

Code: 
Public  Sub Haversine(lat1 As Double, lon1 As Double, lat2 As Double, lon2 As Double) As Double
  
    Dim lat1Rad As Double = DegreesToRadians(lat1)
    Dim lon1Rad As Double = DegreesToRadians(lon1)
    Dim lat2Rad As Double = DegreesToRadians(lat2)
    Dim lon2Rad As Double = DegreesToRadians(lon2)
 
    ' Haversine formula
    Dim dlon As Double = lon2Rad - lon1Rad
    Dim dlat As Double = lat2Rad - lat1Rad
    Dim a As Double = Power(Sin(dlat / 2), 2) + Cos(lat1Rad) * Cos(lat2Rad) * Power(Sin(dlon / 2), 2)
    Dim c As Double = 2 * ATan2(Sqrt(a), Sqrt(1 - a))
   Dim radiusOfEarth As Double] = 6717407                 '<------------- Meter/Miles/KMS
    Dim distance As Double = radiusOfEarth * c
'    Log(distance)
    Return distance
End Sub

Meter = 6717407
Miles 6378.137
KMS 6371

This works pretty good

Salud
Click to expand...
This seems to the distance between 2 points.
I need the shortest distance between a point and a line.
Or am I mistaken?

RBS
 
Upvote 0
Well I use only have between 2 points.
SALUD
 
Upvote 0
RB Smissaert said:
Looks a lot simpler than your Grok approach.
Shame it doesn't seem to give the right results!
Click to expand...

Did the Grok approach not work? I agree it has more operations, but they're also faster, ie 10-to-20 times faster than trig operations, so it's probably a wash from a difficulty and performance perspective.

Anyway, using trig, and given points a, b and p, forming a triangle, I can see that angles abp and bap both have to be less than 90 degrees, if the closest point on the line is to be within that line segment ie within a and b. Otherwise, the closest point on the line is either a or b, and you'd just calculate those two distances and take the lowest.

Hmm then I thought what if those two distances are the same, which would you choose, a or b? But... if they're equal, then the closest point on the line would be the mid-point of the line ie (a + b) / 2

Actually, now I look at the Grok approach, I'm thinking the double-letter-variables-that-look-like-multiplications are making it look harder than it is.

It uses:

- 4 subtractions (looks like 6, but 2 of them are repeats)
- 4 multiplications
- 2 additions
- 1 division

and then if the closest point on the line is between a and b (aka if t is > 0 and < 1) then:

- no more subtractions (assuming we saved them from first stage)
- 2 more multiplications
- 2 more additions

so if that approach actually comes up with the right answers, then I'm leaning towards the devil we know solution we understand rather than the solution we don't.
 
Last edited:
Upvote 0
ChatGPT came up with a similar approach, so that bodes well:

chatgpt.com

ChatGPT - Closest point formula

Shared via ChatGPT
chatgpt.com chatgpt.com

I like that it explicitly does the repeated subtractions only once, saving them as dx and dy, which makes the subsequent equations look less scary complicated.
 
Upvote 0
CoPilot will even generate a sketch too (albeit with unequal x and y axis scales, wtf?) (although it's nice to see AI can still shoot itself in the foot)
1760168605511.png
 
Last edited:
Upvote 0

Distance measurement of 2 GPS positions with B4J

Hello everyone, I wrote a small program to measure the distance of 2 GPS coordinates. As a B4A version it works perfectly. How do I have to rewrite this for B4J because there is no package android.location? I am attaching both as a zip file. Thanks for your help.
www.b4x.com www.b4x.com
 
Upvote 0
1760173553201.png

Code: 
Sub NearestPointOnline(ax As Double, ay As Double, bx As Double, by As Double, px As Double, py As Double) As Float()
    Dim ABx As Double = bx - ax
    Dim ABy As Double = by - ay
    Dim APx As Double = px - ax
    Dim APy As Double = py - ay
    
    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 = ax + t * ABx
    Dim cy As Double = ay + t * ABy
    
    Return Array As Float(cx, cy)
End Sub
 
Upvote 0
Thanks, will give that a try.
I work with lat/lng types and I take it your code wouldn't work with that, so I take it I have to convert to map points (x, y types) indicating
screen coordinates.

RBS
 
Upvote 0
RB Smissaert said:
have to convert to map points (x, y types) indicating screen coordinates.
Click to expand...

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.
 
Upvote 0
emexes said:
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.
Click to expand...
I think it can work on lat/lng coordinates. When I tried it on a single set (point X and line points A and B, all set by lat/lng types), using the posted code from TILogistic, then
it works and I get the right (shortest) distance from point X to line A_B.
However when I run this in a loop, looping through a recordset, holding lat and lng values, I get an error, which sofar I have been unable to track and explain:

p_map$ResumableSub_GetDistanceFromRoute3resume (java line: 8242)
java.lang.ClassCastException: java.lang.Integer cannot be cast to java.lang.Double
at b4a.exampleljjll.p_map$ResumableSub_GetDistanceFromRoute3.resume(p_map.java:8242)
at anywheresoftware.b4a.BA.checkAndRunWaitForEvent(BA.java:275)
at anywheresoftware.b4a.BA.raiseEvent2(BA.java:215)
at anywheresoftware.b4a.BA.raiseEvent(BA.java:201)
at anywheresoftware.b4a.keywords.Common$15.run(Common.java:1857)
at android.os.Handler.handleCallback(Handler.java:995)
at android.os.Handler.dispatchMessage(Handler.java:103)
at android.os.Looper.loopOnce(Looper.java:273)
at android.os.Looper.loop(Looper.java:363)
at android.app.ActivityThread.main(ActivityThread.java:10060)
at java.lang.reflect.Method.invoke(Native Method)
at com.android.internal.os.RuntimeInit$MethodAndArgsCaller.run(RuntimeInit.java:632)
at com.android.internal.os.ZygoteInit.main(ZygoteInit.java:975)
Starter, LogError, strErrorMsg: java.lang.ClassCastException: java.lang.Integer cannot be cast to java.lang.Double

Code: 
Sub GetDistanceFromRoute3(iIDX As Int, tLLX As TMapLatLng, tLLA As TMapLatLng, tLLB As TMapLatLng) As ResumableSub
    
    Dim rs1 As ResumableSub = NearestPointOnline2(iIDX, tLLX, tLLA, tLLB)
    Wait For (rs1) Complete (tLL_AB As TMapLatLng)
    
    'SetMapLocation(tLL_AB)
    'Log(iIDX & ": " & tLL_AB.fLat & "," & tLL_AB.fLng)
    'Sleep(3000)
    
    Dim rs2 As ResumableSub = MapUtilities.DistVincenty(tLL_AB, tLLX, "GetDistanceFromRoute3")
    Wait For (rs2) Complete (dDistanceLX As Double)

    
'    If iIDX = 34 Or iIDX = 35 Then
'        Log(iIDX & ": " & iIDX & ", distance: " & dDistanceLX)
'    End If
    
    Return dDistanceLX

End Sub

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

Also I get too short distances when the lines are not near the point X (but all less than 1500 meters).
Maybe it does indeed need the lat/lng scaling as you suggested.
Will investigate further.

RBS
 
Upvote 0
RB Smissaert said:
I think it can work on lat/lng coordinates. When I tried it on a single set (point X and line points A and B, all set by lat/lng types), using the posted code from TILogistic, then
it works and I get the right (shortest) distance from point X to line A_B.
However when I run this in a loop, looping through a recordset, holding lat and lng values, I get an error, which sofar I have been unable to track and explain:

p_map$ResumableSub_GetDistanceFromRoute3resume (java line: 8242)
java.lang.ClassCastException: java.lang.Integer cannot be cast to java.lang.Double
at b4a.exampleljjll.p_map$ResumableSub_GetDistanceFromRoute3.resume(p_map.java:8242)
at anywheresoftware.b4a.BA.checkAndRunWaitForEvent(BA.java:275)
at anywheresoftware.b4a.BA.raiseEvent2(BA.java:215)
at anywheresoftware.b4a.BA.raiseEvent(BA.java:201)
at anywheresoftware.b4a.keywords.Common$15.run(Common.java:1857)
at android.os.Handler.handleCallback(Handler.java:995)
at android.os.Handler.dispatchMessage(Handler.java:103)
at android.os.Looper.loopOnce(Looper.java:273)
at android.os.Looper.loop(Looper.java:363)
at android.app.ActivityThread.main(ActivityThread.java:10060)
at java.lang.reflect.Method.invoke(Native Method)
at com.android.internal.os.RuntimeInit$MethodAndArgsCaller.run(RuntimeInit.java:632)
at com.android.internal.os.ZygoteInit.main(ZygoteInit.java:975)
Starter, LogError, strErrorMsg: java.lang.ClassCastException: java.lang.Integer cannot be cast to java.lang.Double

Code: 
Sub GetDistanceFromRoute3(iIDX As Int, tLLX As TMapLatLng, tLLA As TMapLatLng, tLLB As TMapLatLng) As ResumableSub
   
    Dim rs1 As ResumableSub = NearestPointOnline2(iIDX, tLLX, tLLA, tLLB)
    Wait For (rs1) Complete (tLL_AB As TMapLatLng)
   
    'SetMapLocation(tLL_AB)
    'Log(iIDX & ": " & tLL_AB.fLat & "," & tLL_AB.fLng)
    'Sleep(3000)
   
    Dim rs2 As ResumableSub = MapUtilities.DistVincenty(tLL_AB, tLLX, "GetDistanceFromRoute3")
    Wait For (rs2) Complete (dDistanceLX As Double)

   
'    If iIDX = 34 Or iIDX = 35 Then
'        Log(iIDX & ": " & iIDX & ", distance: " & dDistanceLX)
'    End If
   
    Return dDistanceLX

End Sub

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

Also I get too short distances when the lines are not near the point X (but all less than 1500 meters).
Maybe it does indeed need the lat/lng scaling as you suggested.
Will investigate further.

RBS
Click to expand...
Forgot to post the main, initial Sub:

Code: 
Sub GetShortestDistanceFromRoute(iRoute_ID As Int, tLLX As TMapLatLng) As ResumableSub
    
    Dim strSQL As String
    Dim oRS As ResultSet
    Dim tLLA As TMapLatLng
    Dim tLLB As TMapLatLng
    Dim dDistance As Double
    
    strSQL = "select route_idx, latitude, longitude from route_data where route_id = ? order by route_idx asc"
    oRS = cMP.cConn.SQLNon_Clinical.ExecQuery2(strSQL, Array As String(iRoute_ID))
    oRS.Position = 0
    
    'Log("oRS.RowCount: " & oRS.RowCount)
    
    tLLA.fLat = oRS.GetDouble2(1)
    tLLA.fLng = oRS.GetDouble2(2)
    
    dDistance = 1000000
    
    Do While oRS.NextRow
        tLLB.fLat = oRS.GetDouble2(1)
        tLLB.fLng = oRS.GetDouble2(2)
        Dim rs As ResumableSub = GetDistanceFromRoute3(oRS.GetInt2(0), tLLX, tLLB, tLLA)
        'Dim rs As ResumableSub = GetDistanceFromRoute3(oRS.GetInt2(0), tLLX, tLLA, tLLB)
        Wait For (rs) Complete (dDistanceLX As Double)
        Log("dDistanceLX: " & dDistanceLX & ", Route_IDX: " & oRS.GetInt2(0))
        If dDistanceLX < dDistance Then
            dDistance = dDistanceLX
            'Log("dDistance: " & dDistance & ", Route_IDX: " & oRS.GetInt2(0))
        End If
        tLLA.fLat = oRS.GetDouble2(1)
        tLLA.fLng = oRS.GetDouble2(2)
    Loop
    
    oRS.Close
    
    Return dDistance
    
End Sub

RBS
 
Upvote 0
RB Smissaert said:
p_map$ResumableSub_GetDistanceFromRoute3resume (java line: 8242)
java.lang.ClassCastException: java.lang.Integer cannot be cast to java.lang.Double
Click to expand...

Hmm well "cannot be cast to" looks like a type mismatch, but... I am with you, also can't see anything in your code that would cause that.

First step would be to add some logging to try isolate which ResumableSub "return" or "resume" is the problem.

Also, see if the Java line 8242 is within the B4A transcompiler output Java source code in or under the project's Objects folder.

RB Smissaert said:
Dim rs2 As ResumableSub = MapUtilities.DistVincenty(tLL_AB, tLLX, "GetDistanceFromRoute3")
Wait For (rs2) Complete (dDistanceLX As Double)
Click to expand...

Which library is MapUtilities in? Do you have a download link or page for it?
 
Last edited:
Upvote 0
and, not that it matters much, but is there a reason for retrieving data from the resultset a second time here:

RB Smissaert said:
Code: 
    End If
    tLLA.fLat = oRS.GetDouble2(1)
    tLLA.fLng = oRS.GetDouble2(2)
Loop
Click to expand...

rather than just reusing what was already retrieved eg:

Code: 
    End If
    tLLA = tLLB    'start of next segment = end of this segment
Loop
 
Upvote 0
emexes said:
Hmm well "cannot be cast to" looks like a type mismatch, but... I am with you, also can't see anything in your code that would cause that.

First step would be to add some logging to try isolate which ResumableSub "return" or "resume" is the problem.

Also, see if the Java line 8242 is within the B4A transcompiler output Java source code in or under the project's Objects folder.



Which library is MapUtilities in? Do you have a download link or page for it?
Click to expand...
I found the cause of the type mismatch in Sub DistVincenty.
In 2 places it could return 0 (zero).

This is the adjusted Sub and that gives no errors now.
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.

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

Code: 
Public Sub DistVincenty(tLL1 As TMapLatLng, tLL2 As TMapLatLng) As ResumableSub
        
    'Log(tLL1.fLat & "," & tLL1.fLng & " _ " & tLL2.fLat & "," & tLL2.fLng)
    '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
    Dim d As Double = 0 '>>> Just to avoid error: java.lang.ClassCastException: java.lang.Integer cannot be cast to java.lang.Double
    
    'Added for simplicity
    '--------------------
    If tLL1.fLat = tLL2.fLat Then
        If tLL1.fLng = tLL2.fLng Then
            Return d
        End If
    End If
    
'    If strCallingSub <> "GetDistanceFromRoute3" Then
'        Log("Early return DistVincenty")
'        Return s
'    End If

    '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 d  '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 d
    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

>> Which library is MapUtilities in?
This is part of an OSM map library. The Vincenty distance calculation though I added myself from code from the net.
 
Upvote 0
You must log in or register to reply here.

Similar Threads

P
Another implementation of Dijkstra's shortest distance algorithm
Replies
2
Views
4K
carlos7000
carlos7000
Pedro Caldeira
  • Question Question
Android Question Google Maps shortest distance
Replies
6
Views
4K
TILogistic
TILogistic
D
  • Article Article
Android Code Snippet Measure distance between 2 LatLng points
Replies
3
Views
4K
klaus
klaus
Beja
  • Question Question
Android Question GPS: shortest route
Replies
7
Views
4K
Beja
Beja
F
  • Question Question
Measurement of Distances
Replies
14
Views
14K
IslamQabel
I
Back
Top