GeoHEIF: Fixed respecting CRS axis order (e.g., for EPSG:4326)

[jerstlouis]

What does this PR do?

Fixes an axis-ordering bug in GeoHEIF driver, in particular for EPSG:4326 CRS.

What are related issues/pull requests?

https://gitlab.ogc.org/ogc/T20-GIMI/-/issues/59

Tasklist

• ADD YOUR TASKS HERE
• Make sure code is correctly formatted (cf pre-commit configuration)
• Add test case(s)
• [N/A] Add documentation
• [N/A] Updated Python API documentation (swig/include/python/docs/)
• Review
• Adjust for comments
• All CI builds and checks have passed

Environment

N/A

[rouault]

I don't understand that comment. "const" means that it does not change the observable state of the object (that is the GDALHEIFDataset instance). What observable side effects does GetSpatialRef() cause ?

[jerstlouis]

I'm not sure how the "observable state of the object" is defined in C++ and whether this "set things up the first time" is an OK practice (assuming it's done in a thread-safe manner).

The main issue is that this function might not be thread-safe, since if it's called for the first time by two different threads, they both might be trying to modify the geoHEIF member of the instance at the same time.

But as it turns out and was the source of the problem I chased for the last few hours, there was a very observable side-effect causing problems. Calling GDALGetGeoTransform() before and after calling GetSpatialRef() yielded a completely different result since geoHEIF.m_oSRS was null before and non-null after.

[rouault]

The main issue is that this function might not be thread-safe,

that's by design. No GDAL objects are thread-safe unless otherwise stated: cf https://gdal.org/en/stable/user/multithreading.html#gdal-api-re-entrant-but-generally-not-thread-safe

Calling GDALGetGeoTransform() before and after calling GetSpatialRef() yielded a completely different result since geoHEIF.m_oSRS was null before and non-null after.

That's a different thing, an implementation bug. I believe the comment in GetGeoTransform() is sufficient, and the one in GetSpatialRef() can be removed

Actually looking further at that code, I believe that GDALHEIFDataset::GetGeoTransform() should not call GDALHEIFDataset::GetSpatialRef() itself. But rather GeoHEIF::GetGeoTransform() should make sure that its m_oSRS member is properly initialized if it needs it to access it, since that's the one that needs the SRS. Perhaps the right thing to do is to set m_oSRS as mutable in the GeoHEIF class itself, and make GeoHEIF::GetGeoTransform() call GeoHEIF::GetSpatialRef()

[rouault]

Perhaps the right thing to do is to set m_oSRS as mutable in the GeoHEIF class itself, and make GeoHEIF::GetGeoTransform() call GeoHEIF::GetSpatialRef()

yes we definitely want to do that since the GeoHEIF class is also used by the AVIF driver, and thus we don't want to do the same hack in both drivers.

[jerstlouis]

Wouldn't it be simpler / preferable to automatically set up the SRS when initializing the dataset, before either of these methods can get called?

No GDAL objects are thread-safe unless otherwise stated

Even after reading this, I would have been very surprised that a method declared const behaves this way.

[rouault]

Wouldn't it be simpler / preferable to automatically set up the SRS when initializing the dataset, before either of these methods can get called?

Actually I do see my suggestion of calling GeoHEIF::GetSpatialRef() from GeoHEIF::GetGeoTransform() woudn't work since the geoHEIF.extractSRS() must be called with the relevant metadata first. So, I'm backing away my latest suggestion. I believe you're current fix is OK. (and the AVIF driver doesn't need a fix since it has a method processProperties() that feeds all needed metadata to the GeoHEIF instance)

[jerstlouis]

I've removed my // REVIEW: comment, and I'll let you or @bradh figure out the best way to address the GetSpatialRef() behavior.

[rouault]

This would deserve a new test case in autotest/gdrivers/heif.py reading a file with EPSG:4326 cheking that:
assert srs.GetDataAxisToSRSAxisMapping() == [2, 1]
and that the ds.GetGeoTransform() is such that its coefficient at index 0 is a longitude
to use the "GIS friendly order" used by all other GDAL raster drivers

[rouault]

CC @bradh

[rouault]

[coveralls]

coverage: 74.178% (+0.006%) from 74.172%
when pulling 87f0de1 on jerstlouis:master
into 362584f on OSGeo:master.

[bradh]

I think it would be cleaner to just do the HEIF properties the same way I did the AVIF properties.

Then we can look at the axis order problem separately.

[jerstlouis]

@rouault , @desruisseaux mentioned:

The axes[0] != 1 check should be abs(axes[0]) != 1, because GDAL use negative numbers as a convention meaning that axis direction should be reversed (e.g. south -> north).

is this the case that GetDataAxisToSRSAxisMapping() could return negative values here?

Thanks!

[rouault]

is this the case that GetDataAxisToSRSAxisMapping() could return negative values here?

in theory yes, but I'm not smart enough to have figured out what the OAMS_TRADITIONAL_GIS_ORDER should be for the weird southing,westing or westing,southing CRS in the context of raster data should be (at least in the GeoTIFF context, people use the CRS variants that use northing-easting or easting-northing axis order to avoid those headaches. southing,westing or westing,southing are left only for purely geodesy non-mapping use cases), so in practice you won't get negative values, unless you manually set them with SetDataAxisToSRSAxisMapping(). I would probably not worry too much about that...

[jerstlouis]

I've converted this to a non-draft PR since @desruisseaux 's recent GDAL implementation guidance comments / pseudo-code only relate to the rare cases of southing / westing for pure geodesy, and multi-dimensional data which is not yet supported by the GeoHEIF driver anyways.

If you don't mind I will let any of you @bradh or @rouault to add the test case -- please feel free to add yourselves as co-author of this commit. Also fine to fix reading the properties as @bradh suggested first, and then we wouldn't need that hackish call to GetSpatialRef().
Also feel free to derive test data from our GNOSIS Map Server that now generates GeoHEIF in multiple CRSs e.g.:

https://maps.gnosis.earth/ogcapi/collections/blueMarble/map?crs=[OGC:CRS84]&f=heif

https://maps.gnosis.earth/ogcapi/collections/blueMarble/map?crs=[EPSG:4326]&f=heif

https://maps.gnosis.earth/ogcapi/collections/blueMarble/map?crs=[EPSG:3395]&f=heif

https://maps.gnosis.earth/ogcapi/collections/blueMarble/map?crs=[EPSG:3857]&f=heif

Currently, this is only for /map but we should add support for /coverage and /tiles this coming week, as well as the band information and OGC branding that's still missing.

I've tested all these outputs with gdalinfo and QGIS and everything seems fine.

Thanks all!