18#include <gdalwarper.h>
19#include <ogr_srs_api.h>
30static double ceil_with_tolerance(
double value )
32 if ( std::fabs( value - std::round( value ) ) < 1e-6 )
33 return std::round( value );
35 return std::ceil( value );
38static double floor_with_tolerance(
double value )
40 if ( std::fabs( value - std::round( value ) ) < 1e-6 )
41 return std::round( value );
43 return std::floor( value );
46static double fmod_with_tolerance(
double num,
double denom )
48 return num - floor_with_tolerance( num / denom ) * denom;
51static QgsRectangle transform_to_extent(
const double *geotransform,
double xSize,
double ySize )
55 geotransform[0] + geotransform[1] * xSize,
56 geotransform[3] + geotransform[5] * ySize );
61static int CPL_STDCALL _progress(
double dfComplete,
const char *pszMessage,
void *pProgressArg )
63 Q_UNUSED( pszMessage )
67 return handler->
progress( dfComplete );
72static CPLErr rescalePreWarpChunkProcessor(
void *pKern,
void *pArg )
74 GDALWarpKernel *kern = ( GDALWarpKernel * ) pKern;
75 const double cellsize = ( (
double * )pArg )[0];
77 for (
int nBand = 0; nBand < kern->nBands; ++nBand )
79 float *bandData = (
float * ) kern->papabySrcImage[nBand];
80 for (
int nLine = 0; nLine < kern->nSrcYSize; ++nLine )
82 for (
int nPixel = 0; nPixel < kern->nSrcXSize; ++nPixel )
84 bandData[nPixel] /= cellsize;
86 bandData += kern->nSrcXSize;
92static CPLErr rescalePostWarpChunkProcessor(
void *pKern,
void *pArg )
94 GDALWarpKernel *kern = ( GDALWarpKernel * ) pKern;
95 const double cellsize = ( (
double * )pArg )[1];
97 for (
int nBand = 0; nBand < kern->nBands; ++nBand )
99 float *bandData = (
float * ) kern->papabyDstImage[nBand];
100 for (
int nLine = 0; nLine < kern->nDstYSize; ++nLine )
102 for (
int nPixel = 0; nPixel < kern->nDstXSize; ++nPixel )
104 bandData[nPixel] *= cellsize;
106 bandData += kern->nDstXSize;
123 for (
int i = 0; i < 6; ++i )
155 if ( customCRSWkt.isEmpty() || customCRSWkt == rasterInfo.
crs() )
160 if ( !customCellSize.isValid() )
171 if ( customGridOffset.x() < 0 || customGridOffset.y() < 0 )
173 if ( !customCellSize.isValid() )
183 mGridOffsetX = fmod_with_tolerance( rasterInfo.
origin().x(), customCellSize.width() );
184 mGridOffsetY = fmod_with_tolerance( rasterInfo.
origin().y(), customCellSize.height() );
199 mCrsWkt = QStringLiteral(
"_error_" );
207 if ( !customCellSize.isValid() )
218 if ( customGridOffset.x() < 0 || customGridOffset.y() < 0 )
237 if (
mCrsWkt == QLatin1String(
"_error_" ) )
245 mErrorMessage = QObject::tr(
"Cell size must not be zero." );
252 double finalExtent[4] = { 0, 0, 0, 0 };
255 for (
int i = 0; i <
mRasters.count(); ++i )
265 mErrorMessage = QString(
"Failed to get suggested warp output.\n\n"
267 "Source WKT:\n%2\n\nDestination WKT:\n%3" )
276 if ( finalExtent[0] == 0 && finalExtent[1] == 0 && finalExtent[2] == 0 && finalExtent[3] == 0 )
287 if ( extent.
xMinimum() > finalExtent[0] ) finalExtent[0] = extent.
xMinimum();
288 if ( extent.
yMinimum() > finalExtent[1] ) finalExtent[1] = extent.
yMinimum();
289 if ( extent.
xMaximum() < finalExtent[2] ) finalExtent[2] = extent.
xMaximum();
290 if ( extent.
yMaximum() < finalExtent[3] ) finalExtent[3] = extent.
yMaximum();
305 if ( clipX0 > finalExtent[0] ) finalExtent[0] = clipX0;
306 if ( clipY0 > finalExtent[1] ) finalExtent[1] = clipY0;
307 if ( clipX1 < finalExtent[2] ) finalExtent[2] = clipX1;
308 if ( clipY1 < finalExtent[3] ) finalExtent[3] = clipY1;
317 const int xSize = floor_with_tolerance( ( finalExtent[2] - originX ) /
mCellSizeX );
318 const int ySize = floor_with_tolerance( ( finalExtent[3] - originY ) /
mCellSizeY );
320 if ( xSize <= 0 || ySize <= 0 )
322 mErrorMessage = QObject::tr(
"No common intersecting area." );
362 const auto constMRasters =
mRasters;
363 for (
const Item &r : constMRasters )
374 qDebug(
"---ALIGN------------------" );
375 qDebug(
"wkt %s",
mCrsWkt.toLatin1().constData() );
377 qDebug(
"transform" );
382 qDebug(
"extent %s", e.
toString().toLatin1().constData() );
388 double bestCellArea = INFINITY;
398 const auto constMRasters =
mRasters;
399 for (
const Item &raster : constMRasters )
404 const double cellArea = cs.width() * cs.height();
405 if ( cellArea < bestCellArea )
407 bestCellArea = cellArea;
419 GDALDriverH hDriver = GDALGetDriverByName(
"GTiff" );
422 mErrorMessage = QStringLiteral(
"GDALGetDriverByName(GTiff) failed." );
436 const int bandCount = GDALGetRasterCount( hSrcDS.get() );
437 const GDALDataType eDT = GDALGetRasterDataType( GDALGetRasterBand( hSrcDS.get(), 1 ) );
441 bandCount, eDT,
nullptr ) );
449 GDALSetProjection( hDstDS.get(),
mCrsWkt.toLatin1().constData() );
453 GDALColorTableH hCT = GDALGetRasterColorTable( GDALGetRasterBand( hSrcDS.get(), 1 ) );
455 GDALSetRasterColorTable( GDALGetRasterBand( hDstDS.get(), 1 ), hCT );
461 psWarpOptions->hSrcDS = hSrcDS.get();
462 psWarpOptions->hDstDS = hDstDS.get();
464 psWarpOptions->nBandCount = GDALGetRasterCount( hSrcDS.get() );
465 psWarpOptions->panSrcBands = (
int * ) CPLMalloc(
sizeof(
int ) * psWarpOptions->nBandCount );
466 psWarpOptions->panDstBands = (
int * ) CPLMalloc(
sizeof(
int ) * psWarpOptions->nBandCount );
467 for (
int i = 0; i < psWarpOptions->nBandCount; ++i )
469 psWarpOptions->panSrcBands[i] = i + 1;
470 psWarpOptions->panDstBands[i] = i + 1;
473 psWarpOptions->eResampleAlg =
static_cast< GDALResampleAlg
>( raster.
resampleMethod );
476 psWarpOptions->pfnProgress = _progress;
477 psWarpOptions->pProgressArg =
this;
480 psWarpOptions->pTransformerArg =
481 GDALCreateGenImgProjTransformer( hSrcDS.get(), GDALGetProjectionRef( hSrcDS.get() ),
482 hDstDS.get(), GDALGetProjectionRef( hDstDS.get() ),
484 psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
486 double rescaleArg[2];
491 psWarpOptions->pfnPreWarpChunkProcessor = rescalePreWarpChunkProcessor;
492 psWarpOptions->pfnPostWarpChunkProcessor = rescalePostWarpChunkProcessor;
493 psWarpOptions->pPreWarpProcessorArg = rescaleArg;
494 psWarpOptions->pPostWarpProcessorArg = rescaleArg;
496 psWarpOptions->eWorkingDataType = GDT_Float32;
500 GDALWarpOperation oOperation;
501 oOperation.Initialize( psWarpOptions.get() );
504 GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
514 void *hTransformArg = GDALCreateGenImgProjTransformer( info.
mDataset.get(), info.
mCrsWkt.toLatin1().constData(),
nullptr, destWkt.toLatin1().constData(), FALSE, 0, 1 );
515 if ( !hTransformArg )
519 double adfDstGeoTransform[6];
521 int nPixels = 0, nLines = 0;
523 eErr = GDALSuggestedWarpOutput2( info.
mDataset.get(),
524 GDALGenImgProjTransform, hTransformArg,
525 adfDstGeoTransform, &nPixels, &nLines, extents, 0 );
526 GDALDestroyGenImgProjTransformer( hTransformArg );
528 if ( eErr != CE_None )
531 const QSizeF cs( std::fabs( adfDstGeoTransform[1] ), std::fabs( adfDstGeoTransform[5] ) );
534 *rect =
QgsRectangle( extents[0], extents[1], extents[2], extents[3] );
538 *
gridOffset = QPointF( fmod_with_tolerance( adfDstGeoTransform[0], cs.width() ),
539 fmod_with_tolerance( adfDstGeoTransform[3], cs.height() ) );
549 mDataset.reset( GDALOpen( layerpath.toUtf8().constData(), GA_ReadOnly ) );
559 mCrsWkt = QString::fromLatin1( GDALGetProjectionRef(
mDataset.get() ) );
587 qDebug(
"---RASTER INFO------------------" );
588 qDebug(
"wkt %s",
mCrsWkt.toLatin1().constData() );
593 qDebug(
"extent %s", r.
toString().toLatin1().constData() );
595 qDebug(
"transform" );
602 GDALRasterBandH hBand = GDALGetRasterBand( mDataset.get(), 1 );
608 float *pafScanline = (
float * ) CPLMalloc(
sizeof(
float ) );
609 const CPLErr err = GDALRasterIO( hBand, GF_Read, px, py, 1, 1,
610 pafScanline, 1, 1, GDT_Float32, 0, 0 );
611 const double value = err == CE_None ? pafScanline[0] : std::numeric_limits<double>::quiet_NaN();
612 CPLFree( pafScanline );
@ PreferredGdal
Preferred format for conversion of CRS to WKT for use with the GDAL library.
QgsAlignRaster takes one or more raster layers and warps (resamples) them so they have the same:
void dump() const
write contents of the object to standard error stream - for debugging
int mYSize
Computed raster grid height.
bool setParametersFromRaster(const RasterInfo &rasterInfo, const QString &customCRSWkt=QString(), QSizeF customCellSize=QSizeF(), QPointF customGridOffset=QPointF(-1, -1))
Set destination CRS, cell size and grid offset from a raster file.
QSizeF cellSize() const
Gets output cell size.
QgsRectangle clipExtent() const
Gets clipping extent (region of interest).
QPointF gridOffset() const
double mGridOffsetX
Destination grid offset - expected to be in interval <0,cellsize)
QgsRectangle alignedRasterExtent() const
Returns the expected extent of the resulting aligned raster.
bool run()
Run the alignment process.
double mClipExtent[4]
Optional clip extent: sets "requested area" which be extended to fit the raster grid.
List mRasters
List of rasters to be aligned (with their output files and other options)
bool createAndWarp(const Item &raster)
Internal function for processing of one raster (1. create output, 2. do the alignment)
QString mCrsWkt
Destination CRS - stored in well-known text (WKT) format.
void setClipExtent(double xmin, double ymin, double xmax, double ymax)
Configure clipping extent (region of interest).
int mXSize
Computed raster grid width.
int suggestedReferenceLayer() const
Returns the index of the layer which has smallest cell size (returns -1 on error)
bool checkInputParameters()
Determine destination extent from the input rasters and calculate derived values.
double mCellSizeX
Destination cell size.
double mGeoTransform[6]
Computed geo-transform.
static bool suggestedWarpOutput(const RasterInfo &info, const QString &destWkt, QSizeF *cellSize=nullptr, QPointF *gridOffset=nullptr, QgsRectangle *rect=nullptr)
Determine suggested output of raster warp to a different CRS. Returns true on success.
QSize alignedRasterSize() const
Returns the expected size of the resulting aligned raster.
QString mErrorMessage
Last error message from run()
This class represents a coordinate reference system (CRS).
QString toWkt(Qgis::CrsWktVariant variant=Qgis::CrsWktVariant::Wkt1Gdal, bool multiline=false, int indentationWidth=4) const
Returns a WKT representation of this CRS.
A rectangle specified with double values.
QString toString(int precision=16) const
Returns a string representation of form xmin,ymin : xmax,ymax Coordinates will be truncated to the sp...
double xMinimum() const
Returns the x minimum value (left side of rectangle).
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
double xMaximum() const
Returns the x maximum value (right side of rectangle).
double yMaximum() const
Returns the y maximum value (top side of rectangle).
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
std::unique_ptr< GDALWarpOptions, GDALWarpOptionsDeleter > warp_options_unique_ptr
Scoped GDAL warp options.
Definition of one raster layer for alignment.
Qgis::GdalResampleAlgorithm resampleMethod
resampling method to be used
bool rescaleValues
rescaling of values according to the change of pixel size
double srcCellSizeInDestCRS
used for rescaling of values (if necessary)
QString inputFilename
filename of the source raster
QString outputFilename
filename of the newly created aligned raster (will be overwritten if exists already)
Helper struct to be sub-classed for progress reporting.
virtual bool progress(double complete)=0
Method to be overridden for progress reporting.
Utility class for gathering information about rasters.
QString crs() const
Returns the CRS in WKT format.
int mYSize
raster grid size Y
gdal::dataset_unique_ptr mDataset
handle to open GDAL dataset
void dump() const
Write contents of the object to standard error stream - for debugging.
RasterInfo(const QString &layerpath)
Construct raster info with a path to a raster file.
double mGeoTransform[6]
geotransform coefficients
QgsRectangle extent() const
Returns the extent of the raster.
QPointF gridOffset() const
Returns the grid offset.
QSizeF cellSize() const
Returns the cell size in map units.
double identify(double mx, double my)
Gets raster value at the given coordinates (from the first band)
int mXSize
raster grid size X
QPointF origin() const
Returns the origin of the raster.
int mBandCnt
number of raster's bands
QString mCrsWkt
CRS stored in WKT format.