gdal grib driver gets the origin and corner coords wrong

[winkey]

rush@winkey ~ $ gdalinfo gfs.t06z.pgrb2f00
Driver: GRIB/GRIdded Binary (.grb)
Files: gfs.t06z.pgrb2f00
Size is 720, 361
Coordinate System is:
GEOGCS["Coordinate System imported from GRIB file",
    DATUM["unknown",
        SPHEROID["Sphere",6371229,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (0.000000000000000,270.000000000000000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Corner Coordinates:
Upper Left  (       0.000,     270.000) (  0d 0'0.01"E,270d 0'0.00"N)
Lower Left  (   0.0000000,  89.5000000) (  0d 0'0.01"E, 89d30'0.00"N)
Upper Right (     360.000,     270.000) (360d 0'0.00"E,270d 0'0.00"N)
Lower Right (     360.000,      89.500) (360d 0'0.00"E, 89d30'0.00"N)
Center      (     180.000,     179.750) (180d 0'0.00"E,179d45'0.00"N)

degrib gfs.t06z.pgrb2f00 -out stdout -C -Csv | less
GDS | Number of Points | 259920
GDS | Projection Type | 0 (Latitude/Longitude)
GDS | Shape of Earth | sphere
GDS | Radius | 6371.229000 (km)
GDS | Nx (Number of points on parallel) | 720
GDS | Ny (Number of points on meridian) | 361
GDS | Lat1 | 90.000000
GDS | Lon1 | 0.000000
GDS | u/v vectors relative to | easterly/northerly
GDS | Lat2 | -90.000000
GDS | Lon2 | 359.500000
GDS | Dx | 0.500000 (degrees)
GDS | Dy | 0.500000 (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

the dataset can be found at http://atmos.ucsd.edu/gfs.t06z.pgrb2f00
i did not attach it because of its size
-rw-r--r-- 1 rush users 42806506 May 19 02:28 gfs.t06z.pgrb2f00

[winkey]

sorry here the link again ​http://atmos.ucsd.edu/gfs.t06z.pgrb2f00

[Even Rouault]

Reproduced with GDAL 1.6.0dev

[warmerdam]

I believe this may be closely related to a problem privately reported by Erick Jones, so I'll include his analysis here:

Frank,

We're having trouble displaying GRIB2 data from GFS using GDAL. What follows is everything I've been able to determine about this problem. I apologize for the lengthy email; I wasn't sure if I should send this to you directly, or if I should submit it as a bug report (which I'm not entirely certain how to do).

The image is correct, but the position is incorrect: it's supposed to be defined for latitudes 90 S to 90 N (the entire globe), but GDAL thinks it ranges from 90 N to 270 N, which of course makes no sense. The file itself does not appear to be corrupt. If we use "wgrib2" or "degrib", the results are correct, but any GDAL based tool (including gdal_translate and fwtools) interprets it incorrectly.

I just checked out the most recent version of GDAL from subversion and the problem is still there.

I have attached a sample GRIB2 file exhibiting this problem, gfstmp.grib2. Here is the output when I run gdalinfo against this file:

Driver: GRIB/GRIdded Binary (.grb)
Files: gfstmp.grib2
Size is 360, 181
Coordinate System is:
GEOGCS["Coordinate System imported from GRIB file",
    DATUM["unknown",
        SPHEROID["Sphere",6371229,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
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)
Band 1 Block=360x1 Type=Float64, ColorInterp=Undefined
  Description = 50000[Pa] ISBL="Isobaric surface"
  Metadata:
    GRIB_UNIT=[K]
    GRIB_COMMENT=Temperature [K]
    GRIB_ELEMENT=TMP
    GRIB_SHORT_NAME=50000-ISBL
    GRIB_REF_TIME=  1222257600 sec UTC
    GRIB_VALID_TIME=  1222279200 sec UTC
    GRIB_FORECAST_SECONDS=21600 sec
    GRIB_PDS_PDTN=0
    GRIB_PDS_TEMPLATE_NUMBERS=0 0 2 0 96 0 0 0 1 0 0 0 6 100 0 0 0 195 80 255 0 0 0 0 0 

As you can see, the "upper left" and "upper right" specify 270 as the Y coordinate.

I believe I have an idea of what's causing this. First of all, here is the code that seems to be the immediate problem, starting at line 580 in gribdataset.cpp:

    else
    {
        rMinX = meta->gds.lon1; // longitude in degrees, to be transformed to meters (or degrees in case of latlon)
        rMaxY = meta->gds.lat1; // latitude in degrees, to be transformed to meters

        if (meta->gds.scan == GRIB2BIT_2) // Y is minY, GDAL wants maxY
            rMaxY += (meta->gds.Ny - 1) * meta->gds.Dy; // -1 because we GDAL needs the coordinates of the centre of the pixel
        rPixelSizeX = meta->gds.Dx;
        rPixelSizeY = meta->gds.Dy;
    } 

Since meta->gds.scan is in fact equal to GRIB2BIT_2 (0x40), the addition takes place. However, contrary to expectation, rMaxY was already 90, so after this executes it will become 270!

What doesn't make sense here is that meta->gds.lat1 is the "starting" latitude, and according to the scanning mode, latitude should increase from there -- see Table 3.4 of the Grib documentation, which you can find at the following URL:

​http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table3-4.shtml

In particular, bit 2 is set, so rows scan in the +y direction.

So if these are the values specified in the GRIB file, then it would seem to be corrupt data -- but then, how does "degrib" or "wgrib2" end up showing correct values?

Actually, these are NOT the values present in the data. Something in the GRIB2 API is mangling them. See attached image "grib_hex.png", where I did a hex dump of this data file and looked at the values present. The image is marked up to show where starting latitude, ending latitude, and scan mode are specified. In particular, note that the scan mode is 0, not 0x40. (Scan mode 0 would definitely make more sense with a starting latitude of 90 N and ending latitude of 90 S, since it means that rows scan in the -y direction.)

The mangling occurs in the TransferFloat function in grib2api.c. Here is the code starting at line 515:

      } else {
         for (i = 0; i < ngrdpts; i++) {
            ScanIndex2XY (i, &x, &y, *scan, nx, ny);
            /* ScanIndex returns value as if scan was 0100(0000) */
            curIndex = (x - 1) + (y - 1) * nx;
            myAssert (curIndex < nd2x3);
            ain[curIndex] = fld[i];
         }
      }
      *scan = 64 + (*scan & 0x0f); 

This alters the scan mode to be equal to 64 (0x40), and reorders the rows to match that scan mode. Indeed, the comment above the function indicates that this is the intended behavior.

Unfortunately, what this code does not do is reverse the starting and ending latitudes. Thus, we end up with an inconsistent grid description.

I'm not really sure what would be the best way to correct this. We can't skip the problematic test in gribdataset.cpp, because it would be needed if the original scan mode were in fact 0x40. Perhaps we need to just look at lat1 and lat2, and store whichever is higher in rMaxY? That seems like a bit of a messy hack, though. Perhaps we need to examine how degrib handles this; after all, it uses the same grib2api.c module, so should be affected by this same issue, but it still gives correct results.

[hamish]

Hi,

I too am seeing a shift in GRIB data.

AFAICT it is simply shifted half a cell to the SE from where it should be. I assumed this was due to the usual grid vs. cell-center convention mistake, but that's just a round up the usual suspects guess.

attaching screenshot.

list of commands to recreate my image is here:

​http://grass.osgeo.org/wiki/GRIB

GRIB data in the screenshot downloaded from globalmarinenet.net.

thanks, Hamish

ps- Is there an internal tag for no-data? In the above data it is represented by "9999", which is easy to set with gdal_translate -a_nodata, but in other GRIB data I have it is what I guess is the max for Float64, which is a number 4-6 lines long; too much to type by hand.

[hamish]

grib cells shifted 1/2 to SE

[hamish]

grib cells shifted 1/2 to SE seem to be a different issue, so moved into its own ticket (#2637)

Hamish

[winkey]

The last couple of days i tried a few things including a newer g2clib and studying diffs between gdals degrib code and degrib, this is what i ended up with.

If you get rid of that check around line 580 in gribdataset.cpp:

if (meta->gds.scan == GRIB2BIT_2) // Y is minY, GDAL wants maxY
            rMaxY += (meta->gds.Ny - 1) * meta->gds.Dy; // -1 because we GDAL needs the coordinates of the centre of the pixel

everything seems to work on the data i tried ​http://www.ftp.ncep.noaa.gov/data/nccf/com/gfs/prod/ both sst and gfs ​http://www.ftp.ncep.noaa.gov/data/nccf/com/wave/

Does anyone have any cases where the original code worked properly?

Brian

[winkey]

that also works with ecmwf from the ucar ldm feed.

however it does not work with stuff i downloaded from globalmarinenet.net (hamish's data)

as a side note some data seems to need a half pixel shift, some does not. i think this may be a product dependent issue.

Brian

[grazianogiuliani]

+1 for the modification to the code. I have tested it with GFS, ECMWF, NAM, RUC, DWD, grib_api produced and w3lib produced GRIB files. In all cases simply commenting the two lines above does the trick, but I have added the code used by degrib itself in mymapf.c: it just checks for the image to be oriented consistently with the order forced in grib2api.c. About the pixel shift, I have tested the patch proposed in #2637, and built a comprehensive patch to solve that problem, this problem, and incorporate latest versions of degrib (1.94) and g2_clib (1.1.9). I have not any Microsoft OS to test it, and I can grant it working only on Linux x86_64.

[grazianogiuliani]

Omnibus patch to solve #2550, #2637 and upgrade grib driver.

[hamish]

Replying to winkey:

that also works with ecmwf from the ucar ldm feed.

however it does not work with stuff i downloaded from globalmarinenet.net (hamish's data)

as a side note some data seems to need a half pixel shift, some does not. i think this may be a product dependent issue.

is there some authoritative & well documented dataset to test with that we can blindly trust? globalmarinenet.net seems to do some 3rd party reprocessing so obviously not the most trustworthy source.

some e.g. NOAA source produced directly from MM5 model output would be ideal. is this it: ?

Brian wrote:

everything seems to work on the data i tried ​http://www.ftp.ncep.noaa.gov/data/nccf/com/gfs/prod/ both sst and gfs ​http://www.ftp.ncep.noaa.gov/data/nccf/com/wave/

ie are we understanding & solving the root of the problem with this patch, or just throwing in hacks until it seems to work?

thanks, Hamish

[grazianogiuliani]

We have in the case of the GRIB data format the following list of problems:

• The format is WMO standard but it does NOT exist any reference implementation
• The format is built by accretion of modifications and rely on external tables which are compiled by meteo centers and it does not exist any central repository for them.
• It is used to exchange mainly meteo and model data products, and it has been extended for satellite data, mainly geostationary as Meteosat
• Lot of codecs exist in the wild, but many of them are written in Fortran language and are made just for internal use in the data originating center.

Taken this, the two major organizations which distribute GRIB coded data are the NCEP in the US and the ECMWF/Eumetsat in Europe (sorry for the others, but I am talking about global model datasets, not just local ones). Each of them has its reference implementation of the codec, its own internal tables describing data, etc.

• grib_api for ECMWF
• w3lib for GRIB1 and g2lib for GRIB2 for NCEP

National Weather Service in the US distributes also two programs as swiss army knives, wgrib for GRIB1 and wgrib2 for GRIB2, able to easily decode, inventory and dump and for GRIB2 also transform in NetDF. Any of this implementation does not completely support all the documented features of the format, but just the subset which is commonly used: WMO documents are not revised frequently and they just add features, never remove or obsolete anything. The degrib code used in GDAL comes from NOAA NWS Meteorological Development Laboratory, and its aim is to access data stored in the National Digital Forecast Database, so it does support a subset of the format, with just the parts needed to access those data. So, we are talking of a partial implementation which partially supports a subset of the specifications, distributed by one of the data providers. To answer Hamish question in short:

it does not exist an authoritative & well documented dataset and GDAL implementation is not supposed to successfully decode any GRIB message.

Then comes my opinion, in which if the driver is able to access global datasets from NCEP, NOAA, ECMWF, is able to georeference local model data at least of the NAM (actually WRF model GRIB1 encoded data produced by w3lib), we can be happy with this.

About the patch, it just implements what is in another part of the degrib codebase, so it is not "the solution", but "a solution", which is consistent with GDAL incorporated code. And Yes, it implements what is strictly speaking the "a bit of a messy hack" Erick Jones speaks about, but it is taken from the same overall messy hacked codebase of degrib. For the bump in the degrib code proposed, it comes from the fact that both NCEP and ECMWF have made modification to their tables, so we have less "unknown parameter with unkonwn units" in the inventory, and a non checked free have been corrected.

As a conclusion, we are using one of the implementations, in my opinion not one of the most complete ones, and the whole driver is a complete hack. But it is enough for lot of us.

[warmerdam]

I'm not really in a position to review this big patch (4.5MB uncompressed!) and I'm concerned that some custom changes to the degrib library may have been lost in the process.

I'm not sure what to do next, but I'd be interested in feedback from other grib users on the patch.

[hamish]

Replying to warmerdam:

I'm not really in a position to review this big patch (4.5MB uncompressed!) and I'm concerned that some custom changes to the degrib library may have been lost in the process.

I'm not sure what to do next, but I'd be interested in feedback from other grib users on the patch.

step 1: the patch should be split into smaller chunks & reupload: one for each bugfix and one for the new degib library version.

Hamish

[ejones]

GFS grib file for which GDAL exhibits this problem

[ejones]

Binary examination of attached grib file, as mentioned in my writeup

[ejones]

Proposed patch -- just let rMaxY equal larger of (lat1, lat2)

[ejones]

The latest version of GDAL (1.9.0) still seems to have the same problem described in this ticket.

In our display software, we have been using the change shown in "simple_patch.txt" for years now, and it allows GFS data to show up properly. I'd like to revisit the idea of getting this added to the official code base, so we don't have to keep manually patching new versions of GDAL.

I definitely don't think this calls for updating to a whole new version of degrib, along with all the new problems that could introduce. The advantage of this patch is that it clearly only affects the one specific code path that was problematic, so there should not be any side effects.

[warmerdam]

On review I think it is appropriate to go ahead and apply the patch.

[Daniel Morissette]

Committed simple_patch.txt in SVN trunk r23122.

[hamish]

any idea how this affects #2637 (or not)?

thanks, Hamish

[Daniel Morissette]

Unfortunately I don't know about #2637. We'd need someone more familiar with the format to answer.

[winkey]

I can confirm this fixes my issues with the global gfs .5 deg grib files. there used to be a gap at the prime meridian

[wget ​http://www.ftp.ncep.noaa.gov/data/nccf/com/gfs/prod/gfs.2011092906/gfs.t06z.pgrb2bf00]

Band 183 Block=720x1 Type=Float64, ColorInterp=Undefined
  Description = 87500[Pa] ISBL="Isobaric surface"
  Metadata:
    GRIB_COMMENT=Temperature [K]
    GRIB_ELEMENT=TMP
    GRIB_FORECAST_SECONDS=0 sec
    GRIB_REF_TIME=  1317276000 sec UTC
    GRIB_SHORT_NAME=87500-ISBL
    GRIB_UNIT=[K]
    GRIB_VALID_TIME=  1317276000 sec UTC

gdal_contour gfs.t06z.pgrb2bf00 -b 183 -f libkml -i 1.5 gribtest.kml
0...10...20...30...40...50...60...70...80...90...100 - done.

Thanks

[winkey]

I feel i should note that if the degrib and g2clib portions of the driver are never updated. gdal will never learn to handle new grib files

Brian

Replying to ejones:

The latest version of GDAL (1.9.0) still seems to have the same problem described in this ticket.

In our display software, we have been using the change shown in "simple_patch.txt" for years now, and it allows GFS data to show up properly. I'd like to revisit the idea of getting this added to the official code base, so we don't have to keep manually patching new versions of GDAL.

I definitely don't think this calls for updating to a whole new version of degrib, along with all the new problems that could introduce. The advantage of this patch is that it clearly only affects the one specific code path that was problematic, so there should not be any side effects.