diff --git a/autotest/gcore/tiff_srs.py b/autotest/gcore/tiff_srs.py index b5b08bc54c76..b5fe8d96d2af 100755 --- a/autotest/gcore/tiff_srs.py +++ b/autotest/gcore/tiff_srs.py @@ -1289,3 +1289,59 @@ def test_tiff_srs_read_compound_with_VerticalCitationGeoKey_only(): wkt == """COMPD_CS["TestMS",GEOGCS["NAD83(2011)",DATUM["NAD83_National_Spatial_Reference_System_2011",SPHEROID["GRS 1980",6378137,298.257222101004,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","1116"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","6318"]],VERT_CS["NAVD88 height",VERT_DATUM["North American Vertical Datum 1988",2005,AUTHORITY["EPSG","5103"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Up",UP]]]""" ) + + +@pytest.mark.parametrize( + "code", + [ + 7415, # Amersfoort / RD New + NAP height + 9707, # WGS 84 + EGM96 height + ], +) +@pytest.mark.require_proj( + 7, 2 +) # not necessarily the minimum version, but 9707 doesn't exist in PROJ 6.x +def test_tiff_srs_read_compound_with_EPSG_code(code): + + """Test bugfix for https://github.com/OSGeo/gdal/issues/7982""" + + filename = "/vsimem/test_tiff_srs_read_compound_with_EPSG_code.tif" + srs = osr.SpatialReference() + srs.ImportFromEPSG(code) + srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) + ds = gdal.GetDriverByName("GTiff").Create(filename, 1, 1) + ds.SetSpatialRef(srs) + ds = None + ds = gdal.Open(filename) + gdal.ErrorReset() + got_srs = ds.GetSpatialRef() + assert gdal.GetLastErrorMsg() == "", srs.ExportToWkt(["FORMAT=WKT2_2019"]) + assert got_srs.GetAuthorityCode(None) == str(code) + assert got_srs.IsSame(srs) + ds = None + gdal.Unlink(filename) + + +def test_tiff_srs_read_compound_without_EPSG_code(): + + """Test case where identification of code for CompoundCRS (added for + bugfix of https://github.com/OSGeo/gdal/issues/7982) doesn't trigger""" + + filename = "/vsimem/test_tiff_srs_read_compound_without_EPSG_code.tif" + srs = osr.SpatialReference() + # WGS 84 + NAP height, unlikely to have a EPSG code ever + srs.SetFromUserInput("EPSG:4326+5709") + srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) + ds = gdal.GetDriverByName("GTiff").Create(filename, 1, 1) + ds.SetSpatialRef(srs) + ds = None + ds = gdal.Open(filename) + gdal.ErrorReset() + got_srs = ds.GetSpatialRef() + assert gdal.GetLastErrorMsg() == "", srs.ExportToWkt(["FORMAT=WKT2_2019"]) + assert got_srs.GetAuthorityCode(None) is None + assert got_srs.GetAuthorityCode("GEOGCS") == "4326" + assert got_srs.GetAuthorityCode("VERT_CS") == "5709" + assert got_srs.IsSame(srs) + ds = None + gdal.Unlink(filename) diff --git a/frmts/gtiff/gt_wkt_srs.cpp b/frmts/gtiff/gt_wkt_srs.cpp index 488874b1c820..7895aaf6d2f0 100644 --- a/frmts/gtiff/gt_wkt_srs.cpp +++ b/frmts/gtiff/gt_wkt_srs.cpp @@ -1689,6 +1689,8 @@ OGRSpatialReferenceH GTIFGetOGISDefnAsOSR(GTIF *hGTIF, GTIFDefn *psDefn) if (bCanBuildCompoundCRS) { + const bool bHorizontalHasCode = + oSRS.GetAuthorityCode(nullptr) != nullptr; const char *pszHorizontalName = oSRS.GetName(); const std::string osHorizontalName( pszHorizontalName ? pszHorizontalName : "unnamed"); @@ -1735,6 +1737,19 @@ OGRSpatialReferenceH GTIFGetOGISDefnAsOSR(GTIF *hGTIF, GTIFDefn *psDefn) { oSRS.GetRoot()->AddChild(oVertSRS.GetRoot()->Clone()); bNeedManualVertCS = false; + + // GeoTIFF doesn't store EPSG code of CompoundCRS, so + // if we have an EPSG code for both the horizontal and vertical + // parts, check if there's a known CompoundCRS associating + // both + if (bHorizontalHasCode && verticalCSType != KvUserDefined && + verticalCSType > 0) + { + const auto *poSRSMatch = oSRS.FindBestMatch(100); + if (poSRSMatch) + oSRS = *poSRSMatch; + delete poSRSMatch; + } } } }