Revise how SRS methods deal with TOWGS84

When importing from EPSG, GDAL 3.0.0 to 3.0.2 followed the scheme used in GDAL 1.x and GDAL 2.x, that is they tried to attach a TOWGS84 transformation. While the intent was to have some sort of backward compatibility, this may be a pain for the future. So do the following changes: - importFromEPSG(): no longer attach a TOWGS84 transformation, unless the user set the OSR_ADD_TOWGS84_ON_IMPORT_FROM_EPSG=YES configuration option - Add a OGRSpatialReference::AddGuessedTOWGS84() to attach such a transformation, when possible (note: this will do it in a more prudent way than GDAL 1.x/2.x, that is only if a transformation is found for the whole area of use of the CRS) - exportToProj4(): if the SRS has no transformation to WGS84, attach a +towgs84 if the SRS has a EPSG code and AddGuessedToWGS84() succeeds. Can be disabled with the OSR_ADD_TOWGS84_ON_EXPORT_TO_PROJ4 = NO configuration option - exportToWkt() with WKT1 format: Add a OSR_ADD_TOWGS84_ON_EXPORT_TO_WKT1 = YES/NO configuration option, which defaults to NO. If set to YES, then a transformation to WGS84 using AddGuessedToWGS84() logic is researched to add a TOWGS84[] node.
OSGeo · Dec 17, 2019 · cc02dc4 · cc02dc4
1 parent bc83793
commit cc02dc4
Show file tree

Hide file tree

Showing 17 changed files with 296 additions and 82 deletions.
diff --git a/autotest/gcore/tiff_srs.py b/autotest/gcore/tiff_srs.py
@@ -650,6 +650,7 @@ def test_tiff_srs_towgs84_from_epsg_do_not_write_it():
     ds = gdal.GetDriverByName('GTiff').Create(filename, 1, 1)
     srs_in = osr.SpatialReference()
     srs_in.ImportFromEPSG(31468)
+    srs_in.AddGuessedTOWGS84()
     assert srs_in.HasTOWGS84()
     ds.SetSpatialRef(srs_in)
     ds = None
@@ -666,6 +667,7 @@ def test_tiff_srs_towgs84_from_epsg_force_write_it():
     ds = gdal.GetDriverByName('GTiff').Create(filename, 1, 1)
     srs_in = osr.SpatialReference()
     srs_in.ImportFromEPSG(31468)
+    srs_in.AddGuessedTOWGS84()
     assert srs_in.HasTOWGS84()
     with gdaltest.config_option('GTIFF_WRITE_TOWGS84', 'YES'):
         ds.SetSpatialRef(srs_in)

diff --git a/autotest/gcore/tiff_write.py b/autotest/gcore/tiff_write.py
@@ -2004,7 +2004,7 @@ def test_tiff_write_64():
     wkt = ds.GetProjection()
     ds = None
 
-    expected_wkt = """GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]"""
+    expected_wkt = """GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]"""
 
     assert wkt == expected_wkt, 'coordinate system does not exactly match.'
 

diff --git a/autotest/gdrivers/dted.py b/autotest/gdrivers/dted.py
@@ -148,8 +148,7 @@ def test_dted_7():
 
     assert gdal.GetLastErrorMsg() is not None, 'An expected warning was not emitted'
 
-    assert prj == 'GEOGCS["WGS 72",DATUM["World_Geodetic_System_1972",SPHEROID["WGS 72",6378135,298.26]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4322"]]', \
-        ('Projection does not match expected:\n%s' % prj)
+    assert prj.startswith('GEOGCS["WGS 72"')
 
 ###############################################################################
 # Test a file whose checksum is corrupted

diff --git a/autotest/ogr/ogr_shape.py b/autotest/ogr/ogr_shape.py
@@ -3821,7 +3821,6 @@ def test_ogr_shape_99():
         DATUM["CH1903",
             SPHEROID["Bessel 1841",6377397.155,299.1528128,
                 AUTHORITY["EPSG","7004"]],
-            TOWGS84[674.374,15.056,405.346,0,0,0,0],
             AUTHORITY["EPSG","6149"]],
         PRIMEM["Greenwich",0,
             AUTHORITY["EPSG","8901"]],

diff --git a/autotest/osr/osr_compd.py b/autotest/osr/osr_compd.py
@@ -178,7 +178,6 @@ def test_osr_compd_5():
             DATUM["North_American_Datum_1983",
                 SPHEROID["GRS 1980",6378137,298.257222101,
                     AUTHORITY["EPSG","7019"]],
-                TOWGS84[0,0,0,0,0,0,0],
                 AUTHORITY["EPSG","6269"]],
             PRIMEM["Greenwich",0,
                 AUTHORITY["EPSG","8901"]],
@@ -211,7 +210,7 @@ def test_osr_compd_5():
         print('warning they are equivalent, but not completely the same')
         print(wkt)
 
-    exp_proj4 = '+proj=utm +zone=11 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +vunits=m +no_defs'
+    exp_proj4 = '+proj=utm +zone=11 +datum=NAD83 +units=m +vunits=m +no_defs'
     proj4 = srs.ExportToProj4()
     assert proj4 == exp_proj4, ('Did not get expected proj.4 string, got:' + proj4)
 

diff --git a/autotest/osr/osr_ct_proj.py b/autotest/osr/osr_ct_proj.py
@@ -131,7 +131,17 @@
      'EPSG:4326', (49.9988573027651,9.99881145557889, 0.0), 1e-8,
      'DHDN -> WGS84 using BETA2007', None, 'GRID:BETA2007.gsb'),
 
-    ('EPSG:4314', (50, 10, 0.0), 1e-8,
+    ("""GEOGCS["DHDN",
+    DATUM["Deutsches_Hauptdreiecksnetz",
+        SPHEROID["Bessel 1841",6377397.155,299.1528128,
+            AUTHORITY["EPSG","7004"]],
+        TOWGS84[598.1,73.7,418.2,0.202,0.045,-2.455,6.7],
+        AUTHORITY["EPSG","6314"]],
+    PRIMEM["Greenwich",0,
+        AUTHORITY["EPSG","8901"]],
+    UNIT["degree",0.0174532925199433,
+        AUTHORITY["EPSG","9122"]],
+    AUTHORITY["EPSG","4314"]]""", (50, 10, 0.0), 1e-8,
      'EPSG:4326', (49.9988572643058,9.99881392529464,0), 1e-8,
      'DHDN -> WGS84 using TOWGS84 automatically set', 'OSR_CT_USE_DEFAULT_EPSG_TOWGS84=YES', None),
 

diff --git a/autotest/osr/osr_epsg.py b/autotest/osr/osr_epsg.py
@@ -47,14 +47,13 @@ def test_osr_epsg_1():
     assert srs.GetAuthorityCode(None) == '3003'
 
 ###############################################################################
-# Check that EPSG:4312 lookup has the towgs84 values set properly
-# from gcs.override.csv.
-
+# Check that EPSG:4312 w.r.t towgs84 values
 
 def test_osr_epsg_2():
 
     srs = osr.SpatialReference()
-    srs.ImportFromEPSG(4312)
+    with gdaltest.config_option('OSR_ADD_TOWGS84_ON_IMPORT_FROM_EPSG', 'YES'):
+        srs.ImportFromEPSG(4312)
 
     if float(srs.GetAttrValue('TOWGS84', 6)) != pytest.approx(2.4232, abs=0.0005):
         print(srs.ExportToPrettyWkt())
@@ -71,6 +70,7 @@ def test_osr_epsg_3():
     for epsg in [3120, 2172, 2173, 2174, 2175, 3328]:
         srs = osr.SpatialReference()
         srs.ImportFromEPSG(epsg)
+        srs.AddGuessedTOWGS84()
 
         expected_towgs84 = [33.4, -146.6, -76.3, -0.359, -0.053, 0.844, -0.84]
 

diff --git a/gdal/frmts/gtiff/geotiff.cpp b/gdal/frmts/gtiff/geotiff.cpp
@@ -12714,29 +12714,12 @@ void GTiffDataset::LookForProjection()
 
         if( GTIFGetDefn( hGTIF, psGTIFDefn ) )
         {
-            char* pszProjection = GTIFGetOGISDefn( hGTIF, psGTIFDefn );
-            if( pszProjection )
-            {
-                m_oSRS.SetFromUserInput(pszProjection);
-                double adfTOWGS84[7];
-                bool bHasTOWGS84 = m_oSRS.GetTOWGS84(&adfTOWGS84[0], 7) == OGRERR_NONE;
-                const char* pszCode = m_oSRS.GetAuthorityCode(nullptr);
-                if( pszCode )
-                {
-                    m_oSRS.importFromEPSG(atoi(pszCode));
-                    if( bHasTOWGS84 )
-                    {
-                        m_oSRS.SetTOWGS84(adfTOWGS84[0],
-                                        adfTOWGS84[1],
-                                        adfTOWGS84[2],
-                                        adfTOWGS84[3],
-                                        adfTOWGS84[4],
-                                        adfTOWGS84[5],
-                                        adfTOWGS84[6]);
-                    }
-                }
+            OGRSpatialReferenceH hSRS = GTIFGetOGISDefnAsOSR( hGTIF, psGTIFDefn );
+            if( hSRS )
+            {
+                m_oSRS = *(OGRSpatialReference::FromHandle(hSRS));
+                OSRDestroySpatialReference(hSRS);
             }
-            CPLFree(pszProjection);
 
             if( m_oSRS.IsCompound() )
             {

diff --git a/gdal/frmts/gtiff/gt_wkt_srs.cpp b/gdal/frmts/gtiff/gt_wkt_srs.cpp
@@ -305,10 +305,10 @@ int GDALGTIFKeyGetDOUBLE( GTIF *hGTIF, geokey_t key,
 }
 
 /************************************************************************/
-/*                          GTIFGetOGISDefn()                           */
+/*                      GTIFGetOGISDefnAsOSR()                          */
 /************************************************************************/
 
-char *GTIFGetOGISDefn( GTIF *hGTIF, GTIFDefn * psDefn )
+OGRSpatialReferenceH GTIFGetOGISDefnAsOSR( GTIF *hGTIF, GTIFDefn * psDefn )
 
 {
     OGRSpatialReference oSRS;
@@ -335,32 +335,27 @@ char *GTIFGetOGISDefn( GTIF *hGTIF, GTIFDefn * psDefn )
         && psDefn->Model != ModelTypeGeographic
         && psDefn->Model != ModelTypeGeocentric )
     {
-        char *pszWKT = nullptr;
         char szPeStr[2400] = { '\0' };
 
         /** check if there is a pe string citation key **/
         if( GDALGTIFKeyGetASCII( hGTIF, PCSCitationGeoKey, szPeStr,
                                  0, sizeof(szPeStr) ) &&
             strstr(szPeStr, "ESRI PE String = " ) )
         {
-            pszWKT = CPLStrdup( szPeStr + strlen("ESRI PE String = ") );
+            const char* pszWKT = szPeStr + strlen("ESRI PE String = ");
+            oSRS.importFromWkt(pszWKT);
 
             if( strstr( pszWKT,
                         "PROJCS[\"WGS_1984_Web_Mercator_Auxiliary_Sphere\"" ) )
             {
-                oSRS.importFromWkt(pszWKT);
                 oSRS.SetExtension(
                     "PROJCS", "PROJ4",
                     "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 "
                     "+x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null "
                     "+wktext  +no_defs" );  // TODO(schwehr): Why 2 spaces?
-
-                CPLFree(pszWKT);
-                pszWKT = nullptr;
-                oSRS.exportToWkt(&pszWKT);
             }
 
-            return pszWKT;
+            return OGRSpatialReference::ToHandle(oSRS.Clone());
         }
         else
         {
@@ -408,9 +403,7 @@ char *GTIFGetOGISDefn( GTIF *hGTIF, GTIFDefn * psDefn )
 
                 GTIFFreeMemory( pszUnitsName );
             }
-            oSRS.exportToWkt( &pszWKT );
-
-            return pszWKT;
+            return OGRSpatialReference::ToHandle(oSRS.Clone());
         }
     }
 
@@ -582,10 +575,7 @@ char *GTIFGetOGISDefn( GTIF *hGTIF, GTIFDefn * psDefn )
         if( CheckCitationKeyForStatePlaneUTM( hGTIF, psDefn, &oSRS,
                                               &linearUnitIsSet ) )
         {
-            oSRS.morphFromESRI();
-            char *pszWKT = nullptr;
-            if( oSRS.exportToWkt( &pszWKT ) == OGRERR_NONE )
-                return pszWKT;
+            return OGRSpatialReference::ToHandle(oSRS.Clone());
         }
 
         /* Handle ESRI PE string in citation */
@@ -733,15 +723,14 @@ char *GTIFGetOGISDefn( GTIF *hGTIF, GTIFDefn * psDefn )
                               &(psDefn->PMLongToGreenwich), 0, 1 );
     }
 
-    bool aUnitGot = false;
     if( !pszAngularUnits )
     {
 #if LIBGEOTIFF_VERSION >= 1600
         GTIFGetUOMAngleInfoEx( projContext,
 #else
         GTIFGetUOMAngleInfo(
 #endif
-            psDefn->UOMAngle, &pszAngularUnits, nullptr );
+            psDefn->UOMAngle, &pszAngularUnits, &psDefn->UOMAngleInDegrees );
         if( pszAngularUnits == nullptr )
             pszAngularUnits = CPLStrdup("unknown");
         else
@@ -753,7 +742,6 @@ char *GTIFGetOGISDefn( GTIF *hGTIF, GTIFDefn * psDefn )
         if( GDALGTIFKeyGetDOUBLE(hGTIF, GeogAngularUnitSizeGeoKey, &dfRadians,
                                  0, 1) )
         {
-            aUnitGot = true;
             psDefn->UOMAngleInDegrees = dfRadians / CPLAtof(SRS_UA_DEGREE_CONV);
         }
     }
@@ -1015,14 +1003,6 @@ char *GTIFGetOGISDefn( GTIF *hGTIF, GTIFDefn * psDefn )
         for( ; i < 10; i++ )
             adfParm[i] = 0.0;
 
-        if(!aUnitGot)
-        {
-            adfParm[0] *= psDefn->UOMAngleInDegrees;
-            adfParm[1] *= psDefn->UOMAngleInDegrees;
-            adfParm[2] *= psDefn->UOMAngleInDegrees;
-            adfParm[3] *= psDefn->UOMAngleInDegrees;
-        }
-
 /* -------------------------------------------------------------------- */
 /*      Translation the fundamental projection.                         */
 /* -------------------------------------------------------------------- */
@@ -1382,14 +1362,41 @@ char *GTIFGetOGISDefn( GTIF *hGTIF, GTIFDefn * psDefn )
             ->AddChild( new OGR_SRSNode( "UP" ) );
     }
 
-/* ==================================================================== */
-/*      Return the WKT serialization of the object.                     */
-/* ==================================================================== */
+    // Hack for tiff_read.py:test_tiff_grads so as to normalize angular
+    // parameters to grad
+    if( psDefn->UOMAngleInDegrees != 1.0 )
+    {
+        char *pszWKT = nullptr;
+        const char* const apszOptions[] = {
+            "FORMAT=WKT1", "ADD_TOWGS84_ON_EXPORT_TO_WKT1=NO", nullptr };
+        if( oSRS.exportToWkt( &pszWKT, apszOptions ) == OGRERR_NONE )
+        {
+            oSRS.importFromWkt(pszWKT);
+        }
+        CPLFree(pszWKT);
+    }
+
+    return OGRSpatialReference::ToHandle(oSRS.Clone());
+}
+
+
+/************************************************************************/
+/*                          GTIFGetOGISDefn()                           */
+/************************************************************************/
+
+char *GTIFGetOGISDefn( GTIF *hGTIF, GTIFDefn * psDefn )
+{
+    OGRSpatialReferenceH hSRS = GTIFGetOGISDefnAsOSR( hGTIF, psDefn );
 
     char *pszWKT = nullptr;
-    if( oSRS.exportToWkt( &pszWKT ) == OGRERR_NONE )
+    if( hSRS &&
+        OGRSpatialReference::FromHandle(hSRS)->exportToWkt( &pszWKT ) == OGRERR_NONE )
+    {
+        OSRDestroySpatialReference(hSRS);
         return pszWKT;
+    }
     CPLFree(pszWKT);
+    OSRDestroySpatialReference(hSRS);
 
     return nullptr;
 }
@@ -2770,11 +2777,14 @@ int GTIFSetFromOGISDefnEx( GTIF * psGTIF, const char *pszOGCWKT,
             CPLErrorHandlerPusher oErrorHandler(CPLQuietErrorHandler);
             OGRSpatialReference oRefSRS;
             double adfRefTOWGS84[7] = { 0.0 };
-            if( oRefSRS.importFromEPSG(atoi(pszAuthCode)) == OGRERR_NONE &&
-                oRefSRS.GetTOWGS84(adfRefTOWGS84) == OGRERR_NONE &&
-                memcmp(adfRefTOWGS84, adfTOWGS84, sizeof(adfTOWGS84)) == 0 )
+            if( oRefSRS.importFromEPSG(atoi(pszAuthCode)) == OGRERR_NONE )
             {
-                bUseReferenceTOWGS84 = true;
+                oRefSRS.AddGuessedTOWGS84();
+                if( oRefSRS.GetTOWGS84(adfRefTOWGS84) == OGRERR_NONE &&
+                    memcmp(adfRefTOWGS84, adfTOWGS84, sizeof(adfTOWGS84)) == 0 )
+                {
+                    bUseReferenceTOWGS84 = true;
+                }
             }
         }
         const char* pszWriteTOWGS84 =

diff --git a/gdal/frmts/gtiff/gt_wkt_srs.h b/gdal/frmts/gtiff/gt_wkt_srs.h
@@ -32,6 +32,7 @@
 #define GT_WKT_SRS_H_INCLUDED
 
 #include "cpl_port.h"
+#include "ogr_srs_api.h"
 
 #include "geo_normalize.h"
 #include "geotiff.h"
@@ -53,6 +54,8 @@ typedef enum
     GEOTIFF_VERSION_1_1
 } GeoTIFFVersionEnum;
 
+OGRSpatialReferenceH GTIFGetOGISDefnAsOSR( GTIF *, GTIFDefn * );
+
 int GTIFSetFromOGISDefnEx( GTIF *, const char *, GTIFFKeysFlavorEnum,
                            GeoTIFFVersionEnum );
 

diff --git a/gdal/ogr/ogr_spatialref.h b/gdal/ogr/ogr_spatialref.h
@@ -367,6 +367,7 @@ class CPL_DLL OGRSpatialReference
                             double = 0.0, double = 0.0, double = 0.0,
                             double = 0.0 );
     OGRErr      GetTOWGS84( double *padfCoef, int nCoeff = 7 ) const;
+    OGRErr      AddGuessedTOWGS84();
 
     double      GetSemiMajor( OGRErr * = nullptr ) const;
     double      GetSemiMinor( OGRErr * = nullptr ) const;

diff --git a/gdal/ogr/ogr_srs_api.h b/gdal/ogr/ogr_srs_api.h
@@ -569,6 +569,7 @@ OGRErr CPL_DLL OSRSetTOWGS84( OGRSpatialReferenceH hSRS,
                               double, double, double,
                               double, double, double, double );
 OGRErr CPL_DLL OSRGetTOWGS84( OGRSpatialReferenceH hSRS, double *, int );
+OGRErr CPL_DLL OSRAddGuessedTOWGS84( OGRSpatialReferenceH hSRS);
 
 OGRErr CPL_DLL OSRSetCompoundCS( OGRSpatialReferenceH hSRS,
                                  const char *pszName,

diff --git a/gdal/ogr/ogrct.cpp b/gdal/ogr/ogrct.cpp
@@ -841,10 +841,11 @@ int OGRProjCT::Initialize( const OGRSpatialReference * poSourceIn,
     else if( !bWebMercatorToWGS84LongLat )
     {
         const auto CanUseAuthorityDef = [](const OGRSpatialReference* poSRS1,
-                                           const OGRSpatialReference* poSRS2,
+                                           OGRSpatialReference* poSRSFromAuth,
                                            const char* pszAuth)
         {
-            if( EQUAL(pszAuth, "EPSG") )
+            if( EQUAL(pszAuth, "EPSG") &&
+                CPLTestBool(CPLGetConfigOption("OSR_CT_USE_DEFAULT_EPSG_TOWGS84", "NO")) )
             {
                 // We don't want by default to honour 'default' TOWGS84 terms that come with the EPSG code
                 // because there might be a better transformation from that
@@ -855,10 +856,12 @@ int OGRProjCT::Initialize( const OGRSpatialReference * poSourceIn,
                 // OSR_CT_USE_DEFAULT_EPSG_TOWGS84 configuration option to YES
                 double adfTOWGS84_1[7];
                 double adfTOWGS84_2[7];
+
+                poSRSFromAuth->AddGuessedTOWGS84();
+
                 if( poSRS1->GetTOWGS84(adfTOWGS84_1) == OGRERR_NONE &&
-                    poSRS2->GetTOWGS84(adfTOWGS84_2) == OGRERR_NONE &&
-                    memcmp(adfTOWGS84_1, adfTOWGS84_2, sizeof(adfTOWGS84_1)) == 0 &&
-                    CPLTestBool(CPLGetConfigOption("OSR_CT_USE_DEFAULT_EPSG_TOWGS84", "NO")) )
+                    poSRSFromAuth->GetTOWGS84(adfTOWGS84_2) == OGRERR_NONE &&
+                    memcmp(adfTOWGS84_1, adfTOWGS84_2, sizeof(adfTOWGS84_1)) == 0 )
                 {
                     return false;
                 }