cannot retrieve EGM96 heights for coordinates very near cell boundaries

[ravenAtSafe]

Coordinates sufficiently close to EGM96 cell boundaries result in the wrong cell being chosen, and later validation aborts because of that.

Here are two sample points that cause issues:

lng=0x40467fffffffffff/44.99999999999999 + any lat
lat=0x4046800000000001/45.00000000000001 + any lng

The following is the hack I am using to work around this:

Index: C:/svn/trunk1-x86/3rd/redistributable/free_use/csmap/Source/CS_egm96.c
===================================================================
--- C:/svn/trunk1-x86/3rd/redistributable/free_use/csmap/Source/CS_egm96.c	(revision 78538)
+++ C:/svn/trunk1-x86/3rd/redistributable/free_use/csmap/Source/CS_egm96.c	(revision 78539)
@@ -396,6 +396,32 @@
 	tt = (lclLng - cellNW [LNG]) / __This->density [LNG];
 	uu = (cellNW [LAT] - lclLat) / __This->density [LAT];
 
+   // Safe Software Inc.
+   // PR:31774
+   // SECTION_START: Support coordinates very near cell boundaries
+   //
+   // The 1.0E-12 padding above will have pushed us into the wrong cell for
+   // coordinates very close to cell boundaries.  Examples:
+   //    lng=0x40467fffffffffff/44.99999999999999 + any lat
+   //    lat=0x4046800000000001/45.00000000000001 + any lng
+   //
+   //    http://babbage.cs.qc.edu/IEEE-754/Decimal.html is helpful
+   //
+   // Actual impact of being in the wrong cell is likely extremely small (thanks
+   // to the weighting), but since we will fail the following validation we might
+   // as well massage tt and uu into valid ranges as that will give very slightly
+   // better results than just disabling the validation.  Hopefully there isn't
+   // any way to get truly bad values here because it isn't obvious to me what a
+   // reasonable epsilon is here.  And of course it would be great if we had just
+   // chosen the correct cell to begin with, but potential floating point
+   // shenanigans will have to be considered far more thoroughly than I am
+   // willing to right now.
+   if (tt < 0.0) tt = 0.0;
+   if (tt > 1.0) tt = 1.0;
+   if (uu < 0.0) uu = 0.0;
+   if (uu > 1.0) uu = 1.0;
+   // SECTION_END: Support coordinates very near cell boundaries
+
 	/* Again, some defensive stuff. */
 	if (tt < 0.0 || tt > 1.0 || uu < 0.0 || uu > 1.0)
 	{

[Norm Olsen]

Sorry for the delay in responding. I did indeed see this one some months ago,m but couldn't deal with it at the time, and then forgot all about it. I need to check the Trac more regularly.

I believe the true solution to this problem is what the Canadians use in their National Transformation. That is, detect when the point to be caluclated is exactly (within a certain range) on a grid cell line, and then switch the calculatiuon method to one that does a linear interpolation bwteen the two appropriate points which define the grid cell line.

In this case, the tholerance would be about 5.0E-12. The precision of IEEE real numbers is 14+ digits. Since we're dealing with numbers which routinely have 3 digits to the left of the decimal point, 5.0E-12 should do it. In rough numbers, 5.0E-012 degrees equals approximately 0.001 millimeters, so that value is certainly small enough.

This is obviously not a trivial fix, so it will have to wait a couple of weeks.

[Norm Olsen]

Att version 2460, this was corrected bey removing the additional constant added to the real portion of the eleNbr and recNbr values. The additional constant is appropriate in the case where the grid cell density is not a nice number. For example, a grid cell density of 5 minutes would produce a grid cell size of 0.016777777777777 degrees. Not a nice number.

In the case of EGN96, the grid cell denisty is 15 minutes, which in degrees is 0.25000000 (a very nice number), and the extra constant is not necessary. We make this change with the fact that encountering geoid height data in this format with a grid cell density other than 15 minutes is very unlikely.

Submission produced versionnnn 2460.