grib1 & grib2 : pixel size precision introduces error for corner coordinates

[mprecise]

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.

[mprecise]

Hex dump analysis

[mprecise]

sample grib1 file exhibiting the error

[mprecise]

potential fix for the error

[aboudreault]

I don't know grib format at all... so it's a little hard for me to agree on anything. After applied the patch, I have this output:

Origin = (0.042000000000000,89.957999999999998)
Pixel Size = (0.083313888888889,-0.083294444444444)
Corner Coordinates:
Upper Left  (   0.0420000,  89.9580000) (  0d 2'31.20"E, 89d57'28.80"N)
Lower Left  (   0.0420000, -89.9580000) (  0d 2'31.20"E, 89d57'28.80"S)
Upper Right (     359.958,      89.958) (359d57'28.80"E, 89d57'28.80"N)
Lower Right (     359.958,     -89.958) (359d57'28.80"E, 89d57'28.80"S)
Center      ( 180.0000000,   0.0000000) (180d 0' 0.00"E,  0d 0' 0.01"N)

I assume it to be correct, according to the GRIB output. The patch passed the gdal autotest. FrankW, rouault, what do you think of this bug and fix? Also, do you see any important issue regarding the consequence?

[mprecise]

Alan, I'm confirming that what you posted is the output I expect for the sample file. As you mentioned, I'm also unsure of any consequences.

[Even Rouault]

I don't know the GRIB world, so I have no authoritative opinion on this. The idea behind the patch looks reasonable. Apart from testing on a larger set of data (or asking interested testers to help doing it) to check how well it works, difficult to predict the effects.

Alan, I've not verified what the autotest actually tests on this and don't know how representative it is. Anyway, in the case, it is just a matter of checking the output of gdalinfo before/after the patch.

[aboudreault]

Alright. I do not see any bad effects neither. I've committed the patch in r23236 so it can be tested using trunk. bcrosby, ejones ... since you seem to know the format GRIB better than us. Feel free to test this with your data and comment.

[bcrosby]

I've been able to confirm this bug with my data (the lower right coordinates are off by 1 degree, since my pixel size is 1.0 degrees.

Output from degrib:

GDS | Number of Points | 65160
GDS | Projection Type | 0 (Latitude/Longitude)
GDS | Shape of Earth | sphere
GDS | Radius | 6371.229000 (km)
GDS | Nx (Number of points on parallel) | 360
GDS | Ny (Number of points on meridian) | 181
GDS | Lat1 | 90.000000
GDS | Lon1 | 0.000000
GDS | u/v vectors relative to | easterly/northerly
GDS | Lat2 | -90.000000
GDS | Lon2 | 359.000000
GDS | Dx | 1.000000 (degrees)
GDS | Dy | 1.000000 (degrees)
GDS | Input GRIB2 grid, scan mode | 64 (0100)
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

gdalinfo:

Origin = (0.000000000000000,270.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Corner Coordinates:
Upper Left  (       0.000,     270.000) (  0d 0'0.01"E,270d 0'0.00"N)
Lower Left  (   0.0000000,  89.0000000) (  0d 0'0.01"E, 89d 0'0.00"N)
Upper Right (     360.000,     270.000) (360d 0'0.00"E,270d 0'0.00"N)
Lower Right (     360.000,      89.000) (360d 0'0.00"E, 89d 0'0.00"N)
Center      (     180.000,     179.500) (180d 0'0.00"E,179d30'0.00"N)

I haven't test the patched version of gdal on my data, but looking at the diff, this doesn't seem like it would pose a problem for me.

[mprecise]

bcrosby, I think your data shows a different problem, as it doesn't seem to be related to precision in the value of pixel size since you're using a nice round number dx=1, dy =-1.

Your degrib and gdalinfo output did lead me to look more closely at my patch. We're still investigating, but there may be something else going on in the grib driver besides the pixel precision I reported.

[ejones]

Yeah, I think bcrosby's example is actually an instance of ticket #2637.

[mprecise]

I have an amendment to my previous patch. brosby's degrib output with nice round numbers made me realize that I should have been subtracting 1 from Dx and Dy in the calculations used in the other patch.

90 to -90 is 181 unique points, including point 0 through the equator

0 to 359 is 360 unique points

My first patch would have given this for bcrosby's file:

(max_lat - min_lat) / Ny => (90 - (-90)) / 181 = 0.9944751381215
(lon2 - lon1) / Nx => (359 - 0) / 360 = 0.99722222222

But by subtracting 1 from Ny and Nx, we get the pixel sizes as reported by degrib:

(max_lat - min_lat) / (Ny - 1) => (90 - (-90)) / (181 - 1) = 1.0
(lon2 - lon1) / (Nx - 1) => (359 - 0) / (360 - 1) = 1.0

The 2nd patch I'm attaching corrects this defect. I've also added some sanity checks to attempt to handle any GRIB files that the two combined patches shouldn't be applied to. If a sanity check fails, the corresponding pixel width or pixel height will get its value from the GRIB file, as it did before.

I've taken the patch attached to #2637 in consideration with my changes and have assumed it will be applied. They should function together without issue.

[mprecise]

A patch to fix an error in pixelsize.patch.

[aboudreault]

Thanks for the patch correction mprecise! Fixed and committed in r23302.