Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
no GDAL resampling with unknown data type (#53200) 
* no GDAL resampling with unknown data type

* Update src/core/qgsgdalutils.cpp

Co-authored-by: Even Rouault <even.rouault@spatialys.com>

* Update src/core/raster/qgsrasterlayerrenderer.cpp

Co-authored-by: Nyall Dawson <nyall.dawson@gmail.com>

---------

Co-authored-by: Even Rouault <even.rouault@spatialys.com>
Co-authored-by: Nyall Dawson <nyall.dawson@gmail.com>
  • Loading branch information
@vcloarec @rouault @nyalldawson
3 people committed May 24, 2023
1 parent 01c10ae commit eab811f
Show file tree
Hide file tree
Showing 2 changed files with 45 additions and 41 deletions.
7 changes: 4 additions & 3 deletions src/core/qgsgdalutils.cpp
View file
Expand Up @@ -272,10 +272,11 @@ static bool resampleSingleBandRasterStatic( GDALDatasetH hSrcDS, GDALDatasetH hD
psWarpOptions->papszWarpOptions = CSLSetNameValue( psWarpOptions-> papszWarpOptions, "INIT_DEST", "NO_DATA" );

// Initialize and execute the warp operation.
bool retVal = false;
GDALWarpOperation oOperation;
oOperation.Initialize( psWarpOptions.get() );

const bool retVal { oOperation.ChunkAndWarpImage( 0, 0, GDALGetRasterXSize( hDstDS ), GDALGetRasterYSize( hDstDS ) ) == CE_None };
CPLErr initResult = oOperation.Initialize( psWarpOptions.get() );
if ( initResult != CE_Failure )
retVal = oOperation.ChunkAndWarpImage( 0, 0, GDALGetRasterXSize( hDstDS ), GDALGetRasterYSize( hDstDS ) ) == CE_None;
GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
return retVal;
}
Expand Down
79 changes: 41 additions & 38 deletions src/core/raster/qgsrasterlayerrenderer.cpp
View file
Expand Up @@ -483,45 +483,48 @@ void QgsRasterLayerRenderer::drawElevationMap()

Qgis::DataType dataType = dataProvider->dataType( mElevationBand );

// we need extra pixels on border to avoid effect border with resampling (at least 2 pixels band for cubic alg)
int sourceWidth = viewWidth + 4;
int sourceHeight = viewHeight + 4;
viewExtentInLayerCoordinate = QgsRectangle(
viewExtentInLayerCoordinate.xMinimum() - xLayerResol * 2,
viewExtentInLayerCoordinate.yMinimum() - yLayerResol * 2,
viewExtentInLayerCoordinate.xMaximum() + xLayerResol * 2,
viewExtentInLayerCoordinate.yMaximum() + yLayerResol * 2 );

// Now we can do the resampling
std::unique_ptr<QgsRasterBlock> sourcedata( dataProvider->block( mElevationBand, viewExtentInLayerCoordinate, sourceWidth, sourceHeight, mFeedback ) );
gdal::dataset_unique_ptr gdalDsInput =
QgsGdalUtils::blockToSingleBandMemoryDataset( viewExtentInLayerCoordinate, sourcedata.get() );


elevationBlock.reset( new QgsRasterBlock( dataType,
outputWidth,
outputHeight ) );

elevationBlock->setNoDataValue( dataProvider->sourceNoDataValue( mElevationBand ) );

gdal::dataset_unique_ptr gdalDsOutput =
QgsGdalUtils::blockToSingleBandMemoryDataset( mRasterViewPort->mDrawnExtent, elevationBlock.get() );

// For coordinate transformation, we try to obtain a coordinate operation string from the transform context.
// Depending of the CRS, if we can't we use GDAL transformation directly from the source and destination CRS
QString coordinateOperation;
const QgsCoordinateTransformContext &transformContext = renderContext()->transformContext();
if ( transformContext.mustReverseCoordinateOperation( mRasterViewPort->mDestCRS, mRasterViewPort->mSrcCRS ) )
coordinateOperation = transformContext.calculateCoordinateOperation( mRasterViewPort->mSrcCRS, mRasterViewPort->mDestCRS );
else
coordinateOperation = transformContext.calculateCoordinateOperation( mRasterViewPort->mDestCRS, mRasterViewPort->mSrcCRS );
if ( dataType != Qgis::DataType::UnknownDataType ) // resampling data by GDAL is not compatible with unknown data type
{
// we need extra pixels on border to avoid effect border with resampling (at least 2 pixels band for cubic alg)
int sourceWidth = viewWidth + 4;
int sourceHeight = viewHeight + 4;
viewExtentInLayerCoordinate = QgsRectangle(
viewExtentInLayerCoordinate.xMinimum() - xLayerResol * 2,
viewExtentInLayerCoordinate.yMinimum() - yLayerResol * 2,
viewExtentInLayerCoordinate.xMaximum() + xLayerResol * 2,
viewExtentInLayerCoordinate.yMaximum() + yLayerResol * 2 );

// Now we can do the resampling
std::unique_ptr<QgsRasterBlock> sourcedata( dataProvider->block( mElevationBand, viewExtentInLayerCoordinate, sourceWidth, sourceHeight, mFeedback ) );
gdal::dataset_unique_ptr gdalDsInput =
QgsGdalUtils::blockToSingleBandMemoryDataset( viewExtentInLayerCoordinate, sourcedata.get() );

if ( coordinateOperation.isEmpty() )
canRenderElevation = QgsGdalUtils::resampleSingleBandRaster( gdalDsInput.get(), gdalDsOutput.get(), alg,
mRasterViewPort->mSrcCRS, mRasterViewPort->mDestCRS );
else
canRenderElevation = QgsGdalUtils::resampleSingleBandRaster( gdalDsInput.get(), gdalDsOutput.get(), alg,
coordinateOperation.toUtf8().constData() );

elevationBlock.reset( new QgsRasterBlock( dataType,
outputWidth,
outputHeight ) );

elevationBlock->setNoDataValue( dataProvider->sourceNoDataValue( mElevationBand ) );

gdal::dataset_unique_ptr gdalDsOutput =
QgsGdalUtils::blockToSingleBandMemoryDataset( mRasterViewPort->mDrawnExtent, elevationBlock.get() );

// For coordinate transformation, we try to obtain a coordinate operation string from the transform context.
// Depending of the CRS, if we can't we use GDAL transformation directly from the source and destination CRS
QString coordinateOperation;
const QgsCoordinateTransformContext &transformContext = renderContext()->transformContext();
if ( transformContext.mustReverseCoordinateOperation( mRasterViewPort->mDestCRS, mRasterViewPort->mSrcCRS ) )
coordinateOperation = transformContext.calculateCoordinateOperation( mRasterViewPort->mSrcCRS, mRasterViewPort->mDestCRS );
else
coordinateOperation = transformContext.calculateCoordinateOperation( mRasterViewPort->mDestCRS, mRasterViewPort->mSrcCRS );

if ( coordinateOperation.isEmpty() )
canRenderElevation = QgsGdalUtils::resampleSingleBandRaster( gdalDsInput.get(), gdalDsOutput.get(), alg,
mRasterViewPort->mSrcCRS, mRasterViewPort->mDestCRS );
else
canRenderElevation = QgsGdalUtils::resampleSingleBandRaster( gdalDsInput.get(), gdalDsOutput.get(), alg,
coordinateOperation.toUtf8().constData() );
}
}

if ( canRenderElevation )
Expand Down

0 comments on commit eab811f

Please sign in to comment.