Beta Test / Technical

[manimaul]

interesting, aren't we using some bits of code from proj.4 already? i'm interested to hear Dave's thoughts on this when he gets back in 2-3 weeks. it will be nice to get this cleaned up so we can start rolling opencpn into Debian.

[HamishB]

I just did a fn by fn review of georef.c in the sf bug tracker ticket.
a lot of it is ok from public domain sources. but great circle stuff will have to be replaced AFAICT as the originating Aviation calc website does not state an explicit license. (we could always email the author and ask him to put up a statement...)

as mentioned geod.c and geodesic.c from proj4/gerald are a good choice.

also, some other sources:

from GpsDrive: calc_wpdist()
http://gpsdrive.svn.sourceforge.net/viewvc/gpsdrive/trunk/src/gpsmisc.c

from GRASS:
https://trac.osgeo.org/grass/browser/grass/trunk/lib/gis/geodesic.c

roll your own:
http://trac.osgeo.org/proj/wiki/GeodesicCalculations



grep doesn't show any pj_*() functions from PROJ.4 outside of src/mygdal/, but GDAL uses PROJ.4 heavily. (Frank W. is the maintainer of both projects)


Hamish

[cagney]

Originally Posted by PjotrC

I think Mr. Mercator may have a problem with small values.
The original 1.3.6 too, but it is easier to see here.
Until anchor chain lengths came into the picture it was no issue...

Piotr

I think that you notice problems with short distances especially in an E - W direction. This is due to the implementation.
There is some problems implementing Mercartor course/distance near the E - W line. To get away from this Dave implemented GreatCircle sailing for everything, then..I came along with my GreatCircle GPX script, and what looked like Mercator distance was equal to the GreatCircle distance, ..a bit confusing. This was in last December before 1.3.6. So Dave reverted to pure Mercator Sailing and implemented a hack to get away from the E - W problem. This is showing up now when the question of measuring very short distances has popped up. Hopefully changing to one of the libraries that Hamish mentions, will solve this.
I don't think that 1.3.6 uses GC calculations at all, any code is probably just left overs from the earlier implementation. But I've been wrong before .

Thomas

[PjotrC]

Originally Posted by cagney

I think that you notice problems with short distances especially in an E - W direction.

Thomas,

Thanks for the explanation. Yes, it concerns E-W direction and very short distances. I have nearby a speed-measurement range and the beacons are obviously in place.

I was wondering why in the automatic anchorage marking a planar formula was used, and although it was OK there, it was obviously unsuitable for anchor watch check.

Piotr

[brak]

Originally Posted by PjotrC

I'm at latitude 54 N. I think that the difference in measurement in exactly 270 bearing differs too much from that of 283. The bearing 001 looks OK. With longer distances, say 5 miles, it is not noticeable, but even before, looking at tracks in close zoom, I had the impression of some unexplained differences between the chart and textual position representation.

Therefore my vague request for similar experience.

True scale is 800, zoom 14x (overzoom), but it is necessary, if anchor watch and tracking should be useful. In this situation, I am almost not interested in the chart, just the behaviour of the boat at anchor.

Piotr

Could it be this bit from the DistanceBearingMercator?

// This idea comes from Thomas(Cagney)
// We simply require the dlat to be (slightly) non-zero, and carry on.

mlat0 = lat0;
if(fabs(lat1 - lat0) < 1e-4)
mlat0 += .001;

Adding .001 to a latitude is 3.6 seconds. That'd be something like 360 feet. That can put you quite a bit off.

If this method is used, a *very* small value needs to be added (anything on the order of 0E-4 is huge). More importantly, as practice had shown, there is no need to do anything like that. You can safely check for exact latitude equality and for cases where lat1 != lat2 use regular math with virtually infinitely small difference. Doubles and double-precision calculations have enough precision to put you within a fraction of an inch of a desired location if you simply let CPU do it's job.

You also do not need to switch to Great Circle when measuring distances E/W, there is a simple trigonometric formula for that.

[HamishB]

Originally Posted by brak

Could it be this bit from the DistanceBearingMercator?
Adding .001 to a latitude is 3.6 seconds. That'd be something like 360 feet. That can put you quite a bit off.

If this method is used, a *very* small value needs to be added (anything on the order of 0E-4 is huge).

1e-9 degree works out to be about 0.1mm (at worst), or about one order of magnitude better than the best RTK GPS will give you.

More importantly, as practice had shown, there is no need to do anything like that. You can safely check for exact latitude equality and for cases where lat1 != lat2 use regular math with virtually infinitely small difference. Doubles and double-precision calculations have enough precision to put you within a fraction of an inch of a desired location if you simply let CPU do it's job.

I run into code about once a month that is buggy due to not using
if( fabs(a-b) < epsilon ) for FP comparisons. It's a good idea.

What Every Computer Scientist Should Know About Floating-Point Arithmetic


Hamish

[brak]

Originally Posted by HamishB

1e-9 degree works out to be about 0.1mm (at worst), or about one order of magnitude better than the best RTK GPS will give you.

Oh, so degrees are stored * 10E6? Didn't know that. Well, never mind then.
They seemed to me to be treated as double values directly storing degrees (in which case 0.001 is 1E-3, isn't it?) (Edit: Ok, so right before that I see code like "lon0x += 360.", which to me would indicate that they are degrees stored as double).

I run into code about once a month that is buggy due to not using
if( fabs(a-b) < epsilon ) for FP comparisons. It's a good idea.

What Every Computer Scientist Should Know About Floating-Point Arithmetic

It is in some places, it's not in others. It comes down to knowing your domain, possible values and required precision. There isn't a single universal answer.

[HamishB]

Originally Posted by brak

Oh, so degrees are stored * 10E6? Didn't know that. Well, never mind then.

not that I know of, I assume they are just doubles AFAIK.

They seemed to me to be treated as double values directly storing degrees (in which case 0.001 is 1E-3, isn't it?)

1852m/naut.mi * 60nm/degree gives you 111km/degree at the equator.
so 1e-9 degree is roughly about a millimeter there.

It is in some places, it's not in others. It comes down to knowing your domain, possible values and required precision. There isn't a single universal answer.

yeah but variables end up getting extended and abused and writen out to files and data gets moved between compilers and CPU architectures and endian conversions, so IMO it's best to just always do it the paranoid way so you don't have to worry.

another possible reason is to protect against NaNs when it gets used by atan2() for the bearing calc.


Hamish

[brak]

I was talking about that 0.001 value - I am not quite sure which 1E-9 you mean. In any case, I decided to test that (sorry, I am somewhat obsessive with numeric precision ).

So here is what I see - if I look at east/north values returned by toSM_ECC(), they are correct until latitude difference gets to 0.001 (precisely that additional value). Once it gets to that value, it never goes below - not surprising, since latitude stays above 0.001 Which means that mercator north remains high, which in turn results in large errors when latitude difference is under 0.001 (this means bearings between 85 - 95 degrees, for example).

Dumping this term restores precision of Mercator bearing/distance calculation.
It's easy to test. Use DistanceBearingMercator() to calculate distance between
lat/lon 1:1 and 1.0001:1.001. Then comment out the " if(fabs(lat1 - lat0) < 1e-4) mlat0 += .001;" and try again. The distance returned will be about 165 meters in the first case (which is wrong) and 112 meters in the second (which is correct).

Edit: that means commenting out lines 797-799 of the current codebase.

[HamishB]

Originally Posted by brak

I was talking about that 0.001 value - I am not quite sure which 1E-9 you mean.

the 1e-9 is what the 0.001 should be.

presumably the idea is to nudge the position by 1mm. So 0.001 would come from other code/an assumption that the units of that variable are in meters or feet. But these are in degrees.

In any case, I decided to test that

excellent.

So here is what I see - if I look at east/north values returned by toSM_ECC(), they are correct until latitude difference gets to 0.001 (precisely that additional value). Once it gets to that value, it never goes below - not surprising, since latitude stays above 0.001 Which means that mercator north remains high, which in turn results in large errors when latitude difference is under 0.001 (this means bearings between 85 - 95 degrees, for example).
Dumping this term restores precision of Mercator bearing/distance calculation.

try this patch:

Index: src/georef.c
===================================================================
RCS file: /cvsroot/opencpn/opencpn/src/georef.c,v
retrieving revision 1.21
diff -u -r1.21 georef.c
--- src/georef.c        17 Dec 2009 02:50:09 -0000      1.21
+++ src/georef.c        22 Feb 2010 04:32:43 -0000
@@ -795,8 +795,8 @@
       //    This idea comes from Thomas(Cagney)
       //    We simply require the dlat to be (slightly) non-zero, and carry on.
       mlat0 = lat0;
-      if(fabs(lat1 - lat0) < 1e-4)
-            mlat0 += .001;
+      if(fabs(lat1 - lat0) < 1e-9)
+            mlat0 += 1e-9;
 
       toSM_ECC(lat1, lon1x, mlat0, lon0x, &east, &north);

also, I notice that the inverse flattening term for the WGS84 ellipsoid is wrong in the {to|from}SM_ECC() functions:

Index: src/georef.c
===================================================================
RCS file: /cvsroot/opencpn/opencpn/src/georef.c,v
retrieving revision 1.21
diff -u -r1.21 georef.c
--- src/georef.c        17 Dec 2009 02:50:09 -0000      1.21
+++ src/georef.c        22 Feb 2010 04:42:37 -0000
@@ -443,7 +443,7 @@
       double test;
       double ypy;
 
-      double f = 1.0 / 298.;    // flattening .003355
+      double f = 1.0 / WGSinvf;       // WGS84 ellipsoid flattening parameter
       double e2 = 2 * f - f * f;      // eccentricity^2  .006700
       double e = sqrt(e2);
 
@@ -498,7 +498,7 @@
       double z, s0, lon1;
       double falsen, t, xi, esf;
 
-      double f = 1.0 / 298.;    // flattening .003355
+      double f = 1.0 / WGSinvf;       // WGS84 ellipsoid flattening parameter
       double es = 2 * f - f * f;      // eccentricity^2  .006700
       double e = sqrt(es);

regards,
Hamish

[HamishB]

Originally Posted by HamishB

yeah but variables end up getting extended and abused and writen out to files and data gets moved between compilers and CPU architectures and endian conversions, so IMO it's best to just always do it the paranoid way so you don't have to worry.

... after thinking about it a bit more, I think it's fairly safe to assume that after conversion to radians and back, and perhaps to Mercator and back, or a great-circle distance/bearing calc with distance=0, that FP precision will not be preserved to the least significant bit.


Hamish

[brak]

Originally Posted by HamishB

the 1e-9 is what the 0.001 should be.

try this patch:

I am sure it will work and if/when that makes it to the release, the OP will be happy too.

[Psyches]

Hamish, brak: georef.c fixed per Hamish's patch - checked in as v1.22;

Any issues, someone please let me know. Also, I haven't had much time to scan posts...any other such patches (e.g. clear bugs with low risk fixes), feel free to email me as well, to catch my attention.

Mark

[cagney]

The first Mac release candidate for 1.3.6 is now available on opencpn.org here: http://opencpn.org/sites/default/fil...3.6-mac-rc.zip

Thank you "Jallum", for providing the package!

All Mac users ....please test this pre-release.



Thomas

[CarinaPDX]

The 1.3.6 RC has a couple of very obvious issues when I execute it on 10.6.2, which I would like to report or at least check to see if they are being tracked already. Where are Mac-specific bugs being tracked? For that matter, is the bug tracker on Sourceforge the place for general bugs?

There are two or three threads that are currently dealing with Mac development, which is confusing to say the least. I suspect that one of the problems I have noted was already fixed by a patch Erwin developed and referred to on a different thread. Wouldn't it be best to have a single, sticky thread dedicated to Mac development?

I would like to be a beta tester, if nothing else, for the Mac versions. Is there a mailing list to alert testers to new versions needing testing?

I'm not being critical here - I just want to help move this great project forward.

Greg