[gdal-dev] GRIB format GeoTransform error

[Graziano Giuliani]

Hi,

I am playing in this days with the GRIB IO layer, using the API for reading GRIB
into my application.

I have find a possible bug (I am referring to V 1.6.2, but is still there in
SVN) on the Metadata read from file in SetGribMetaData, at lines 604-605 of the
gribdataset.cpp source file.

           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

If I do not comment those two lines, I have GeoTransform wrong. Seems rotating
what is already rotated.... commenting the two lines I get correct number for
rMaxY. If I want to reverse again from South to North, shouldn't that equality
check be changed to check if different?

I have no time now to correctly investigate this, but please note that the
following code can be used to test.

GDALDataset *Dataset;
Dataset = (GDALDataset *) GDALOpen( filename, GA_ReadOnly );

static const char llpjwks[] = "GEOGCS[\"WGS 84\", DATUM[\"WGS_1984\",
SPHEROID[\"WGS 84\",6378137,298.257223563, AUTHORITY[\"EPSG\",\"7030\"]],
TOWGS84[0,0,0,0,0,0,0], AUTHORITY[\"EPSG\",\"6326\"]], PRIMEM[\"Greenwich\",0,
AUTHORITY[\"EPSG\",\"8901\"]], UNIT[\"degree\",0.0174532925199433,
AUTHORITY[\"EPSG\",\"9108\"]], AXIS[\"Lat\",NORTH], AXIS[\"Long\",EAST],
AUTHORITY[\"EPSG\",\"4326\"]]";

char *lpjinfo = strdup(llpjwks);
OGRSpatialReference LLSRS;
LLSRS.importFromWkt(&lpjinfo);

OGRSpatialReference INSRS;
char *ipjinfo = strdup(Dataset->GetProjectionRef());
INSRS.importFromWkt(&ipjinfo);
OGRCoordinateTransformation *LLIJ = OGRCreateCoordinateTransformation(&LLSRS,
&INSRS);
OGRCoordinateTransformation *IJLL = OGRCreateCoordinateTransformation(&INSRS,
&LLSRS);

double adfGeoTransform[6];
double fdaGeoTransform[6];
Dataset->GetGeoTransform(adfGeoTransform);
GDALInvGeoTransform(adfGeoTransform, fdaGeoTransform);

x = lon;
y = lat;
LLIJ->Transform(1, &x, &y);
i = fdaGeoTransform[0] + x*fdaGeoTransform[1] + y*fdaGeoTransform[2];
j = fdaGeoTransform[3] + x*fdaGeoTransform[4] + y*fdaGeoTransform[5];

On a GFS regular lat/lon grid 0.5 as from NOMADS
(http://nomads.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/gfs.2009111300/gfs.t00z.pgrb2bf00
for example), I have adfGeoTransform[3] (rMaxY) at 270. Should be 90, and
commenting out the two lines in gribdataset.cpp I have the correct answer (and
everything works as expected, i.e. for example for latitude 45, longitude 45 I
have I = 90, J = 90.

Thanks.

Graziano.