Opened 8 years ago

Closed 5 years ago

#2637 closed defect (fixed)

grib driver: grid vs cell-center convention wrong?

Reported by: hamish Owned by: warmerdam
Priority: normal Milestone: 1.9.0
Component: GDAL_Raster Version: svn-trunk
Severity: normal Keywords: grib
Cc: Kyle Shannon, bfraser, dmorissette, mprecise, Even Rouault

Description (last modified by 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.

Attachments (3)

grib_half_cell_shift.png (112.4 KB) - added by hamish 8 years ago.
grib cells shifted 1/2 to SE
gdal_grib_node.patch (618 bytes) - added by tbouffon 8 years ago.
gdal_grib_node_june2010_trunk.patch (588 bytes) - added by hamish 7 years ago.
updated patch vs. latest svn trunk for grid registration error

Download all attachments as: .zip

Change History (24)

Changed 8 years ago by hamish

Attachment: grib_half_cell_shift.png added

grib cells shifted 1/2 to SE

comment:1 Changed 8 years ago by hamish

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

Hamish

comment:2 Changed 8 years ago by Kyle Shannon

Cc: Kyle Shannon added

comment:3 Changed 8 years ago by tbouffon

Resolution: worksforme
Status: newclosed

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 !

Changed 8 years ago by tbouffon

Attachment: gdal_grib_node.patch added

comment:4 Changed 8 years ago by tbouffon

Resolution: worksforme
Status: closedreopened

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

Changed 7 years ago by hamish

updated patch vs. latest svn trunk for grid registration error

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

comment:6 Changed 7 years ago by hamish

sample data in #3629 (a duplicate)

comment:7 in reply to:  6 Changed 6 years ago by rhijmans

Description: modified (diff)

Replying to hamish:

sample data in #3629 (a duplicate)

comment:8 Changed 6 years ago by hamish

Description: modified (diff)

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

comment:9 in reply to:  8 ; Changed 6 years ago by 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.

comment:10 in reply to:  9 ; Changed 6 years ago by 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

comment:11 in reply to:  10 Changed 6 years ago by 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.

comment:12 Changed 6 years ago by Kyle Shannon

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

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

comment:14 Changed 6 years ago by bfraser

Cc: bfraser added

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

comment:16 Changed 6 years ago by dmorissette

Cc: dmorissette added

comment:17 Changed 5 years ago by mprecise

Cc: mprecise added

comment:18 Changed 5 years ago by aboudreault

Cc: Even Rouault added

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.

comment:19 Changed 5 years ago by Even Rouault

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

comment:20 Changed 5 years ago by warmerdam

Milestone: 1.6.41.9.0

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

comment:21 Changed 5 years ago by Even Rouault

Resolution: fixed
Status: reopenedclosed

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

Note: See TracTickets for help on using tickets.