grib driver: grid vs cell-center convention wrong?

[hamish]

Hi,

I'd previously sent this to ticket #2550 but as it seems to be mostly unrelated to the problem described there I think it deserves its own one.

I 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.

a list of commands to recreate my image is here:

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

GRIB data in the screenshot downloaded from globalmarinenet.net. I can't guarantee the integrity of globalmarinenet's data (I have no idea about it other than the link); so it would be good to verify that this is happening with other datasets (ie something directly from the WMO) and in the spec. In the PDF which speaks of triangulation method it talks about data being at the vertices, fwtw.

thanks, Hamish

ps- Is there an internal tag in the format 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]

I'd note that the cell size for these files are typically 1-1.25 degrees, so the data is offset > 50km.

Hamish

[tbouffon]

Hi, I haven't had time to check 1.6.2, but since the bug is not marked as closed, the attached patch seems to shift the area back. I only tested it for one day, so use with care !

[tbouffon]

Sorry I didn't mean to close the bug, just wanted to says my patch works for me. By the way, The bug is still in 1.6.2, and the patch still works fine

[hamish]

updated patch vs. latest svn trunk for grid registration error

[hamish]

recut patch compatible with the latest svn trunk added.

apparent confirmation of the grid-registration error diagnosis:

​http://thread.gmane.org/gmane.comp.gis.grass.user/35145/focus=35661 ​http://thread.gmane.org/gmane.comp.gis.grass.user/35145/focus=35666

Hamish

[hamish]

sample data in #3629 (a duplicate)

[rhijmans]

Replying to hamish:

sample data in #3629 (a duplicate)

[hamish]

It would be really nice if this error gets fixed. It has been open for a long time.

hey you edited the wrong box! :)

(easy to do, I filed a bug upstream at trac HQ about a year ago to make it fold away into a little [-] expanding [+] box. maybe just open that for the submitter or sysops to have access to change?)

description restored..

[Kyle Shannon]

Replying to hamish:

It would be really nice if this error gets fixed. It has been open for a long time.

hey you edited the wrong box! :)

(easy to do, I filed a bug upstream at trac HQ about a year ago to make it fold away into a little [-] expanding [+] box. maybe just open that for the submitter or sysops to have access to change?)

description restored..

Have you tried version 1.8? According to your link, you are just using gdal_translate to create tif files. RFC 33: GTiff - Correcting PixelIsPoint Interpretation (Adopted) fixed how GTiff handles the origin registration. On my system after a translate, I get the same coordinates for all corners. 1.8.0 is available for download.

[hamish]

Replying to kyle:

Have you tried version 1.8?

not yet

According to your link, you are just using gdal_translate to create tif files.

well really I'm trying to do much much more than that, conversion to geotiff using gdal command line tools is just used as a most basic example. Really we're linking to libgdal directly.

RFC 33: GTiff - Correcting PixelIsPoint Interpretation (Adopted) fixed how GTiff handles the origin registration.

Interesting, it will be worth a try rebuilding GRASS against the new GDAL version and seeing if r.in.gdal shows the same problem. I suspect it may because GRIB driver -> GRASS doesn't pass through GeoTiff, but maybe my earlier test did have a gdal_translate step thrown in? So long ago now, I can't remember. will need to retest.

also because the data source for my test data is not by some authoritative source such as the WMO, and just some random website, there's a good chance that the source was using GDAL too, any maybe they had a geotiff pass-through step at some point? no way of knowing. but if someone knows a link to an authoritative data source for known-good GRIB files it would be useful.

thanks for the tip, Hamish

[hijmans]

Replying to hamish:

Replying to kyle:

Have you tried version 1.8? According to your link, you are just using gdal_translate to create tif files.

well really I'm trying to do much much more than that, conversion to geotiff using gdal command line tools is just used as a most basic example. Really we're linking to libgdal directly.

Same here. I am getting a wrong georeference (half a pixel shift) from gdalinfo. This is not related to geoTiff files.

also because the data source for my test data is not by some authoritative source such as the WMO, and just some random website, there's a good chance that the source was using GDAL too, any maybe they had a geotiff pass-through step at some point? no way of knowing. but if someone knows a link to an authoritative data source for known-good GRIB files it would be useful.

Same here, but I opened the the file I had with degrib ( www.nws.noaa.gov/mdl/degrib/ ) and exported it to a shapefile and to a csv file, and the results were good, not shifted. So if degrib (from NOAA) gets it right, it seems clear that gdal gets it wrong.

[Kyle Shannon]

Okay, thanks. I misunderstood the issue. I will take a quick look at this.

[hamish]

Hi,

for the record I get the same bounds from gdal 1.8.0's gdalinfo as I do with r.in.gdal in GRASS 6.5svn using libgdal 1.6.3.

(testfile: NARR grib file in LCC from NOAA, downloaded from: ​http://nomads.ncdc.noaa.gov/cgi-bin/ncdc-ui/define-collection.pl?model_sys=narr&model_name=narr-a&grid_name=221 )

no land mask on that so it's hard to check for alignment.

---

more sample data here, with a land mask applied: ​http://gribs.ocens.net/NewZealand.wind.grb.bz2

others here:

​http://gribs.ocens.net/

no idea about the quality of the ocens.net source's data, process, or land masking.

but ZyGrib (zygrib.org) shows the land and the grids lined up seemingly correctly for this ocens.net data, while gdal 1.6.3 (debian/stable) seemingly misaligns it, as in the earlier attached screenshots in this ticket.

next stop in the testing: ​http://www.nws.noaa.gov/mdl/degrib/

Hamish

[ejones]

I've been examining this code carefully, and I believe hamish is correct -- this is a bug in the GDAL Grib driver.

From my understanding of Grib format, it's clear to me that it references cell centers. Consider a GFS file (for example, the sample "gfstmp.grib2" that I attached on ticket #2550). Here's an excerpt of what I get when querying the grid definition using degrib:

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)

This says the file extends from 90N to 90S. There are two possible interpretations of that:

(1) The first row might extend from 90N to 89N, and the last row from 89S to 90S. This would also imply a row from 1N to the equator, and one from the equator down to 1S, but no row centered on the equator. Thus, 180 total rows.

(2) The first row is "centered" at 90N. With one row centered on each latitude above and below the equator, and also one row centered on the equator itself, this yields 181 total rows.

Since the header also explicitly says there are 181 rows, then the latter interpretation must be correct.

However, GDAL apparently follows the convention that the coordinates reference the cell corner. For example, the documentation for GDALDataset::GetGeoTransform states the following: "The upper left corner of the upper left pixel is at position (padfTransform[0],padfTransform[3])". So after reading a Grib file, the Grib driver needs to translate the corner by a half pixel; but, as has been pointed out here, it does not do so. So if I query gfstmp.grib2 using gdalinfo, I get the following:

Origin = (0.000000000000000,90.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Corner Coordinates:
Upper Left  (   0.0000000,  90.0000000) (  0d 0' 0.01"E, 90d 0' 0.00"N)
Lower Left  (       0.000,     -91.000) (  0d 0' 0.01"E, 91d 0' 0.00"S)
Upper Right (     360.000,      90.000) (360d 0' 0.00"E, 90d 0' 0.00"N)
Lower Right (     360.000,     -91.000) (360d 0' 0.00"E, 91d 0' 0.00"S)
Center      ( 180.0000000,  -0.5000000) (180d 0' 0.00"E,  0d30' 0.00"S)

The values shown for "origin" and "pixel size" are what it puts in the corresponding geo transform fields. As you can see, gdalinfo then assumes that these are cell-corner references, and going by that assumption it calculates the bottom of the grid to be 91S, which is not correct.

It looks to me like the patch submitted here is the correct way to fix this. My only concern is that it could potentially break software that has already been designed to work around this. For example, when writing our own display, I initially interpreted GetGeoTransform as referencing cell centers, since that's what it appeared to be doing (and because I was only working with Grib data). Once this is fixed, we can easily modify our display routines accordingly; but I do wonder if this will affect some users who are not aware of the issue.

[aboudreault]

I think that this issue had to be fixed in gdal anyway, even if some users were using it with a workaround. We should add a note in the next release document about that change. What do you think EvenR?

Thanks hamish for the patch and others for their useful comments! Fixed and committed in r23301.

[Even Rouault]

Yes, this could be mentionned in the "Backward compatibility issues" section of the NEWS

[warmerdam]

Is this addressed in the NEWS file? I notice the changes are in trunk so I've updated the milestone to 1.9.0.

[Even Rouault]

This was mentionned in the NEWS, but not in "Backward compatibility issues" section. Done in r23622