Conversion between QgsCoordinateReferenceSystem and OGRSpatialReferen…

…ce: use authId as much as possible This helps preserving metadata, like datum code, which is generally absent from WKT2 output. Or when importing from WKT1, extent is missing. This is an alternative to the fix of OSGeo/gdal#5218 for OSGeo/gdal#5217.
qgis · Feb 3, 2022 · 7bd931c · 7bd931c
1 parent 71ab6ec
commit 7bd931c
Show file tree

Hide file tree

Showing 2 changed files with 76 additions and 6 deletions.
diff --git a/src/core/qgsogrutils.cpp b/src/core/qgsogrutils.cpp
@@ -962,24 +962,76 @@ QgsCoordinateReferenceSystem QgsOgrUtils::OGRSpatialReferenceToCrs( OGRSpatialRe
   if ( wkt.isEmpty() )
     return QgsCoordinateReferenceSystem();
 
+  const char *authorityName = OSRGetAuthorityName( srs, nullptr );
+  const char *authorityCode = OSRGetAuthorityCode( srs, nullptr );
+  QgsCoordinateReferenceSystem res;
+  if ( authorityName && authorityCode )
+  {
+    QString authId = QString( authorityName ) + ':' + QString( authorityCode );
+    OGRSpatialReferenceH ogrSrsTmp = OSRNewSpatialReference( nullptr );
+    // Check that the CRS build from authId and the input one are the "same".
+    if ( OSRSetFromUserInput( ogrSrsTmp, authId.toUtf8().constData() ) != OGRERR_NONE &&
+         OSRIsSame( srs, ogrSrsTmp ) )
+    {
+      res = QgsCoordinateReferenceSystem();
+      res.createFromUserInput( authId );
+    }
+    OSRDestroySpatialReference( ogrSrsTmp );
+  }
+  if ( !res.isValid() )
+    res = QgsCoordinateReferenceSystem::fromWkt( wkt );
+
 #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3,4,0)
-  QgsCoordinateReferenceSystem res = QgsCoordinateReferenceSystem::fromWkt( wkt );
   const double coordinateEpoch = OSRGetCoordinateEpoch( srs );
   if ( coordinateEpoch > 0 )
     res.setCoordinateEpoch( coordinateEpoch );
-  return res;
-#else
-  return QgsCoordinateReferenceSystem::fromWkt( wkt );
 #endif
+  return res;
 }
 
 OGRSpatialReferenceH QgsOgrUtils::crsToOGRSpatialReference( const QgsCoordinateReferenceSystem &crs )
 {
   if ( crs.isValid() )
   {
-    const QString srsWkt = crs.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED_GDAL );
+    OGRSpatialReferenceH ogrSrs = nullptr;
 
-    if ( OGRSpatialReferenceH ogrSrs = OSRNewSpatialReference( srsWkt.toLocal8Bit().constData() ) )
+    // First try instantiating the CRS from its authId. This will give a
+    // more complete representation of the CRS for GDAL. In particular it might
+    // help a few drivers to get the datum code, that would be missing in WKT-2.
+    // See https://github.com/OSGeo/gdal/pull/5218
+    const QString authId = crs.authid();
+    const QString srsWkt = crs.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED_GDAL );
+    if ( !authId.isEmpty() )
+    {
+      ogrSrs = OSRNewSpatialReference( nullptr );
+      if ( OSRSetFromUserInput( ogrSrs, authId.toUtf8().constData() ) == OGRERR_NONE )
+      {
+        // Check that the CRS build from WKT and authId are the "same".
+        OGRSpatialReferenceH ogrSrsFromWkt = OSRNewSpatialReference( srsWkt.toUtf8().constData() );
+        if ( ogrSrsFromWkt )
+        {
+          if ( !OSRIsSame( ogrSrs, ogrSrsFromWkt ) )
+          {
+            OSRDestroySpatialReference( ogrSrs );
+            ogrSrs = ogrSrsFromWkt;
+          }
+          else
+          {
+            OSRDestroySpatialReference( ogrSrsFromWkt );
+          }
+        }
+      }
+      else
+      {
+        OSRDestroySpatialReference( ogrSrs );
+        ogrSrs = nullptr;
+      }
+    }
+    if ( !ogrSrs )
+    {
+      ogrSrs = OSRNewSpatialReference( srsWkt.toUtf8().constData() );
+    }
+    if ( ogrSrs )
     {
       OSRSetAxisMappingStrategy( ogrSrs, OAMS_TRADITIONAL_GIS_ORDER );
 #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3,4,0)

diff --git a/tests/src/core/testqgsogrutils.cpp b/tests/src/core/testqgsogrutils.cpp
@@ -727,6 +727,15 @@ void TestQgsOgrUtils::ogrCrsConversion()
     const QgsCoordinateReferenceSystem crs1( QStringLiteral( "EPSG:3111" ) );
     OGRSpatialReferenceH srs = QgsOgrUtils::crsToOGRSpatialReference( crs1 );
     QVERIFY( srs );
+
+    // Check that OGRSpatialReferenceH object built has all information preserved
+    const char *authName = OSRGetAuthorityName( srs, "DATUM" );
+    QVERIFY( authName );
+    QCOMPARE( QString( authName ), "EPSG" );
+    const char *authCode = OSRGetAuthorityCode( srs, "DATUM" );
+    QVERIFY( authCode );
+    QCOMPARE( QString( authCode ), "6283" );
+
     const QgsCoordinateReferenceSystem crs2( QgsOgrUtils::OGRSpatialReferenceToCrs( srs ) );
     // round trip should be lossless
     QCOMPARE( crs1, crs2 );
@@ -737,6 +746,15 @@ void TestQgsOgrUtils::ogrCrsConversion()
 #endif
   }
 
+  {
+    OGRSpatialReferenceH srs = OSRNewSpatialReference( "PROJCS[\"GDA94 / Vicgrid\",GEOGCS[\"GDA94\",DATUM[\"Geocentric_Datum_of_Australia_1994\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6283\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4283\"]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"latitude_of_origin\",-37],PARAMETER[\"central_meridian\",145],PARAMETER[\"standard_parallel_1\",-36],PARAMETER[\"standard_parallel_2\",-38],PARAMETER[\"false_easting\",2500000],PARAMETER[\"false_northing\",2500000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"3111\"]]" );
+    // Check that we used EPSG:3111 to instantiate the CRS, and thus get the
+    // extent from PROJ
+    const QgsCoordinateReferenceSystem crs( QgsOgrUtils::OGRSpatialReferenceToCrs( srs ) );
+    OSRRelease( srs );
+    QVERIFY( !crs.bounds().isEmpty() );
+  }
+
 #if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3,4,0)
   {
     // test conversion with a coordinate epoch, should work on GDAL 3.4+