Skip to content

Commit

Permalink
GeoTIFF SRS reader: try to retrieve the EPSG code of a CompoundCRS if…
Browse files Browse the repository at this point in the history
… the one of its horizontal and vertical part is known (fixes #7982)
  • Loading branch information
rouault committed Jun 20, 2023
1 parent 483fc3b commit bc36ae3
Show file tree
Hide file tree
Showing 2 changed files with 71 additions and 0 deletions.
56 changes: 56 additions & 0 deletions autotest/gcore/tiff_srs.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
15 changes: 15 additions & 0 deletions frmts/gtiff/gt_wkt_srs.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand Down Expand Up @@ -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;
}
}
}
}
Expand Down

0 comments on commit bc36ae3

Please sign in to comment.