GRIB Driver breaks when built on GDAL 3.4

[runette]

WHen MDAL is built on GDAL 3.4 - the GRIB tests fail with the following message

11: C:\Users\runes\Documents\GitHub\MDAL\tests\test_gdal_grib.cpp(214): error: Expected equality of these values: 11: -0.818756103515625 11: value 11: Which is: -0.535552978515625 This is consistent across all platforms that I have tried.

I believe that this might be related to https://github.com/OSGeo/gdal/issues/4524 since - if you run gdalinfo on the test GRIB file you get different results for different versions of GDAL - note the origin and corner coordinates. I do not know enough about the GRIB driver to know what effect this change would have.

IN 3.3.3 you get Warning: Inside GRIB2Inventory, Message # 3 ERROR: Couldn't find 'GRIB' or 'TDLP' Driver: GRIB/GRIdded Binary (.grb, .grb2) Files: wind_only_u_component.grib Size is 480, 241 Coordinate System is: GEOGCRS["Coordinate System imported from GRIB file", DATUM["unnamed", ELLIPSOID["Sphere",6367470,0, LENGTHUNIT["metre",1, ID["EPSG",9001]]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433, ID["EPSG",9122]]], CS[ellipsoidal,2], AXIS["latitude",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433, ID["EPSG",9122]]], AXIS["longitude",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433, ID["EPSG",9122]]]] Data axis to CRS axis mapping: 2,1 Origin = (-0.375000000000000,90.375000000000000) Pixel Size = (0.750000000000000,-0.750000000000000) Corner Coordinates: Upper Left ( -0.3750000, 90.3750000) ( 0d22'30.00"W, 90d22'30.00"N) Lower Left ( -0.3750000, -90.3750000) ( 0d22'30.00"W, 90d22'30.00"S) Upper Right ( 359.625, 90.375) (359d37'30.00"E, 90d22'30.00"N) Lower Right ( 359.625, -90.375) (359d37'30.00"E, 90d22'30.00"S) Center ( 179.6250000, 0.0000000) (179d37'30.00"E, 0d 0' 0.01"N) Band 1 Block=480x1 Type=Float64, ColorInterp=Undefined Description = 0[-] SFC (Ground or water surface) Metadata: GRIB_COMMENT=10 metre u wind component [m/s] GRIB_ELEMENT=10U GRIB_FORECAST_SECONDS=43200 sec GRIB_REF_TIME= 1538352000 sec UTC GRIB_SHORT_NAME=0-SFC GRIB_UNIT=[m/s] GRIB_VALID_TIME= 1538395200 sec UTC Band 2 Block=480x1 Type=Float64, ColorInterp=Undefined Description = 0[-] SFC (Ground or water surface) Metadata: GRIB_COMMENT=10 metre u wind component [m/s] GRIB_ELEMENT=10U GRIB_FORECAST_SECONDS=43200 sec GRIB_REF_TIME= 1538395200 sec UTC GRIB_SHORT_NAME=0-SFC GRIB_UNIT=[m/s] GRIB_VALID_TIME= 1538438400 sec UTC In 3.4 you get Warning: Inside GRIB2Inventory, Message # 3 ERROR: Couldn't find 'GRIB' or 'TDLP' Driver: GRIB/GRIdded Binary (.grb, .grb2) Files: wind_only_u_component.grib Size is 480, 241 Coordinate System is: GEOGCRS["Coordinate System imported from GRIB file", DATUM["unnamed", ELLIPSOID["Sphere",6367470,0, LENGTHUNIT["metre",1, ID["EPSG",9001]]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433, ID["EPSG",9122]]], CS[ellipsoidal,2], AXIS["latitude",north, ORDER[1], ANGLEUNIT["degree",0.0174532925199433, ID["EPSG",9122]]], AXIS["longitude",east, ORDER[2], ANGLEUNIT["degree",0.0174532925199433, ID["EPSG",9122]]]] Data axis to CRS axis mapping: 2,1 Origin = (-180.375000000000000,90.375000000000000) Pixel Size = (0.750000000000000,-0.750000000000000) Corner Coordinates: Upper Left (-180.3750000, 90.3750000) (180d22'30.00"W, 90d22'30.00"N) Lower Left (-180.3750000, -90.3750000) (180d22'30.00"W, 90d22'30.00"S) Upper Right ( 179.6250000, 90.3750000) (179d37'30.00"E, 90d22'30.00"N) Lower Right ( 179.6250000, -90.3750000) (179d37'30.00"E, 90d22'30.00"S) Center ( -0.3750000, 0.0000000) ( 0d22'30.00"W, 0d 0' 0.01"N) Band 1 Block=480x1 Type=Float64, ColorInterp=Undefined Description = 0[-] SFC (Ground or water surface) Metadata: GRIB_COMMENT=10 metre u wind component [m/s] GRIB_ELEMENT=10U GRIB_FORECAST_SECONDS=43200 GRIB_REF_TIME=1538352000 GRIB_SHORT_NAME=0-SFC GRIB_UNIT=[m/s] GRIB_VALID_TIME=1538395200 Band 2 Block=480x1 Type=Float64, ColorInterp=Undefined Description = 0[-] SFC (Ground or water surface) Metadata: GRIB_COMMENT=10 metre u wind component [m/s] GRIB_ELEMENT=10U GRIB_FORECAST_SECONDS=43200 GRIB_REF_TIME=1538395200 GRIB_SHORT_NAME=0-SFC GRIB_UNIT=[m/s] GRIB_VALID_TIME=1538438400

[vcloarec]

@runette , good catch, big chance this could be this change on GDAL. iIt seems GDAL does no handle the grid on the same way geographic coordinates with GRIB, and indexing in MDAL is not the same now. I think the issue in the test comes from an offset in the indexes. I don't have the last version of GDAL install right now, can you try to change the line 213 of the related test file with: double value = getValue( ds, 1880 );

[runette]

Nope - 1880 did not work. it got a different value (obviously) but not the value expected.

[PeterPetrik]

could be caused by https://github.com/OSGeo/gdal/issues/4524

[vcloarec]

1880 did not work

oops, it could be 1840

[runette]

1840 worked

On Fri, 26 Nov 2021 at 12:19, Vincent Cloarec @.***> wrote:

1880 did not work

oops, it could be 1840

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/lutraconsulting/MDAL/issues/391#issuecomment-979938287, or unsubscribe https://github.com/notifications/unsubscribe-auth/AARC2MZCJ2GHLN3NNJGV3XDUN53NDANCNFSM5IZPBZUA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

[vcloarec]

1840 worked

great! So we got it, this is GDAL change. Could be fixed with a version macro:

#if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3,4,0)
double value = getValue( ds, 1840 );
#else
double value = getValue( ds, 1600 );
#endif

[PeterPetrik]

this would fix the test, but is the result still valid? (so only change is renumbering of the vertices/faces, but results are the same??)

[runette]

I managed to get v0.9.0 building on GDAL 3.4 using the GDAL config switch before running the tests

GRIB_ADJUST_LONGITUDE_RANGE=NO

as shown here https://github.com/OSGeo/gdal/pull/4526/files/f1ccd2b72e57d957bcc39ab71373559b2b4c9a97#diff-36ccca1f7d38c764146f7d197fdc78aad21e04193b6e2256daf0d3b916c376fa.

I propose :

1 - that we accept the GDAL change as the correct interpretation. However - we note in the GRIB driver page that this change has happened and reference the config switch needed to return to the old behaviour

2 - that we leave v0.9.0 (and 0.9.1 I guess) as they are - the release notes of both should be changed to note that the tests will fail but that is something only the repo owners can change.

3 - that future versions should have the GRIB tests changed as per Vince's change so that the tests pass successfully - not totally sure how to do that so that the tests work on both GDAL v3.3 and GDAL v3.4 ... ?

[vcloarec]

I propose 3 after checking if it is vizualisation OK under QGIS with GDAL 3.4?

[vcloarec]

I hadd checked the visualizationdifference in QGIS between GDAL version <3.4 and GDAL 3.4 when opening a GRIB file covering the whole globe. There is a little difference coming from a different approach for treating with geographical coordinate over 180.

Before 3.4, GDAL does not shift 0/360° (leading to the fix in 3.4) to -180/180 and MDAL does it. After 3.4, GDAL shifts the dataset but not exactly the same way than MDAL.

The difference is tiny but exists. For example, the latitude extent of a mesh covering whole world (pixel width 0.75°):

• MDAL with GDAl <3.4 : -179.25 to 180.0
• MDAL with GDAL 3.4 : -180.0 to 179.25

The absolute positions of value are the same, there is only a swap between last/first column

So here the problem:

• If we don't care of the inconsistency of MDAL results before and after GDAL 3.4 (only the extent change), then there is no problem with new GDAL 3.4 , and we can also adapt MDAL to have the same shift as GDAL if needed.

• If we care, well, we can disable the shift for GDAL,returning to the old behavior, only MDAL makes the shift with its way, so we keep consistency for MDAL.

I have a preference for the first one, but not a strong opinion.

[runette]

I have no real opinion either way.

I am thinking that keeping consistency with GDAL and QGIS is easier to explain

[PeterPetrik]

+1 to consistency

[vcloarec]

+1 to consistency

consistency MDAL/GDAL or consistency MDAL old/new ?