Fix transformCoordinates() for rasters with RPC georeferencing (#38350)

Unfortunately for raster with RPC georeferencing, the transform used in GDAL provider's transformCoordinates() was using RPC_HEIGHT at zero elevation, while rendering pipeline was using RPC_HEIGHT at the correct elevation from HEIGHT_OFF metadata. This was causing shift in what was rendered in map canvas vs what transform function returned. This is now fixed by also the transform function using the correct elevation.
qgis · Aug 21, 2020 · 418bcf3 · 418bcf3
1 parent b98cfe1
commit 418bcf3
Show file tree

Hide file tree

Showing 3 changed files with 30 additions and 1 deletion.
diff --git a/src/core/providers/gdal/qgsgdalprovider.cpp b/src/core/providers/gdal/qgsgdalprovider.cpp
@@ -2920,6 +2920,7 @@ void QgsGdalProvider::initBaseDataset()
       mGdalDataset =
         QgsGdalUtils::rpcAwareAutoCreateWarpedVrt( mGdalBaseDataset, nullptr, nullptr,
             GRA_NearestNeighbour, 0.2, nullptr );
+      mGdalTransformerArg = QgsGdalUtils::rpcAwareCreateTransformer( mGdalBaseDataset );
     }
     else
     {
@@ -2943,7 +2944,10 @@ void QgsGdalProvider::initBaseDataset()
     mGdalDataset = mGdalBaseDataset;
   }
 
-  mGdalTransformerArg = GDALCreateGenImgProjTransformer( mGdalBaseDataset, nullptr, nullptr, nullptr, TRUE, 1.0, 0 );
+  if ( !mGdalTransformerArg )
+  {
+    mGdalTransformerArg = GDALCreateGenImgProjTransformer( mGdalBaseDataset, nullptr, nullptr, nullptr, TRUE, 1.0, 0 );
+  }
 
   if ( !hasGeoTransform )
   {

diff --git a/src/core/qgsgdalutils.cpp b/src/core/qgsgdalutils.cpp
@@ -490,6 +490,21 @@ GDALDatasetH QgsGdalUtils::rpcAwareAutoCreateWarpedVrt(
   return GDALAutoCreateWarpedVRTEx( hSrcDS, pszSrcWKT, pszDstWKT, eResampleAlg, dfMaxError, psOptionsIn, opts );
 }
 
+void *QgsGdalUtils::rpcAwareCreateTransformer( GDALDatasetH hSrcDS, GDALDatasetH hDstDS, char **papszOptions )
+{
+  char **opts = CSLDuplicate( papszOptions );
+  if ( GDALGetMetadata( hSrcDS, "RPC" ) )
+  {
+    // well-behaved RPC should have height offset a good value for RPC_HEIGHT
+    const char *heightOffStr = GDALGetMetadataItem( hSrcDS, "HEIGHT_OFF", "RPC" );
+    if ( heightOffStr )
+      opts = CSLAddNameValue( opts, "RPC_HEIGHT", heightOffStr );
+  }
+  void *transformer = GDALCreateGenImgProjTransformer2( hSrcDS, hDstDS, opts );
+  CSLDestroy( opts );
+  return transformer;
+}
+
 #ifndef QT_NO_NETWORKPROXY
 void QgsGdalUtils::setupProxy()
 {

diff --git a/src/core/qgsgdalutils.h b/src/core/qgsgdalutils.h
@@ -116,6 +116,16 @@ class CORE_EXPORT QgsGdalUtils
       double dfMaxError,
       const GDALWarpOptions *psOptionsIn );
 
+    /**
+     * This is a wrapper around GDALCreateGenImgProjTransformer2() that takes into account RPC
+     * georeferencing (it sets RPC_HEIGHT in GDALCreateGenImgProjTransformer2 based on HEIGHT_OFF).
+     * By default GDAL would assume that the imagery has zero elevation - if that is not the case,
+     * the image would not be shown in the correct location.
+     *
+     * \since QGIS 3.16
+     */
+    static void *rpcAwareCreateTransformer( GDALDatasetH hSrcDS, GDALDatasetH hDstDS = nullptr, char **papszOptions = nullptr );
+
 #ifndef QT_NO_NETWORKPROXY
     //! Sets the gdal proxy variables
     static void setupProxy();