[Gdal-dev] precision pitfall

[Ray Gardener]

I don't know if this has been discussed earlier, but I thought I'd 
mention it as I've come across it. It's not a GDAL error (although it 
might explain some things people might have come across, so it might be 
good to check GDAL):

Files that include georef data may not have enough precision to prevent 
subtle coordinate errors downstream unless workarounds are used.

Example: GTOPO30 .HDR files include XDIM/YDIM members that are typically 
listed as 0.00833333333333 (1/4800th of 40 degrees).

Doing the following math with double-precision numbers on x86:

     double xdim = atof("0.00833333333333");
     double x = 20; // dem's west edge, 20 deg east.
     double width = xdim * 4800;
     x += width;

     printf(x);        // Reports 60.00000, and even
     printf(x * 3600); // arcseconds will be reported as 216000,

                   // but...

     printf("%d deg, %d', %d\"",
        (int)d,
        ((int)(d * 60)) % 60,
        ((int)(d * 3600)) % 60); // We get: 59 deg, 59' , 59"


Using (int)x or floor(x) to get the whole-number degrees will be 59. 
Turns out there are thirteen double-precision numbers that printf() (and 
the VC debugger) will report as 60 that are slightly less than 60.

If the GTOPO30 HDR file had four more digits of precision, the error 
goes away. Epsilon math will fix this, but just to show where a 
contributing source of the error can be.

In a strict sense, there's actually no error if we take the HDR file at 
face value and treat the extent width as being a hair under 40 degrees. 
But we know this isn't the DEM's width, it is intended to be 40. Being 
strict would cause annoying off-by-one-arcsrcond coordinate readouts to 
appear, and may cause worse errors if such values are used in 
computations downstream (or data exchange roundtripping). A better way 
to fix this might have been for the file spec to say how wide the DEM 
is, or where the extents are absolutely located, instead of how far 
apart its gridposts are.

Ray