Fix transformCoordinates() for rasters with RPC georeferencing

qgis · Sep 2, 2020 · d750708 · d750708
1 parent e60e65c
commit d750708
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
@@ -2680,6 +2680,7 @@ void QgsGdalProvider::initBaseDataset()
       mGdalDataset =
         QgsGdalUtils::rpcAwareAutoCreateWarpedVrt( mGdalBaseDataset, nullptr, nullptr,
             GRA_NearestNeighbour, 0.2, nullptr );
+      mGdalTransformerArg = QgsGdalUtils::rpcAwareCreateTransformer( mGdalBaseDataset );
     }
     else
     {
@@ -2703,7 +2704,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();