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:
And for the bearing (line direction) I have this Sub (thanks to Klaus):
Tried various coding, for example:
But as said, not got the right results yet.
Any suggestion how this should be coded?
RBS
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:
B4X:
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 SubAnd for the bearing (line direction) I have this Sub (thanks to Klaus):
B4X:
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 SubTried various coding, for example:
B4X:
'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 SubBut as said, not got the right results yet.
Any suggestion how this should be coded?
RBS
