Opened 5 years ago

Closed 5 years ago

#4287 closed defect (fixed)

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

Reported by: mprecise Owned by: warmerdam
Priority: normal Milestone:
Component: default Version: 1.8.1
Severity: normal Keywords: grib, grib1, grib2
Cc: mprecise, Even Rouault, aboudreault, ejones, bcrosby, dmorissette

Description

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.

Attachments (4)

eng5min.20111006_analysis.png (21.3 KB) - added by mprecise 5 years ago.
Hex dump analysis
eng5min.20111006.gz (357.6 KB) - added by mprecise 5 years ago.
sample grib1 file exhibiting the error
pixelsize.patch (923 bytes) - added by mprecise 5 years ago.
potential fix for the error
pixelsize2.patch (1.2 KB) - added by mprecise 5 years ago.
A patch to fix an error in pixelsize.patch.

Download all attachments as: .zip

Change History (14)

Changed 5 years ago by mprecise

Hex dump analysis

Changed 5 years ago by mprecise

Attachment: eng5min.20111006.gz added

sample grib1 file exhibiting the error

Changed 5 years ago by mprecise

Attachment: pixelsize.patch added

potential fix for the error

comment:1 Changed 5 years ago by aboudreault

Cc: Even Rouault aboudreault added

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?

comment:2 in reply to:  1 Changed 5 years ago by 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.

comment:3 Changed 5 years ago by 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.

comment:4 Changed 5 years ago by aboudreault

Cc: ejones bcrosby added

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.

comment:5 Changed 5 years ago by 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.

comment:6 Changed 5 years ago by 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.

comment:7 Changed 5 years ago by ejones

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

comment:8 Changed 5 years ago by 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.

Changed 5 years ago by mprecise

Attachment: pixelsize2.patch added

A patch to fix an error in pixelsize.patch.

comment:9 Changed 5 years ago by dmorissette

Cc: dmorissette added

comment:10 Changed 5 years ago by aboudreault

Resolution: fixed
Status: newclosed

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

Note: See TracTickets for help on using tickets.