Solved: what is wrong with this LatLong->Distance code? A Trigonometry "feature". - StephenG1CMZ - 10-12-2016 06:20 PM
Don't be put off by the code - it includes three different implementations all exhibiting the same bug, so you only need to check one...Each one is just a dozen or so lines of trigonometry.
I have coded up three routines for converting a pair of LatLong points to a distance.
HAVERS implements a Haversine spherical Earth calculation
VINCENTYS implements a simplification of Vincenty's formula for a spherical Earth.
These two return practically identical results (trivial rounding in micrometres).
Z-VINCENTY is the normal ellipsoidal formula, perhaps a few km away from the other results.
All 3 exhibit the same bug: Reporting a distance of 1km for an expected distance of 54 km, for example. (The expected result here is from Geoscience Australia).
Since all three implementations are wrong, I assume their is some common error.
I have tried changing my atan2 algorithm and YY XX parameter sequence.
I have tried calculating the longitude difference as abs(LON2-LON1) whereas other sources just subtract. I am aware that VINCENTY has a few superfluous brackets compared to other online sources, but that would not cause all 3 to fail.
I have stared at the code long enough that I need a break before having another look...
Can anyone else see something I have misunderstood in all 3 implementations?
Or recommend a worked example with intermediate values that I can check my code against?
Or is it simply that this application demands more precision than the Prime has?
Code:
LOCAL SP:=" ";
LOCAL EPSILON_V:=3ᴇ−11;
LOCAL LoopSafe:=140;
LOCAL LOOPCOUNT,VALIDITY;
LOCAL WGS84RAm:=6378137.0;
LOCAL WGS84FLATTENINV:=298.257223563;
LOCAL WGS84FLATTEN:=1/WGS84FLATTENINV;
LOCAL RA:=WGS84RAm; //RADIUS EQ
//LOCAL RA:=6378388;//TEST
LOCAL FLATTEN:=WGS84FLATTEN;
LOCAL RB:=IFTE(RA,(1-FLATTEN)*RA,0);//RADIUS POLAR
LOCAL ESS:={RA,RB,FLATTEN};//DEFINE ELLIPSOID SIZE AND SHAPE
LOCAL ESS_UNITS:="m";
D2R(DEGS)
BEGIN
RETURN DEGS*(π/180);
END;
R2D(RADS)
BEGIN
RETURN RADS*(180/π);
END;
ATAN2_UNUSEDYX(YY,XX)
BEGIN //MY ATAN2
IF XX==0 AND YY==0 THEN
RETURN 0; //CONVENIENT;
END;
//-π..π
RETURN ARG(XX+*YY);//BUILT-IN HANDLES MOST
END;
ATAN2YX(YY,XX)
BEGIN
LOCAL RESULT;
IF XX>0 THEN RETURN ATAN(YY/XX); END;
IF YY>0 THEN
RETURN (π/2 - ATAN(XX/YY));
END;
IF YY<0 THEN
RETURN (−π/2 - ATAN(XX/YY));
END;
IF XX<0 THEN
RETURN (ATAN(YY/XX) + π); //PLUSMINUS TBD
END;
IF YY==0 AND XX==0 THEN RETURN (0); END;
//NEVER REACH
PRINT("ATAN2 FAULTY");
END;
HAV(THETA)
//HAVERSINE
BEGIN
LOCAL HS:=SIN(THETA/2);
RETURN HS^2;
END;
REPORT_SHAPE()
BEGIN
LOCAL SP:=" ";
PRINT();
PRINT("RA: "+ESS(1)+SP+ESS_UNITS);
PRINT("RB: "+ESS(2)+SP+ESS_UNITS);
PRINT("f: "+ESS(3));
IF ESS(3)≠0 THEN
PRINT("1/f: "+1/ESS(3));
END;
PRINT("Esc");
WAIT;
END;
EXPORT Z_VINCENTY(LAT1,LON1,LAT2,LON2)
// GEODESY INVERSE PROBLEM ELLIPSOIDAL
// RETURNS DIST IN SAME UNITS AS R
// LAT AND LON IN DEGREES
BEGIN
LOCAL UU1,SINUU1,COSUU1; //U1 AUX LATITUDE
LOCAL UU2,SINUU2,COSUU2; //U2 AUX LATITUDE
LOCAL USQ;
LOCAL AAA,BBB;
LOCAL LONDIFFradians; //L=DIFF IN LON
LOCAL PHI1,PHI2;
LOCAL LAMBDA,SINLAMBDA,COSLAMBDA,LAMBDAOLD;
LOCAL SIGMA,SIN_SIGMA,COS_SIGMA,DELTA_SIGMA,SINSQ_SIGMA; //AUX ARC LENGTH
LOCAL SIN_ALPHA,COSSQ_ALPHA;
LOCAL COS_2SIGMAM,COSSQ_2SIGMAM;
LOCAL CC;
LOCAL BRACKET_TERM;
LOCAL DISTS;
LOCAL AZ_FWD,AZ_REV;
LOCAL RA:=ESS(1);
LOCAL RB:=ESS(2);
LOCAL FLATTEN:=ESS(3);
//REPORT_SHAPE;
// CALC L
LONDIFFradians:=D2R(LON2-LON1); //DIFF IN LON
// CALC U1 U2
PHI1:=D2R(LAT1);
PHI2:=D2R(LAT2);
UU1:=ATAN(1-FLATTEN)*TAN(PHI1);//TBC
UU2:=ATAN(1-FLATTEN)*TAN(PHI2);//TBC
//TO CALC LAMBDA:
// PRE-CALC U1 TRIG ONCEONLY
SINUU1:=SIN(UU1); COSUU1:=COS(UU1);
// PRE-CALC U2 TRIG EACH POSITION
SINUU2:=SIN(UU2); COSUU2:=COS(UU2);
//ITERATE CALCULATING LAMBDA
LOOPCOUNT:=0;
LAMBDA:=LONDIFFradians; //V13: DIFF IN LON ON AUX. FIRST APPROX
//PRINT(); PRINT(LAMBDA);
REPEAT
LOOPCOUNT:=LOOPCOUNT+1;
SINLAMBDA:=SIN(LAMBDA); COSLAMBDA:=COS(LAMBDA); //PRECALC
//V14
BRACKET_TERM:=(COSUU1*SINUU2 - SINUU1*COSUU2*COSLAMBDA);
SINSQ_SIGMA:=(COSUU2*SINLAMBDA)^2 + (BRACKET_TERM)^2;
IF SINSQ_SIGMA==0 THEN
//PRINT("SINσ 0: NO DISTANCE");
//0 DISTANCE MIGHT SUGGEST UNDERFLOW
RETURN {1,0,0,0};//VALID ZERO
END;
SIN_SIGMA:=√(SINSQ_SIGMA);
//V15
COS_SIGMA:=SINUU1*SINUU2 + COSUU1*COSUU2*COSLAMBDA;
//V16
SIGMA:=ATAN2YX(SIN_SIGMA,COS_SIGMA);
//V17
SIN_ALPHA:=(COSUU1*COSUU2*SINLAMBDA)/SIN_SIGMA;
COSSQ_ALPHA:=1-(SIN_ALPHA*SIN_ALPHA);
//V18
IF COSSQ_ALPHA==0 THEN
//PRINT("AVOID COSSQα 0");
COS_2SIGMAM:=0;
ELSE
COS_2SIGMAM:=COS_SIGMA - (2*SINUU1*SINUU2/COSSQ_ALPHA);
END;
CC:=(FLATTEN/16)*COSSQ_ALPHA*(4 + FLATTEN*(4 - 3*COSSQ_ALPHA));
LAMBDAOLD:=LAMBDA;
LAMBDA:=LONDIFFradians + (1-CC)*FLATTEN*SIN_ALPHA*(SIGMA + CC*SIN_SIGMA*(COS_2SIGMAM + CC*COS_SIGMA*(−1+2*(COS_2SIGMAM*COS_2SIGMAM))));
//PRINT(LAMBDA);
UNTIL (ABS(LAMBDA-LAMBDAOLD)≤EPSILON_V) OR (LOOPCOUNT≥LoopSafe);
//PRINT("LOOP DUN");
//V19
USQ:=COSSQ_ALPHA*((RA^2 - RB^2)/(RB^2));
//(V3)
AAA:=1 + (USQ/16384)*(4096 + USQ*(−768 + USQ*(320 - 175*USQ)));
//(V4)
BBB:= (USQ/ 1024)*( 256 + USQ*(−128 + USQ*( 74 - 47*USQ)));
//(V6)
COSSQ_2SIGMAM:=(COS_2SIGMAM)^2;
DELTA_SIGMA:=(1/6)*BBB*COS_2SIGMAM*(−3+4*(SIN_ALPHA)^2)*(−3+4*COSSQ_2SIGMAM) ;
DELTA_SIGMA:=(1/4)*BBB*(COS_SIGMA*(−1+2*COSSQ_2SIGMAM)) - DELTA_SIGMA ;
DELTA_SIGMA:=BBB*SIN_SIGMA*(COS_2SIGMAM + DELTA_SIGMA);
DISTS:=AAA*(SIGMA - DELTA_SIGMA);
//PRINT("ANGULAR FRACTIONAL DIST: "+DISTS/π);
DISTS:=RB*DISTS;
//V20:FORWARD AZIMUTH
AZ_FWD:=ATAN2YX((COSUU2*SINLAMBDA),BRACKET_TERM);
AZ_FWD:=IFTE(AZ_FWD<0,AZ_FWD+2*π,AZ_FWD);//BEARINGS ARE POSITIVE
//V21:REVERSE AZIMUTH
AZ_REV:=ATAN2YX((COSUU1*SINLAMBDA),(−SINUU1*COSUU2 + COSUU1*SINUU2*COSLAMBDA));
AZ_REV:=IFTE(AZ_REV<0,AZ_REV+2*π,AZ_REV);//BEARINGS ARE POSITIVE
DISTS:=ROUND(DISTS,0);//VINCENTY CLAIMS mm,
//RETURN DIST,AZIMUTH,REVERSE AZIMUTH,(LOOPCOUNT:DIAGNOSTIC)
RETURN {DISTS,R2D(AZ_FWD),R2D(AZ_REV),LOOPCOUNT};
END;
VINCENTYS(RA,RB,LAT1,LON1,LAT2,LON2)
//SPHERICAL SIMPLIFICATION: WIKIPEDIA
//USE RR=(2RA+RB)/3
BEGIN
LOCAL NUM,DEN;
LOCAL PHI1,PHI2;
LOCAL SINPHI1,SINPHI2,COSPHI1,COSPHI2;
LOCAL COSDELTALAMBDA,DELTALAMBDA;//DELTALONG rad
LOCAL ANGDIST,DIST;
LOCAL RR:=(2*RA+RB)/3;
PHI1:=D2R(LAT1);
PHI2:=D2R(LAT2);
DELTALAMBDA:=D2R((LON2-LON1));//SHOULD WE ABS?
COSDELTALAMBDA:=COS(DELTALAMBDA);
SINPHI1:=SIN(PHI1); COSPHI1:=COS(PHI1);
SINPHI2:=SIN(PHI2); COSPHI2:=COS(PHI2);
NUM:=√(((COSPHI2*SIN(DELTALAMBDA))^2) + (COSPHI1*SINPHI2 - SINPHI1*COSPHI2*COSDELTALAMBDA)^2);
DEN:=(SINPHI1*SINPHI2 + COSPHI1*COSPHI2*COSDELTALAMBDA);
ANGDIST:=ATAN2YX(NUM,DEN);
//PRINT("ANGDIST "+ANGDIST);
DIST:=RR*ANGDIST;
END;
//HAVERSINE AND VINCENTYSPHERICAL RETURN IDENTICAL RESULTS
//APART FROM VERY TRIVIAL ROUNDING
//SO SERVE TO CHECK EACH OTHER
//(TO BE SURE THAT IS ALWAYS SO:CHECK MIN1 GUARD CONDITIONS ETC MATCH)
HAVERS(LAT1,LON1,LAT2,LON2)
BEGIN
LOCAL RR:=(2*ESS(1)+ESS(2))/3;
LOCAL DLON:=D2R(LON2-LON1);
LOCAL DLAT:=D2R(LAT2-LAT1);
LOCAL AA:=HAV(DLAT) + (COS(D2R(LAT1))*COS(D2R(LAT2))*HAV(DLON));
LOCAL CC:=2*ASIN(MIN(1,√(AA)));//MIN 1 AVOIDS ROUNDING OUTRANGE
LOCAL DD:=RR*CC;
RETURN DD;
END;
ATAN2TEST()
BEGIN
LOCAL AA,BB,RR,RRR;
WHILE 1 DO
AA:=RANDOM();
BB:=RANDOM();
RR:=ATAN2YX(AA,BB);
RRR:=ATAN2YX(AA,BB);
IF ABS(RR-RRR)>1ᴇ−11 THEN
PRINT({RR,RRR});
END;
IF ABS(RR)>π THEN
PRINT(RR/π);
END;
END;
END;
A_VINCENTY(LAT1,LON1,LAT2,LON2)
BEGIN
LOCAL TM,RESULT,DISTS;
PRINT(" ");
PRINT("From LAT "+HMS→(LAT1)+" "+→HMS(LAT1));
PRINT("From LON "+HMS→(LON1)+" "+→HMS(LON1));
PRINT("To LAT "+HMS→(LAT2)+" "+→HMS(LAT2));
PRINT("To LON "+HMS→(LON2)+" "+→HMS(LON2));
TM:=TICKS;
RESULT:=HAVERS(LAT1,LON1,LAT2,LON2);
PRINT("HAVERS: "+(RESULT/1000)+" km");
RESULT:=VINCENTYS(ESS(1),ESS(2),LAT1,LON1,LAT2,LON2);
PRINT("SPHERICAL: "+(RESULT/1000)+" km");
RESULT:=Z_VINCENTY(LAT1,LON1,LAT2,LON2);
TM:=TICKS-TM;
//PRINT("Elapsed "+(TM/1000)+" s");
PRINT("VALIDITY "+(LOOPCOUNT<LoopSafe)+" "+RESULT(4));//BOOLEAN, LOOPCOUNT (LOWER=>BETTER CONVERGENCE)
DISTS:=RESULT(1);
PRINT("DISTANCE: "+(DISTS/1000)+SP+" km");
//PRINT(" AZ: "+→HMS(RESULT(2)));
//PRINT(" AZ: "+→HMS(RESULT(3)));
END;
EXPORT STEPPED()
BEGIN
LOCAL II;
FOR II FROM 1 TO 89 STEP 30 DO
A_VINCENTY(−II,0,II,0);
END;
END;
EXPORT VINCENTY()
BEGIN
//
//ATAN2TEST();
//PRINT("TESTED"); WAIT;
REPORT_SHAPE();
A_VINCENTY(50°03′58.76″,−005°42′53.10″,58°38′38.48″,−003°04′12.34″);//50PERCENT
PRINT("EXPECT 969 954_m");
A_VINCENTY(−37°57′3.72030″,144°25′29.52440″,−37°39′10.15610″,143°55′35.38390″);//EXPECTS 54972m
PRINT("EXPECTED: 54972m");//VINCENTY AUSTRALIA
A_VINCENTY(31.8300167,35.0662833,31.83000000,35.0708167);
PRINT("PYTHON: CIRCLE 428.4 VIN 429.1677 m");
//A_VINCENTY(0,0,0.5,179.5);
//A_VINCENTY(0,0,0.5,179.7);//WILL NOT CONVERGE
//A_VINCENTY(0,0,0,90);
//A_VINCENTY(0,0,0.5,0.5);
//A_VINCENTY(2,2,2,2);
STEPPED();
//PRINT("SEMI: "+π*ESS(1));//MAX
END;
RE: Can you see what is wrong with this LatLong->Distance code? - Guenter Schink - 10-12-2016 08:14 PM
It's rather difficult for me to decipher your code entirely. Therefore I've just taken the Haversine part and adjusted it a little bit. That is, I've done the HAV-function and the D2R in place. Further I replaced RR directly with the value of 6378.137 Km, I think that is correct because Haversine refers to the sphere and not to the ellipsoid.
Code:
EXPORT Haversine(LAT1,LON1,LAT2,LON2)
BEGIN
LOCAL RR:=6378.137;
LOCAL DLON:=(LON2-LON1)/180*Pi;
LOCAL DLAT:=(LAT2-LAT1)/180*Pi;
LOCAL AA:=(sin(DLAT/2))^2 +
COS((LAT1)/180*Pi)*
COS((LAT2)/180*Pi)*(sin(DLON/2))^2;
LOCAL CC:=2*ASIN(MIN(1,√(AA)));//MIN 1 AVOIDS ROUNDING OUTRANGE
LOCAL DD:=RR*CC;
RETURN DD;
END;
This works for me with known values. But I must admit, for rather long distances. If it works for you also then you should investigate the method of achieving "RR" with your "ESS{}". There might be the glitch.
HTH Günter
RE: Can you see what is wrong with this LatLong->Distance code? - StephenG1CMZ - 10-12-2016 09:46 PM
Thank you for taking the time to look at this.
When I first saw your code I was thinking the lower-case Sin was exact and avoiding a rounding error - but it doesn't fix this example, which is still 1 km rather than 54 km.
For spherical calculations I usually use an average of polar and equatorial radius, (actually. (2a+b)/3), or repeat the calculations for both - but that detail isn't significant.
Example data:
http://www.ga.gov.au/geodesy/datums/vincenty_inverse.jsp
Update: a missing atan2 call is added - the reported distance is still 1 km. Actually, that is probably just equivalent to the previous line.
Code:
EXPORT Haversine(LAT1,LON1,LAT2,LON2)
BEGIN
LOCAL RR:=6378.137;
LOCAL DLON:=(LON2-LON1)/180*Pi;
LOCAL DLAT:=(LAT2-LAT1)/180*Pi;
LOCAL AA:=(sin(DLAT/2))^2 +
COS((LAT1)/180*Pi)*
COS((LAT2)/180*Pi)*(sin(DLON/2))^2;
LOCAL CC:=2*ASIN(MIN(1,√(AA)));//MIN 1 AVOIDS ROUNDING OUTRANGE
//INSERT LINES
LOCAL DD;
CC:=2*ATAN2YX(√(AA),√(1-AA^2));
//END INSERT
DD:=RR*CC;//REMOVE LOCAL
RETURN DD;
END;
EXPORT HAVERSINES()
BEGIN
LOCAL DIST;
PRINT();
DIST:=Haversine(−37°57′3.72030″,144°25′29.52440″,−37°39′10.15610″,143°55′35.38390″);//EXPECTS 54972m
PRINT(DIST+" km");
PRINT("EXPECTED: 54972m");//VINCENTY AUSTRALIA
END;
RE: Can you see what is wrong with this LatLong->Distance code? - Guenter Schink - 10-13-2016 10:12 AM
(10-12-2016 09:46 PM)StephenG1CMZ Wrote: Thank you for taking the time to look at this.
When I first saw your code I was thinking the lower-case Sin was exact and avoiding a rounding error - but it doesn't fix this example, which is still 1 km rather than 54 km.
Sorry the lower case "sin" was not intentionally.
But now I think we're coming closer. Don't know what really happens (a bug?). I added commands to print DLON and DLAT. And much to my surprise they showed up in DMS format. Consequently I removed the transformation to radians et voilá the result comes pretty close to what's expected.
Code:
EXPORT Haversine(LAT1,LON1,LAT2,LON2)
BEGIN
LOCAL RR:=6378.137;
LOCAL DLON:=(LON2-LON1);
LOCAL DLAT:=(LAT2-LAT1);
LOCAL AA:=(SIN(DLAT/2))^2 +
COS(LAT1)*
COS(LAT2)*(SIN(DLON/2))^2;
LOCAL CC:=2*ASIN(MIN(1,√(AA)));//MIN 1 AVOIDS ROUNDING OUTRANGE
LOCAL DD:=RR*CC;
PRINT("DLAT " + DLAT);
PRINT("DLON " + DLON);
RETURN ("DD " + DD);
END;
That's really weird. The "Angle Measure" is set to radians! I guess you have it too. The problem seems not to show up when entering the co-ordinates not in DMS but degrees (with decimal degrees)**. On the other hand, I don't understand, at least in the Haversine formula, why a transformation to radians is needed after all.
It will happen on the command line too. Just enter "12°30' * PI / 180" the expression will stay in DMS mode.
HTH, Günter
** That's obviously the reason why I didn't encounter the problem, as I entered the data in degrees format.
Günter
RE: Can you see what is wrong with this LatLong->Distance code? - StephenG1CMZ - 10-13-2016 04:52 PM
Thank you so much for that insight Guenter - I have only briefly read your answer and need to re-read it more thoroughly.
But if I understand you, these two SIN's give different answers, despite - one would have thought - being the same calculation (my home setting is Radians default).
Code:
EXPORT HMSBUG()
BEGIN
IF SIN(12°30′″)≠SIN(HMS→(12°′30″)) THEN
PRINT("THAT SEEMS WRONG");
ELSE
PRINT("SIN WORKS");
END;
END;
My understanding of HMS is that it changed how a value was displayed, whether 12.5 or 12:30 - not that it changed the resulting answers! That is so unexpected to me that I would consider it a bug. (The answers differ completely, not just a rounding of the presentation value).
RE: Can you see what is wrong with this LatLong->Distance code? Trigonometry bug? - Tim Wessman - 10-13-2016 05:31 PM
If you are specifying an input in DMS then your value being input is degrees - essentially you have tagged it with a unit. What would be the expected behavior for these cases?
SIN(1.5)
SIN(1.5_grad)
SIN(1.5_rad)
SIN(1°30')
Seems to me the first would be dependent on your setting and all the rest are essentially units that can be directly handled no matter your system setting to give the correct result. Since DMS is an inherent part of the number, you don't have to do conversions between steps to calculate using DMS. Thus something like 5°15' * 3 is a perfectly valid calculation and returns a DMS result.
Anyway, I think you have your answer about the bug. If you want to be completely independent of your setting, run HMS-> (doesn't matter if you are are in HMS or not, just sets the flag on the real number input) then D->R.
I'm wondering now though if HMS-> and ->HMS should also handle angular unit objects directly so something like ->HMS(1.5_rad) would return 85°56′37.20937″ ... thoughts anyone?
RE: Can you see what is wrong with this LatLong->Distance code? Trigonometry bug? - StephenG1CMZ - 10-13-2016 08:58 PM
I guess where I was wrong in thinking it should "just work", is I wasn't thinking of the DMS format as including a unit. After all, it was degrees, not underscore degrees to be consistent with other units - and the same syntax could also have been hours.
I was thinking of it as more of just a different presentation of a decimal number, and assuming the calculation would be unchanged. Like 1.5 = 1,5.. Like decimal 5 = hex 5.
Seeing your other examples with units, I see it's not as simple as keeping the presentation and calculation separate.
As it stands, this issue seems likely to catch people out unless they are aware of it and code around it. But I am not sure what could or should be changed to make things clearer.
RE: Can you see what is wrong with this LatLong->Distance code? Trigonometry bug? - Guenter Schink - 10-14-2016 09:53 AM
(10-13-2016 05:31 PM)Tim Wessman Wrote: .
Anyway, I think you have your answer about the bug. If you want to be completely independent of your setting, run HMS-> (doesn't matter if you are are in HMS or not, just sets the flag on the real number input) then D->R.
I'm wondering now though if HMS-> and ->HMS should also handle angular unit objects directly so something like ->HMS(1.5_rad) would return 85°56′37.20937″ ... thoughts anyone?
Contrary to my previous statement, with a question mark though, I would no longer consider this to be a bug. In my opinion it is a feature, perhaps needs to be stated clearer somewhere! But where? And a, yet to be implemented, D->R command should work regardless of DMS or HR format, avoiding these ambiguities.
And I think your proposal for ->HMS () would be valuable enhancement.
Günter
Don't forget to implement newRPL
RE: Can you see what is wrong with this LatLong->Distance code? A Trigonometry &quo... - StephenG1CMZ - 10-17-2016 10:17 PM
I have now uploaded corrected spherical calculations to my Geodesy program in the software library. The ellipsoidal Vincenty implementation needs more work, however - it's values are currently less accurate than the spherical calculations.
RE: Solved: what is wrong with this LatLong->Distance code? A Trigonometry "fe... - StephenG1CMZ - 10-23-2016 06:41 PM
Version 0.3 of my Geodesy program now implements Lat Long to distance calculations, including the Vincenty formula. If you have V0.1 or V0.2 please upgrade.
|