Fix misalignment of rasters with RPC (fixes #35796) (#36385)

wonder-sk · web-flow · commit d48c17816f31 · 2020-05-20T17:38:01.000+02:00
* Fix misalignment of rasters with RPC (fixes #35796) 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. Fortunately we can pass RPC_HEIGHT to GDAL to use given value as the fixed elevation for the whole image. We try to use HEIGHT_OFF coefficient ("Geodetic Height Offset") as an estimate for elevation (it seems that ENVI software does that as well). In the future we may want to use also RPC_DEM for an even more precise georeferencing. https://gdal.org/development/rfc/rfc22_rpc.html * Update code to use the new fuction GDALAutoCreateWarpedVRTEx() This function will be available in GDAL >= 3.2 so we have a local copy for the time being for older versions of GDAL.
diff --git a/src/core/providers/gdal/qgsgdalprovider.cpp b/src/core/providers/gdal/qgsgdalprovider.cpp
@@ -2670,9 +2670,18 @@ void QgsGdalProvider::initBaseDataset()
   {
     QgsLogger::warning( QStringLiteral( "Creating Warped VRT." ) );
 
-    mGdalDataset =
-      GDALAutoCreateWarpedVRT( mGdalBaseDataset, nullptr, nullptr,
-                               GRA_NearestNeighbour, 0.2, nullptr );
+    if ( GDALGetMetadata( mGdalBaseDataset, "RPC" ) )
+    {
+      mGdalDataset =
+        QgsGdalUtils::rpcAwareAutoCreateWarpedVrt( mGdalBaseDataset, nullptr, nullptr,
+            GRA_NearestNeighbour, 0.2, nullptr );
+    }
+    else
+    {
+      mGdalDataset =
+        GDALAutoCreateWarpedVRT( mGdalBaseDataset, nullptr, nullptr,
+                                 GRA_NearestNeighbour, 0.2, nullptr );
+    }
 
     if ( !mGdalDataset )
     {
diff --git a/src/core/qgsgdalutils.cpp b/src/core/qgsgdalutils.cpp
@@ -18,6 +18,7 @@
 
 #define CPL_SUPRESS_CPLUSPLUS  //#spellok
 #include "gdal.h"
+#include "gdalwarper.h"
 #include "cpl_string.h"
 
 #include <QString>
@@ -277,3 +278,211 @@ QString QgsGdalUtils::validateCreationOptionsFormat( const QStringList &createOp
     return QStringLiteral( "Failed GDALValidateCreationOptions() test" );
   return QString();
 }
+
+#if GDAL_VERSION_NUM < GDAL_COMPUTE_VERSION(2,3,0)
+// GDAL < 2.3 does not come with GDALWarpInitDefaultBandMapping and GDALWarpInitNoDataReal
+// in the public API so we define them here for rpcAwareAutoCreateWarpedVrt()
+
+static void GDALWarpInitDefaultBandMapping( GDALWarpOptions *psOptionsIn, int nBandCount )
+{
+  if ( psOptionsIn->nBandCount != 0 ) { return; }
+
+  psOptionsIn->nBandCount = nBandCount;
+
+  psOptionsIn->panSrcBands = static_cast<int *>(
+                               CPLMalloc( sizeof( int ) * psOptionsIn->nBandCount ) );
+  psOptionsIn->panDstBands = static_cast<int *>(
+                               CPLMalloc( sizeof( int ) * psOptionsIn->nBandCount ) );
+
+  for ( int i = 0; i < psOptionsIn->nBandCount; i++ )
+  {
+    psOptionsIn->panSrcBands[i] = i + 1;
+    psOptionsIn->panDstBands[i] = i + 1;
+  }
+}
+
+static void InitNoData( int nBandCount, double **ppdNoDataReal, double dDataReal )
+{
+  if ( nBandCount <= 0 ) { return; }
+  if ( *ppdNoDataReal != nullptr ) { return; }
+
+  *ppdNoDataReal = static_cast<double *>(
+                     CPLMalloc( sizeof( double ) * nBandCount ) );
+
+  for ( int i = 0; i < nBandCount; ++i )
+  {
+    ( *ppdNoDataReal )[i] = dDataReal;
+  }
+}
+
+static void GDALWarpInitNoDataReal( GDALWarpOptions *psOptionsIn, double  dNoDataReal )
+{
+  InitNoData( psOptionsIn->nBandCount, &psOptionsIn->padfDstNoDataReal, dNoDataReal );
+  InitNoData( psOptionsIn->nBandCount, &psOptionsIn->padfSrcNoDataReal, dNoDataReal );
+  // older GDAL also requires imaginary values to be set
+  InitNoData( psOptionsIn->nBandCount, &psOptionsIn->padfDstNoDataImag, 0 );
+  InitNoData( psOptionsIn->nBandCount, &psOptionsIn->padfSrcNoDataImag, 0 );
+}
+#endif
+
+#if GDAL_VERSION_NUM < GDAL_COMPUTE_VERSION(3,2,0)
+
+GDALDatasetH GDALAutoCreateWarpedVRTEx( GDALDatasetH hSrcDS, const char *pszSrcWKT, const char *pszDstWKT, GDALResampleAlg eResampleAlg,
+                                        double dfMaxError, const GDALWarpOptions *psOptionsIn, char **papszTransformerOptions )
+{
+  VALIDATE_POINTER1( hSrcDS, "GDALAutoCreateWarpedVRT", nullptr );
+
+  /* -------------------------------------------------------------------- */
+  /*      Populate the warp options.                                      */
+  /* -------------------------------------------------------------------- */
+  GDALWarpOptions *psWO = nullptr;
+  if ( psOptionsIn != nullptr )
+    psWO = GDALCloneWarpOptions( psOptionsIn );
+  else
+    psWO = GDALCreateWarpOptions();
+
+  psWO->eResampleAlg = eResampleAlg;
+
+  psWO->hSrcDS = hSrcDS;
+
+  GDALWarpInitDefaultBandMapping( psWO, GDALGetRasterCount( hSrcDS ) );
+
+  /* -------------------------------------------------------------------- */
+  /*      Setup no data values                                            */
+  /* -------------------------------------------------------------------- */
+  for ( int i = 0; i < psWO->nBandCount; i++ )
+  {
+    GDALRasterBandH rasterBand = GDALGetRasterBand( psWO->hSrcDS, psWO->panSrcBands[i] );
+
+    int hasNoDataValue;
+    double noDataValue = GDALGetRasterNoDataValue( rasterBand, &hasNoDataValue );
+
+    if ( hasNoDataValue )
+    {
+      // Check if the nodata value is out of range
+      int bClamped = FALSE;
+      int bRounded = FALSE;
+      CPL_IGNORE_RET_VAL(
+        GDALAdjustValueToDataType( GDALGetRasterDataType( rasterBand ),
+                                   noDataValue, &bClamped, &bRounded ) );
+      if ( !bClamped )
+      {
+        GDALWarpInitNoDataReal( psWO, -1e10 );
+
+        psWO->padfSrcNoDataReal[i] = noDataValue;
+        psWO->padfDstNoDataReal[i] = noDataValue;
+      }
+    }
+  }
+
+  if ( psWO->padfDstNoDataReal != nullptr )
+  {
+    if ( CSLFetchNameValue( psWO->papszWarpOptions, "INIT_DEST" ) == nullptr )
+    {
+      psWO->papszWarpOptions =
+        CSLSetNameValue( psWO->papszWarpOptions, "INIT_DEST", "NO_DATA" );
+    }
+  }
+
+  /* -------------------------------------------------------------------- */
+  /*      Create the transformer.                                         */
+  /* -------------------------------------------------------------------- */
+  psWO->pfnTransformer = GDALGenImgProjTransform;
+
+  char **papszOptions = nullptr;
+  if ( pszSrcWKT != nullptr )
+    papszOptions = CSLSetNameValue( papszOptions, "SRC_SRS", pszSrcWKT );
+  if ( pszDstWKT != nullptr )
+    papszOptions = CSLSetNameValue( papszOptions, "DST_SRS", pszDstWKT );
+  papszOptions = CSLMerge( papszOptions, papszTransformerOptions );
+  psWO->pTransformerArg =
+    GDALCreateGenImgProjTransformer2( psWO->hSrcDS, nullptr,
+                                      papszOptions );
+  CSLDestroy( papszOptions );
+
+  if ( psWO->pTransformerArg == nullptr )
+  {
+    GDALDestroyWarpOptions( psWO );
+    return nullptr;
+  }
+
+  /* -------------------------------------------------------------------- */
+  /*      Figure out the desired output bounds and resolution.            */
+  /* -------------------------------------------------------------------- */
+  double adfDstGeoTransform[6] = { 0.0 };
+  int nDstPixels = 0;
+  int nDstLines = 0;
+  CPLErr eErr =
+    GDALSuggestedWarpOutput( hSrcDS, psWO->pfnTransformer,
+                             psWO->pTransformerArg,
+                             adfDstGeoTransform, &nDstPixels, &nDstLines );
+  if ( eErr != CE_None )
+  {
+    GDALDestroyTransformer( psWO->pTransformerArg );
+    GDALDestroyWarpOptions( psWO );
+    return nullptr;
+  }
+
+  /* -------------------------------------------------------------------- */
+  /*      Update the transformer to include an output geotransform        */
+  /*      back to pixel/line coordinates.                                 */
+  /*                                                                      */
+  /* -------------------------------------------------------------------- */
+  GDALSetGenImgProjTransformerDstGeoTransform(
+    psWO->pTransformerArg, adfDstGeoTransform );
+
+  /* -------------------------------------------------------------------- */
+  /*      Do we want to apply an approximating transformation?            */
+  /* -------------------------------------------------------------------- */
+  if ( dfMaxError > 0.0 )
+  {
+    psWO->pTransformerArg =
+      GDALCreateApproxTransformer( psWO->pfnTransformer,
+                                   psWO->pTransformerArg,
+                                   dfMaxError );
+    psWO->pfnTransformer = GDALApproxTransform;
+    GDALApproxTransformerOwnsSubtransformer( psWO->pTransformerArg, TRUE );
+  }
+
+  /* -------------------------------------------------------------------- */
+  /*      Create the VRT file.                                            */
+  /* -------------------------------------------------------------------- */
+  GDALDatasetH hDstDS
+    = GDALCreateWarpedVRT( hSrcDS, nDstPixels, nDstLines,
+                           adfDstGeoTransform, psWO );
+
+  GDALDestroyWarpOptions( psWO );
+
+  if ( pszDstWKT != nullptr )
+    GDALSetProjection( hDstDS, pszDstWKT );
+  else if ( pszSrcWKT != nullptr )
+    GDALSetProjection( hDstDS, pszSrcWKT );
+  else if ( GDALGetGCPCount( hSrcDS ) > 0 )
+    GDALSetProjection( hDstDS, GDALGetGCPProjection( hSrcDS ) );
+  else
+    GDALSetProjection( hDstDS, GDALGetProjectionRef( hSrcDS ) );
+
+  return hDstDS;
+}
+#endif
+
+
+GDALDatasetH QgsGdalUtils::rpcAwareAutoCreateWarpedVrt(
+  GDALDatasetH hSrcDS,
+  const char *pszSrcWKT,
+  const char *pszDstWKT,
+  GDALResampleAlg eResampleAlg,
+  double dfMaxError,
+  const GDALWarpOptions *psOptionsIn )
+{
+  char **opts = nullptr;
+  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 );
+  }
+
+  return GDALAutoCreateWarpedVRTEx( hSrcDS, pszSrcWKT, pszDstWKT, eResampleAlg, dfMaxError, psOptionsIn, opts );
+}
diff --git a/src/core/qgsgdalutils.h b/src/core/qgsgdalutils.h
@@ -100,6 +100,22 @@ class CORE_EXPORT QgsGdalUtils
      */
     static gdal::dataset_unique_ptr imageToMemoryDataset( const QImage &image );
 
+    /**
+     * This is a copy of GDALAutoCreateWarpedVRT optimized for imagery using RPC georeferencing
+     * that also 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.14
+     */
+    static GDALDatasetH rpcAwareAutoCreateWarpedVrt(
+      GDALDatasetH hSrcDS,
+      const char *pszSrcWKT,
+      const char *pszDstWKT,
+      GDALResampleAlg eResampleAlg,
+      double dfMaxError,
+      const GDALWarpOptions *psOptionsIn );
+
     friend class TestQgsGdalUtils;
 };
 
