WCS GetCoverage - Problem with GRIDORIGIN and GRIDOFFSETS with latitude/longitude CRS (EPSG:4326 and EPSG:4617)

[jsmoreau]

What is the bug or the crash?

We host a WCS service available here: https://datacube.services.geo.ca/wrapper/ogc/elevation-hrdem-mosaic?version=1.1.1

On QGIS 3.22, we are able to load the layers as expected.

Using recent versions of QGIS (3.32). It does not work for CRS 4617 and 4326. I get an "repair source" icon. If I try to repair, QGIS starts to make request using 3979 projection and then reproject back the result in 4326.

When comparing the requests sent to the WCS server between 3.22 and 3.32 version. I see the following differences:

Request for 3.22 https://datacube.services.geo.ca/wrapper/ogc/elevation-hrdem-mosaic?SERVICE=WCS&VERSION=1.1.1&REQUEST=GetCoverage&FORMAT=image/geotiff&IDENTIFIER=dtm&BOUNDINGBOX=33.74051737904942172,-157.65736828933884794,76.42004209835096162,-26.7718522720007357,urn:ogc:def:crs:EPSG::4326&GRIDBASECRS=urn:ogc:def:crs:EPSG::4326&GRIDCS=urn:ogc:def:cs:OGC:0.0:Grid2dSquareCS&GRIDTYPE=urn:ogc:def:method:WCS:1.1:2dSimpleGrid&GRIDORIGIN=76.42004209835096162,-157.65736828933884794&GRIDOFFSETS=-0.14975271831333875,0.14975459498551269

Here, the GRIDORIGIN and GRIDOFFSETS are properly calculated. If I do, ymax - ymin / y gridoffset value, I get 285 pixels in Y The same calculation for x gives 874 pixels. The GRIDORIGIN seems to be the upper left corner which I guess is fine.

Now for 3.32

The following request is sent to the server: https://datacube.services.geo.ca/wrapper/ogc/elevation-hrdem-mosaic?SERVICE=WCS&VERSION=1.1.1&REQUEST=GetCoverage&FORMAT=image/geotiff&IDENTIFIER=dtm&BOUNDINGBOX=35.80710489177349842,-151.18048205621542479,74.35345458562689203,-33.24873850512415885,urn:ogc:def:crs:EPSG::4326&GRIDBASECRS=urn:ogc:def:crs:EPSG::4326&GRIDCS=urn:ogc:def:cs:OGC:0.0:Grid2dSquareCS&GRIDTYPE=urn:ogc:def:method:WCS:1.1:2dSimpleGrid&GRIDORIGIN=-33.24873850512415885,35.80710489177349842&GRIDOFFSETS=-4.28292774376148877,13.10352706123236111

The GRIDORIGIN coordinate does not fit with the BOUNDINGBOX coordinates. Based on this bbox, it should be 74.35345458562689203,-151.18048205621542479. The GRIDOFFSETS is also wrong. Based on the working request of 3.22, I managed to manually change the 3.32 request to make it work. This modified request now works: https://datacube.services.geo.ca/wrapper/ogc/elevation-hrdem-mosaic?SERVICE=WCS&VERSION=1.1.1&REQUEST=GetCoverage&FORMAT=image/geotiff&IDENTIFIER=dtm&BOUNDINGBOX=35.80710489177349842,-151.18048205621542479,74.35345458562689203,-33.24873850512415885,urn:ogc:def:crs:EPSG::4326&GRIDBASECRS=urn:ogc:def:crs:EPSG::4326&GRIDCS=urn:ogc:def:cs:OGC:0.0:Grid2dSquareCS&GRIDTYPE=urn:ogc:def:method:WCS:1.1:2dSimpleGrid&GRIDORIGIN=74.35345458562689203,-151.18048205621542479&GRIDOFFSETS=-0.13525034980299,0.1349333450241318

Steps to reproduce the issue

1. Set the map canvas CRS to EPSG4326
2. Load the following WCS https://datacube.services.geo.ca/wrapper/ogc/elevation-hrdem-mosaic?version=1.1.1
3. Select the terrain layer and load in canvas

Versions

QGIS version 3.32.0-Lima QGIS code revision 311a8cb8a6 Qt version 5.15.3 Python version 3.9.5 GDAL/OGR version 3.7.0 PROJ version 9.2.1 EPSG Registry database version v10.088 (2023-05-13) GEOS version 3.11.2-CAPI-1.17.2 SQLite version 3.41.1 PDAL version 2.5.3 PostgreSQL client version 15.2 SpatiaLite version 5.0.1 QWT version 6.1.6 QScintilla2 version 2.13.1 OS version Windows 10 Version 2009

Active Python plugins LAStools 1.4 planet_explorer 2.3.1 pluginbuilder3 3.2.1 qgis_stac 1.1.1 QuickWKT 3.1 rasterdataplotting 1.6.3 systeme-prod-plugin-qgis 1.0 zoomtopaste 3.0.3 db_manager 0.1.20 MetaSearch 0.3.6 processing 2.12.99

Supported QGIS version

• [X] I'm running a supported QGIS version according to the roadmap.

New profile

• [x] I tried with a new QGIS profile

Additional context

No response