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@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
--
http://www.spatialys.com
My software is free, but my time generally not.
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
