18 #include <gdalwarper.h> 19 #include <ogr_srs_api.h> 30 static 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 );
38 static 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 );
46 static double fmod_with_tolerance(
double num,
double denom )
48 return num - floor_with_tolerance( num / denom ) * denom;
52 static QgsRectangle transform_to_extent(
const double *geotransform,
double xSize,
double ySize )
56 geotransform[0] + geotransform[1] * xSize,
57 geotransform[3] + geotransform[5] * ySize );
63 static int CPL_STDCALL _progress(
double dfComplete,
const char *pszMessage,
void *pProgressArg )
65 Q_UNUSED( pszMessage );
69 return handler->
progress( dfComplete );
75 static CPLErr rescalePreWarpChunkProcessor(
void *pKern,
void *pArg )
77 GDALWarpKernel *kern = ( GDALWarpKernel * ) pKern;
78 double cellsize = ( (
double * )pArg )[0];
80 for (
int nBand = 0; nBand < kern->nBands; ++nBand )
82 float *bandData = (
float * ) kern->papabySrcImage[nBand];
83 for (
int nLine = 0; nLine < kern->nSrcYSize; ++nLine )
85 for (
int nPixel = 0; nPixel < kern->nSrcXSize; ++nPixel )
87 bandData[nPixel] /= cellsize;
89 bandData += kern->nSrcXSize;
96 static CPLErr rescalePostWarpChunkProcessor(
void *pKern,
void *pArg )
98 GDALWarpKernel *kern = ( GDALWarpKernel * ) pKern;
99 double cellsize = ( (
double * )pArg )[1];
101 for (
int nBand = 0; nBand < kern->nBands; ++nBand )
103 float *bandData = (
float * ) kern->papabyDstImage[nBand];
104 for (
int nLine = 0; nLine < kern->nDstYSize; ++nLine )
106 for (
int nPixel = 0; nPixel < kern->nDstXSize; ++nPixel )
108 bandData[nPixel] *= cellsize;
110 bandData += kern->nDstXSize;
127 for (
int i = 0; i < 6; ++i )
159 if ( customCRSWkt.isEmpty() || customCRSWkt == rasterInfo.
crs() )
164 if ( !customCellSize.isValid() )
175 if ( customGridOffset.x() < 0 || customGridOffset.y() < 0 )
177 if ( !customCellSize.isValid() )
187 mGridOffsetX = fmod_with_tolerance( rasterInfo.
origin().x(), customCellSize.width() );
188 mGridOffsetY = fmod_with_tolerance( rasterInfo.
origin().y(), customCellSize.height() );
203 mCrsWkt = QStringLiteral(
"_error_" );
211 if ( !customCellSize.isValid() )
222 if ( customGridOffset.x() < 0 || customGridOffset.y() < 0 )
241 if (
mCrsWkt == QLatin1String(
"_error_" ) )
249 mErrorMessage = QObject::tr(
"Cell size must not be zero." );
256 double finalExtent[4] = { 0, 0, 0, 0 };
259 for (
int i = 0; i <
mRasters.count(); ++i )
269 mErrorMessage = QString(
"Failed to get suggested warp output.\n\n" 271 "Source WKT:\n%2\n\nDestination WKT:\n%3" )
280 if ( finalExtent[0] == 0 && finalExtent[1] == 0 && finalExtent[2] == 0 && finalExtent[3] == 0 )
291 if ( extent.
xMinimum() > finalExtent[0] ) finalExtent[0] = extent.
xMinimum();
292 if ( extent.
yMinimum() > finalExtent[1] ) finalExtent[1] = extent.
yMinimum();
293 if ( extent.
xMaximum() < finalExtent[2] ) finalExtent[2] = extent.
xMaximum();
294 if ( extent.
yMaximum() < finalExtent[3] ) finalExtent[3] = extent.
yMaximum();
309 if ( clipX0 > finalExtent[0] ) finalExtent[0] = clipX0;
310 if ( clipY0 > finalExtent[1] ) finalExtent[1] = clipY0;
311 if ( clipX1 < finalExtent[2] ) finalExtent[2] = clipX1;
312 if ( clipY1 < finalExtent[3] ) finalExtent[3] = clipY1;
321 int xSize = floor_with_tolerance( ( finalExtent[2] - originX ) /
mCellSizeX );
322 int ySize = floor_with_tolerance( ( finalExtent[3] - originY ) /
mCellSizeY );
324 if ( xSize <= 0 || ySize <= 0 )
326 mErrorMessage = QObject::tr(
"No common intersecting area." );
377 qDebug(
"---ALIGN------------------" );
378 qDebug(
"wkt %s",
mCrsWkt.toLatin1().constData() );
380 qDebug(
"transform" );
385 qDebug(
"extent %s", e.
toString().toLatin1().constData() );
391 double bestCellArea = INFINITY;
399 QString destWkt = destCRS.
toWkt();
406 double cellArea = cs.width() * cs.height();
407 if ( cellArea < bestCellArea )
409 bestCellArea = cellArea;
421 GDALDriverH hDriver = GDALGetDriverByName(
"GTiff" );
424 mErrorMessage = QStringLiteral(
"GDALGetDriverByName(GTiff) failed." );
438 int bandCount = GDALGetRasterCount( hSrcDS.get() );
439 GDALDataType eDT = GDALGetRasterDataType( GDALGetRasterBand( hSrcDS.get(), 1 ) );
443 bandCount, eDT, nullptr ) );
451 GDALSetProjection( hDstDS.get(),
mCrsWkt.toLatin1().constData() );
455 GDALColorTableH hCT = GDALGetRasterColorTable( GDALGetRasterBand( hSrcDS.get(), 1 ) );
457 GDALSetRasterColorTable( GDALGetRasterBand( hDstDS.get(), 1 ), hCT );
463 psWarpOptions->hSrcDS = hSrcDS.get();
464 psWarpOptions->hDstDS = hDstDS.get();
466 psWarpOptions->nBandCount = GDALGetRasterCount( hSrcDS.get() );
467 psWarpOptions->panSrcBands = (
int * ) CPLMalloc(
sizeof(
int ) * psWarpOptions->nBandCount );
468 psWarpOptions->panDstBands = (
int * ) CPLMalloc(
sizeof(
int ) * psWarpOptions->nBandCount );
469 for (
int i = 0; i < psWarpOptions->nBandCount; ++i )
471 psWarpOptions->panSrcBands[i] = i + 1;
472 psWarpOptions->panDstBands[i] = i + 1;
475 psWarpOptions->eResampleAlg =
static_cast< GDALResampleAlg
>( raster.
resampleMethod );
478 psWarpOptions->pfnProgress = _progress;
479 psWarpOptions->pProgressArg =
this;
482 psWarpOptions->pTransformerArg =
483 GDALCreateGenImgProjTransformer( hSrcDS.get(), GDALGetProjectionRef( hSrcDS.get() ),
484 hDstDS.get(), GDALGetProjectionRef( hDstDS.get() ),
486 psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
488 double rescaleArg[2];
493 psWarpOptions->pfnPreWarpChunkProcessor = rescalePreWarpChunkProcessor;
494 psWarpOptions->pfnPostWarpChunkProcessor = rescalePostWarpChunkProcessor;
495 psWarpOptions->pPreWarpProcessorArg = rescaleArg;
496 psWarpOptions->pPostWarpProcessorArg = rescaleArg;
498 psWarpOptions->eWorkingDataType = GDT_Float32;
502 GDALWarpOperation oOperation;
503 oOperation.Initialize( psWarpOptions.get() );
504 oOperation.ChunkAndWarpImage( 0, 0,
mXSize, mYSize );
506 GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
516 void *hTransformArg = GDALCreateGenImgProjTransformer( info.
mDataset.get(), info.
mCrsWkt.toLatin1().constData(),
nullptr, destWkt.toLatin1().constData(), FALSE, 0, 1 );
517 if ( !hTransformArg )
521 double adfDstGeoTransform[6];
523 int nPixels = 0, nLines = 0;
525 eErr = GDALSuggestedWarpOutput2( info.
mDataset.get(),
526 GDALGenImgProjTransform, hTransformArg,
527 adfDstGeoTransform, &nPixels, &nLines, extents, 0 );
528 GDALDestroyGenImgProjTransformer( hTransformArg );
530 if ( eErr != CE_None )
533 QSizeF cs( std::fabs( adfDstGeoTransform[1] ), std::fabs( adfDstGeoTransform[5] ) );
536 *rect =
QgsRectangle( extents[0], extents[1], extents[2], extents[3] );
540 *gridOffset = QPointF( fmod_with_tolerance( adfDstGeoTransform[0], cs.width() ),
541 fmod_with_tolerance( adfDstGeoTransform[3], cs.height() ) );
551 mDataset.reset( GDALOpen( layerpath.toLocal8Bit().constData(), GA_ReadOnly ) );
555 mXSize = GDALGetRasterXSize( mDataset.get() );
556 mYSize = GDALGetRasterYSize( mDataset.get() );
558 ( void ) GDALGetGeoTransform( mDataset.get(),
mGeoTransform );
561 mCrsWkt = QString::fromLatin1( GDALGetProjectionRef( mDataset.get() ) );
563 mBandCnt = GDALGetBandNumber( mDataset.get() );
589 qDebug(
"---RASTER INFO------------------" );
590 qDebug(
"wkt %s",
mCrsWkt.toLatin1().constData() );
595 qDebug(
"extent %s", r.
toString().toLatin1().constData() );
597 qDebug(
"transform" );
604 GDALRasterBandH hBand = GDALGetRasterBand( mDataset.get(), 1 );
610 float *pafScanline = (
float * ) CPLMalloc(
sizeof(
float ) );
611 CPLErr err = GDALRasterIO( hBand, GF_Read, px, py, 1, 1,
612 pafScanline, 1, 1, GDT_Float32, 0, 0 );
613 double value = err == CE_None ? pafScanline[0] : std::numeric_limits<double>::quiet_NaN();
614 CPLFree( pafScanline );
A rectangle specified with double values.
double mCellSizeX
Destination cell size.
QgsRectangle extent() const
Return extent of the raster.
double mGridOffsetX
Destination grid offset - expected to be in interval <0,cellsize)
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.
gdal::dataset_unique_ptr mDataset
handle to open GDAL dataset
bool rescaleValues
rescaling of values according to the change of pixel size
double identify(double mx, double my)
Get raster value at the given coordinates (from the first band)
List mRasters
List of rasters to be aligned (with their output files and other options)
QSizeF cellSize() const
Return cell size in map units.
bool run()
Run the alignment process.
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.
virtual bool progress(double complete)=0
Method to be overridden for progress reporting.
std::unique_ptr< GDALWarpOptions, GDALWarpOptionsDeleter > warp_options_unique_ptr
Scoped GDAL warp options.
QgsRectangle alignedRasterExtent() const
Return expected extent of the resulting aligned raster.
QPointF origin() const
Return origin of the raster.
void setClipExtent(double xmin, double ymin, double xmax, double ymax)
Configure clipping extent (region of interest).
QString outputFilename
filename of the newly created aligned raster (will be overwritten if exists already) ...
int mYSize
Computed raster grid height.
QString crs() const
Return CRS in WKT format.
QSize alignedRasterSize() const
Return expected size of the resulting aligned raster.
Definition of one raster layer for alignment.
bool checkInputParameters()
Determine destination extent from the input rasters and calculate derived values. ...
QString toString(int precision=16) const
Returns a string representation of form xmin,ymin : xmax,ymax Coordinates will be truncated to the sp...
QgsAlignRaster takes one or more raster layers and warps (resamples) them so they have the same: ...
Helper struct to be sub-classed for progress reporting.
double srcCellSizeInDestCRS
used for rescaling of values (if necessary)
QPointF gridOffset() const
Return grid offset.
QString mCrsWkt
CRS stored in WKT format.
void dump() const
write contents of the object to standard error stream - for debugging
std::unique_ptr< void, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
int mXSize
Computed raster grid width.
void dump() const
write contents of the object to standard error stream - for debugging
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 mClipExtent[4]
Optional clip extent: sets "requested area" which be extended to fit the raster grid.
RasterInfo(const QString &layerpath)
Construct raster info with a path to a raster file.
This class represents a coordinate reference system (CRS).
QString toWkt() const
Returns a WKT representation of this CRS.
double xMinimum() const
Returns the x minimum value (left side of rectangle).
double yMaximum() const
Returns the y maximum value (top side of rectangle).
int suggestedReferenceLayer() const
Return index of the layer which has smallest cell size (returns -1 on error)
QString mErrorMessage
Last error message from run()
QgsAlignRaster::ResampleAlg resampleMethod
resampling method to be used
QString inputFilename
filename of the source raster
QgsRectangle clipExtent() const
Get clipping extent (region of interest).
double mGeoTransform[6]
Computed geo-transform.
Utility class for gathering information about rasters.
QString mCrsWkt
Destination CRS - stored in well-known text (WKT) format.
QPointF gridOffset() const
bool createAndWarp(const Item &raster)
Internal function for processing of one raster (1. create output, 2. do the alignment) ...
QSizeF cellSize() const
Get output cell size.