18 #include <gdalwarper.h>
19 #include <ogr_spatialref.h>
33 if ( qAbs( value - qRound( value ) ) < 1e-6 )
34 return qRound( value );
36 return qCeil( value );
41 if ( qAbs( value - qRound( value ) ) < 1e-6 )
42 return qRound( value );
44 return qFloor( value );
57 geotransform[0] + geotransform[1] * xSize,
58 geotransform[3] + geotransform[5] * ySize );
64 static int CPL_STDCALL
_progress(
double dfComplete,
const char* pszMessage,
void* pProgressArg )
66 Q_UNUSED( pszMessage );
70 return handler->
progress( dfComplete );
78 GDALWarpKernel* kern = ( GDALWarpKernel* ) pKern;
79 double cellsize = ((
double* )pArg )[0];
81 for (
int nBand = 0; nBand < kern->nBands; ++nBand )
83 float* bandData = (
float * ) kern->papabySrcImage[nBand];
84 for (
int nLine = 0; nLine < kern->nSrcYSize; ++nLine )
86 for (
int nPixel = 0; nPixel < kern->nSrcXSize; ++nPixel )
88 bandData[nPixel] /= cellsize;
90 bandData += kern->nSrcXSize;
99 GDALWarpKernel* kern = ( GDALWarpKernel* ) pKern;
100 double cellsize = ((
double* )pArg )[1];
102 for (
int nBand = 0; nBand < kern->nBands; ++nBand )
104 float* bandData = (
float * ) kern->papabyDstImage[nBand];
105 for (
int nLine = 0; nLine < kern->nDstYSize; ++nLine )
107 for (
int nPixel = 0; nPixel < kern->nDstXSize; ++nPixel )
109 bandData[nPixel] *= cellsize;
111 bandData += kern->nDstXSize;
120 : mProgressHandler( 0 )
129 for (
int i = 0; i < 6; ++i )
161 if ( customCRSWkt.
isEmpty() || customCRSWkt == rasterInfo.
crs() )
166 if ( !customCellSize.
isValid() )
177 if ( customGridOffset.
x() < 0 || customGridOffset.
y() < 0 )
179 if ( !customCellSize.
isValid() )
213 if ( !customCellSize.
isValid() )
224 if ( customGridOffset.
x() < 0 || customGridOffset.
y() < 0 )
256 for (
int i = 0; i < 6; ++i )
259 double finalExtent[4] = { 0, 0, 0, 0 };
274 "Source WKT:\n%2\n\nDestination WKT:\n%3" )
283 if ( finalExtent[0] == 0 && finalExtent[1] == 0 && finalExtent[2] == 0 && finalExtent[3] == 0 )
294 if ( extent.
xMinimum() > finalExtent[0] ) finalExtent[0] = extent.
xMinimum();
295 if ( extent.
yMinimum() > finalExtent[1] ) finalExtent[1] = extent.
yMinimum();
296 if ( extent.
xMaximum() < finalExtent[2] ) finalExtent[2] = extent.
xMaximum();
297 if ( extent.
yMaximum() < finalExtent[3] ) finalExtent[3] = extent.
yMaximum();
312 if ( clipX0 > finalExtent[0] ) finalExtent[0] = clipX0;
313 if ( clipY0 > finalExtent[1] ) finalExtent[1] = clipY0;
314 if ( clipX1 < finalExtent[2] ) finalExtent[2] = clipX1;
315 if ( clipY1 < finalExtent[3] ) finalExtent[3] = clipY1;
327 if ( xSize <= 0 || ySize <= 0 )
380 qDebug(
"---ALIGN------------------" );
383 qDebug(
"transform" );
394 double bestCellArea = qInf();
410 if ( cellArea < bestCellArea )
412 bestCellArea = cellArea;
424 GDALDriverH hDriver = GDALGetDriverByName(
"GTiff" );
441 int bandCount = GDALGetRasterCount( hSrcDS );
442 GDALDataType eDT = GDALGetRasterDataType( GDALGetRasterBand( hSrcDS, 1 ) );
447 bandCount, eDT, NULL );
460 GDALColorTableH hCT = GDALGetRasterColorTable( GDALGetRasterBand( hSrcDS, 1 ) );
462 GDALSetRasterColorTable( GDALGetRasterBand( hDstDS, 1 ), hCT );
467 GDALWarpOptions* psWarpOptions = GDALCreateWarpOptions();
468 psWarpOptions->hSrcDS = hSrcDS;
469 psWarpOptions->hDstDS = hDstDS;
471 psWarpOptions->nBandCount = GDALGetRasterCount( hSrcDS );
472 psWarpOptions->panSrcBands = (
int * ) CPLMalloc(
sizeof(
int ) * psWarpOptions->nBandCount );
473 psWarpOptions->panDstBands = (
int * ) CPLMalloc(
sizeof(
int ) * psWarpOptions->nBandCount );
474 for (
int i = 0; i < psWarpOptions->nBandCount; ++i )
476 psWarpOptions->panSrcBands[i] = i + 1;
477 psWarpOptions->panDstBands[i] = i + 1;
480 psWarpOptions->eResampleAlg = ( GDALResampleAlg ) raster.
resampleMethod;
484 psWarpOptions->pProgressArg =
this;
487 psWarpOptions->pTransformerArg =
488 GDALCreateGenImgProjTransformer( hSrcDS, GDALGetProjectionRef( hSrcDS ),
489 hDstDS, GDALGetProjectionRef( hDstDS ),
491 psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
493 double rescaleArg[2];
500 psWarpOptions->pPreWarpProcessorArg = rescaleArg;
501 psWarpOptions->pPostWarpProcessorArg = rescaleArg;
503 psWarpOptions->eWorkingDataType = GDT_Float32;
507 GDALWarpOperation oOperation;
508 oOperation.Initialize( psWarpOptions );
509 oOperation.ChunkAndWarpImage( 0, 0,
mXSize, mYSize );
511 GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
512 GDALDestroyWarpOptions( psWarpOptions );
526 if ( !hTransformArg )
530 double adfDstGeoTransform[6];
532 int nPixels = 0, nLines = 0;
534 eErr = GDALSuggestedWarpOutput2( info.
mDataset,
535 GDALGenImgProjTransform, hTransformArg,
536 adfDstGeoTransform, &nPixels, &nLines, extents, 0 );
537 GDALDestroyGenImgProjTransformer( hTransformArg );
539 if ( eErr != CE_None )
542 QSizeF cs( qAbs( adfDstGeoTransform[1] ), qAbs( adfDstGeoTransform[5] ) );
545 *rect =
QgsRectangle( extents[0], extents[1], extents[2], extents[3] );
581 GDALClose( mDataset );
607 qDebug(
"---RASTER INFO------------------" );
615 qDebug(
"transform" );
622 GDALRasterBandH hBand = GDALGetRasterBand( mDataset, 1 );
628 float* pafScanline = (
float * ) CPLMalloc(
sizeof(
float ) );
629 CPLErr err = GDALRasterIO( hBand, GF_Read, px, py, 1, 1,
630 pafScanline, 1, 1, GDT_Float32, 0, 0 );
631 double value = err == CE_None ? pafScanline[0] : std::numeric_limits<double>::quiet_NaN();
632 CPLFree( pafScanline );
QString fromAscii(const char *str, int size)
A rectangle specified with double values.
double mCellSizeX
Destination cell size.
int suggestedReferenceLayer() const
Return index of the layer which has smallest cell size (returns -1 on error)
double mGridOffsetX
Destination grid offset - expected to be in interval <0,cellsize)
QSizeF cellSize() const
Get output cell size.
double yMaximum() const
Get the y maximum value (top side of rectangle)
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.
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)
QString toWkt() const
A helper to get an wkt representation of this srs.
List mRasters
List of rasters to be aligned (with their output files and other options)
QgsRectangle alignedRasterExtent() const
Return expected extent of the resulting aligned raster.
bool run()
Run the alignment process.
static double ceil_with_tolerance(double value)
static double fmod_with_tolerance(double num, double denom)
virtual bool progress(double complete)=0
Method to be overridden for progress reporting.
int mXSize
raster grid size
static CPLErr rescalePreWarpChunkProcessor(void *pKern, void *pArg)
QString tr(const char *sourceText, const char *disambiguation, int n)
QSizeF cellSize() const
Return cell size in map units.
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) ...
static CPLErr rescalePostWarpChunkProcessor(void *pKern, void *pArg)
int mBandCnt
number of raster's bands
void dump() const
write contents of the object to standard error stream - for debugging
int count(const T &value) const
QgsRectangle clipExtent() const
Get clipping extent (region of interest).
double yMinimum() const
Get the y minimum value (bottom side of rectangle)
double xMaximum() const
Get the x maximum value (right side of rectangle)
QPointF gridOffset() const
Return grid offset.
Definition of one raster layer for alignment.
bool checkInputParameters()
Determine destination extent from the input rasters and calculate derived values. ...
const char * constData() const
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)
static double floor_with_tolerance(double value)
QString mCrsWkt
CRS stored in WKT format.
QgsRectangle extent() const
Return extent of the raster.
double mGeoTransform[6]
geotransform coefficients
QByteArray toLocal8Bit() const
int mXSize
Computed raster grid width/height.
static bool suggestedWarpOutput(const RasterInfo &info, const QString &destWkt, QSizeF *cellSize=0, QPointF *gridOffset=0, QgsRectangle *rect=0)
Determine suggested output of raster warp to a different CRS. Returns true on success.
void dump() const
write contents of the object to standard error stream - for debugging
QPointF origin() const
Return origin of the raster.
static QgsRectangle transform_to_extent(const double *geotransform, double xSize, double ySize)
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.
Class for storing a coordinate reference system (CRS)
static int CPL_STDCALL _progress(double dfComplete, const char *pszMessage, void *pProgressArg)
ResampleAlg resampleMethod
resampling method to be used
QSize alignedRasterSize() const
Return expected size of the resulting aligned raster.
QString mErrorMessage
Last error message from run()
QString inputFilename
filename of the source raster
void normalize()
Normalize the rectangle so it has non-negative width/height.
QString crs() const
Return CRS in WKT format.
double mGeoTransform[6]
Computed geo-transform.
QString toString(bool automaticPrecision=false) const
returns string representation of form xmin,ymin xmax,ymax
GDALDatasetH mDataset
handle to open GDAL dataset
QString arg(qlonglong a, int fieldWidth, int base, const QChar &fillChar) const
Utility class for gathering information about rasters.
double xMinimum() const
Get the x minimum value (left side of rectangle)
QString mCrsWkt
Destination CRS - stored in well-known text (WKT) format.
QByteArray toAscii() const
bool createAndWarp(const Item &raster)
Internal function for processing of one raster (1. create output, 2. do the alignment) ...