[gdal-dev] GRIB2 file extent not matchi GRIB2 metadata (grib_dump)

[Even Rouault]

Hi Louis-Philippe,

yes this is a matter of conventions. GRIB indeed uses a center-of-pixel 
registration whereas GDAL uses a top-left corner of pixel ones. So GDAL 
adds a half-pixel shift to expose GRIB georeferencing using its 
convention, which is a feature. What is annoying is that it leads to 
latitudes > 90 for such datasets, but there isn't much that can be done 
in the driver itself.

Except doing post-processing to subset the file to crop the top and 
bottom line (gdal_translate -srcwin 0 1 2400 1199), which will result in 
longitudes in [-90 + 0.15, 90 - 0.15].   That said that file is a bit 
inconsistent. It is quite expected that it uses an odd number of rows 
(1201) to encode a clean (90 - -90) / (1201 - 1) = 0.15 vertical 
resolution. But for the horizontal dimension, it uses an even number of 
columns (2400), which leads to a compute resolution of (360 - 0) / (2400 
- 1) = 0.15006252605252188 resolution.

Or if you're OK to deal with a slight alteration in the georeferencing 
(maximum relative error should be 0.5 / 1200 = 0.04 % at the equator), 
gdal_translate -a_ullr 0 90 360 -90

Another approach would be to crop using non-integer coordinates (to 
remove the top and bottom half-rows) + apply a non-nearest resampling like

gdal_translate -srcwin 0 0.5 2400 1200 -r cubic in.grib2 out.tif

You should get an extent of latitudes in [-90,90].

Even

Le 28/09/2023 à 21:31, Rousseau Lambert, Louis-Philippe (ECCC) via 
gdal-dev a écrit :
>
> Hi,
>
> I am facing an issue with GRIB2 data, and the extent reported by a 
> gdal_info command. Here is a sample data to test: 
> https://drive.google.com/file/d/1URvQs2qHXgRkq3YfC7ZR8VGYK3SKXB59/view?usp=drive_link
>
> Here is what I am testing with:
>
>   * Ubuntu 20.04.6 LTS
>   * GDAL 3.5.2, released 2022/09/02
>   * PROJ Rel. 8.2.0, November 1st, 2021
>
> So, using grib_dump, I see the longitude of the file should be from 0 
> to 360 degress:
>
> /grib_dump out.grib2 | grep shapeOfTheEarth -A 12/
>
> /  shapeOfTheEarth = 6;/
>
> /  Ni = 2400;/
>
> /  Nj = 1201;/
>
> /  iScansNegatively = 0;/
>
> /  jScansPositively = 1;/
>
> /  jPointsAreConsecutive = 0;/
>
> /  alternativeRowScanning = 0;/
>
> /  latitudeOfFirstGridPointInDegrees = -90;/
>
> /*longitudeOfFirstGridPointInDegrees = 0;*/
>
> /  latitudeOfLastGridPointInDegrees = 90;/
>
> /*longitudeOfLastGridPointInDegrees = 360;*/
>
> /  iDirectionIncrementInDegrees = 0.15;/
>
> /  jDirectionIncrementInDegrees = 0.15;/
>
> But using gdal_info I see the extent as:
>
> /gdal_info -proj4 out.grib2 | grep PROJ -A 9 /
>
> /PROJ.4 string is:/
>
> /'+proj=longlat +R=6371229 +no_defs'/
>
> /Origin = (-0.075031263026261,90.075000000000003)/
>
> /Pixel Size = (0.150062526052522,-0.150000000000000)/
>
> /Corner Coordinates:/
>
> /Upper Left  (  -0.0750313,  90.0750000) (  0d 4'30.11"W, 90d 4'30.00"N)/
>
> /Lower Left  (  -0.0750313, -90.0750000) (  0d 4'30.11"W, 90d 4'30.00"S)/
>
> /Upper Right (     360.075,      90.075) (360d 4'30.11"E, 90d 4'30.00"N)/
>
> /Lower Right (     360.075,     -90.075) (360d 4'30.11"E, 90d 4'30.00"S)/
>
> /Center      ( 180.0000000,   0.0000000) (180d 0' 0.00"E,  0d 0' 0.01"N)/
>
> I would expect the gdal_info extent to be also from 0 to 360 degrees 
> (2400 pixels * 0.15 = 360). From reading various source, I see the 
> issue is potentially that the GRIB2 files define the extent using the 
> center of the pixel while gdal assumes it is the top left corner.
>
> My question: Is there a way to change this behavior in gdal, or this 
> expected, or a bug ?
>
> Thanks
>
> LP
>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev

-- 
http://www.spatialys.com
My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20230928/577eb22a/attachment-0001.htm>