mirror of
https://github.com/qgis/QGIS.git
synced 2025-10-17 00:09:36 -04: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.
This commit is contained in:
Martin Dobias
GitHub
parent
7053b78ac6
commit
d48c17816f
@ -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 )
|
||||
{
|
||||
|
||||
@ -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 );
|
||||
}
|
||||
|
||||
@ -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;
|
||||
};
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user