grib1 & grib2 : pixel size precision introduces error for corner coordinates
|Reported by:||mprecise||Owned by:||warmerdam|
|Severity:||normal||Keywords:||grib, grib1, grib2|
|Cc:||mprecise, Even Rouault, aboudreault, ejones, bcrosby, Daniel Morissette|
This problem seems to exist in both grib1 and grib2 files, but according to the files I've tested with, the error due to grib1's limited coordinate precision has the potential to lead to a far greater error.
Grib1 documentation: http://www.nco.ncep.noaa.gov/pmb/docs/on388/
Grib2 documentation: http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc.shtml
From my interpretation of grib1 and grib2 documentation, it seems that two corner coordinates are always given in grib files (Lat1/Lon1 and Lat2/Lon2) and these are not optional. However, GDAL doesn't use or store the Lat2/Lon2 values in any way. When working with a GeoTransform, GDAL calculates this corner using Lat2 = Lat1 + (num_rows * pixel_height), and Lon2 = Lon1 + (num_cols * pixel_width). This can lead to precision problems that varies based on the precision of the pixel width (or pixel height) and the number of rows (or columns).
I'm attaching a grib1 file, eng5min.20111006, that demonstrates the problem. Additional versions of this same test data can be found at ftp://polar.ncep.noaa.gov/pub/cdas, anonymous login. I used degrib 1.96 (2011-02-23) as a comparison - the output of that is below. I've also inspected the file manually - a hex dump image is attached with highlighted areas, eng5min.20111006_analysis.png. What I saw by manually inspecting the file matches up with what degrib reports, so degrib is not in error.
degrib output: GDS | Number of Points | 9331200 GDS | Projection Type | 0 (Latitude/Longitude) GDS | Shape of Earth | sphere GDS | Radius | 6371.200000 (km) GDS | Nx (Number of points on parallel) | 4320 GDS | Ny (Number of points on meridian) | 2160 GDS | Lat1 | 89.958000 GDS | Lon1 | 0.042000 GDS | u/v vectors relative to | easterly/northerly GDS | Lat2 | -89.958000 GDS | Lon2 | -0.042000 GDS | Dx | 0.083000 (degrees) GDS | Dy | 0.083000 (degrees) GDS | Input GRIB2 grid, scan mode | 0 (0000) GDS | Output grid, scan mode | 64 (0100) GDS | (.flt file grid), scan mode | 0 (0000) GDS | Output grid, scan i/x direction | positive GDS | Output grid, scan j/y direction | positive GDS | (.flt file grid), scan j/y direction | negative GDS | Output grid, consecutive points in | i/x direction GDS | Output grid, adjacent rows scan in | same direction
Compare degrib's Lat2 and Lon2 values to the Lower Right corner coordinate reported by gdalinfo (ignoring the fact that GDAL uses the range 0-360 for longitude, that is not relevant to this issue):
Origin = (0.042000000000000,89.957999999999998) Pixel Size = (0.083000000000000,-0.083000000000000) Corner Coordinates: Upper Left ( 0.0420000, 89.9580000) ( 0d 2'31.20"E, 89d57'28.80"N) Lower Left ( 0.0420000, -89.3220000) ( 0d 2'31.20"E, 89d19'19.20"S) Upper Right ( 358.602, 89.958) (358d36' 7.20"E, 89d57'28.80"N) Lower Right ( 358.602, -89.322) (358d36' 7.20"E, 89d19'19.20"S) Center ( 179.3220000, 0.3180000) (179d19'19.20"E, 0d19' 4.80"N)
Calculations GDAL uses to get the Lower Right corner coordinates shown above:
latitude: -0.083 * 2160 = -179.28 89.958 + (-179.28) = -89.322 longitude: 0.083 * 4320 = 358.56 0.042 + 358.56 = 358.602
Notice that latitude is off by 0.636 degrees, and longitude is off by 1.356 degrees vs the Lat2 and Lon2 reported by degrib.
From inspecting the grib1 GDS section 2, it looks like all latitudes and longitudes in a lat/lon grid are limited to 3 decimal places of precision. If we had instead done this:
pixel_height = (Lat1 - Lat2) / Ny => (89.958 - (-89.958)) / 2160 = 0.083294444444444444444444444444444 pixel_width = (360 - (Lon1 - Lon2)) / Nx => (360 - (0.042 - (-0.042))) / 4320 = 0.083313888888888888888888888888889
we see that each pixel's size can more accurately be represented. Lat1/Lat2/Lon1/Lon2 are still limited by 3 decimal precision, but the error is mitigated by the number of columns and rows.
The same problem exists in grib2 files, but to a lesser extent because grib2 uses 6 decimal places for latitudes and longitudes. In the grib2 file I looked at, Lon2 is 310.486386 as reported by degrib, but 310.4865 by GDAL. Similarly, I have a sample file and analysis of it, too, if needed.
This problem could potentially exist for other file formats besides grib1 and grib2, but I'm not familiar with those.
I'm attaching a patch that is a potential solution to this problem. The patch relies on a recent fix in trunk related to #2550. I'm unsure if the way I handle the sign change between Lon1 = 0.042, Lon2 = -0.042 is good enough for a general solution for all grib1 and grib2 files.
One consequence of this patch is that pixel size returned by all other uses of GeoTransform will no longer be what is stored in the file, but in my mind that is the lesser evil.