QGIS API Documentation  3.27.0-Master (bef583a8ef)
qgsalignraster.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsalignraster.cpp
3  --------------------------------------
4  Date : June 2015
5  Copyright : (C) 2015 by Martin Dobias
6  Email : wonder dot sk at gmail dot com
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 
16 #include "qgsalignraster.h"
17 
18 #include <gdalwarper.h>
19 #include <ogr_srs_api.h>
20 #include <cpl_conv.h>
21 #include <limits>
22 
23 #include <QPair>
24 #include <QString>
25 
27 #include "qgsrectangle.h"
28 
29 
30 static double ceil_with_tolerance( double value )
31 {
32  if ( std::fabs( value - std::round( value ) ) < 1e-6 )
33  return std::round( value );
34  else
35  return std::ceil( value );
36 }
37 
38 static double floor_with_tolerance( double value )
39 {
40  if ( std::fabs( value - std::round( value ) ) < 1e-6 )
41  return std::round( value );
42  else
43  return std::floor( value );
44 }
45 
46 static double fmod_with_tolerance( double num, double denom )
47 {
48  return num - floor_with_tolerance( num / denom ) * denom;
49 }
50 
51 
52 static QgsRectangle transform_to_extent( const double *geotransform, double xSize, double ySize )
53 {
54  QgsRectangle r( geotransform[0],
55  geotransform[3],
56  geotransform[0] + geotransform[1] * xSize,
57  geotransform[3] + geotransform[5] * ySize );
58  r.normalize();
59  return r;
60 }
61 
62 
63 static int CPL_STDCALL _progress( double dfComplete, const char *pszMessage, void *pProgressArg )
64 {
65  Q_UNUSED( pszMessage )
66 
67  QgsAlignRaster::ProgressHandler *handler = ( ( QgsAlignRaster * ) pProgressArg )->progressHandler();
68  if ( handler )
69  return handler->progress( dfComplete );
70  else
71  return true;
72 }
73 
74 
75 static CPLErr rescalePreWarpChunkProcessor( void *pKern, void *pArg )
76 {
77  GDALWarpKernel *kern = ( GDALWarpKernel * ) pKern;
78  const double cellsize = ( ( double * )pArg )[0];
79 
80  for ( int nBand = 0; nBand < kern->nBands; ++nBand )
81  {
82  float *bandData = ( float * ) kern->papabySrcImage[nBand];
83  for ( int nLine = 0; nLine < kern->nSrcYSize; ++nLine )
84  {
85  for ( int nPixel = 0; nPixel < kern->nSrcXSize; ++nPixel )
86  {
87  bandData[nPixel] /= cellsize;
88  }
89  bandData += kern->nSrcXSize;
90  }
91  }
92  return CE_None;
93 }
94 
95 
96 static CPLErr rescalePostWarpChunkProcessor( void *pKern, void *pArg )
97 {
98  GDALWarpKernel *kern = ( GDALWarpKernel * ) pKern;
99  const double cellsize = ( ( double * )pArg )[1];
100 
101  for ( int nBand = 0; nBand < kern->nBands; ++nBand )
102  {
103  float *bandData = ( float * ) kern->papabyDstImage[nBand];
104  for ( int nLine = 0; nLine < kern->nDstYSize; ++nLine )
105  {
106  for ( int nPixel = 0; nPixel < kern->nDstXSize; ++nPixel )
107  {
108  bandData[nPixel] *= cellsize;
109  }
110  bandData += kern->nDstXSize;
111  }
112  }
113  return CE_None;
114 }
115 
116 
117 
119 {
120  // parameters
121  mCellSizeX = mCellSizeY = 0;
123  mClipExtent[0] = mClipExtent[1] = mClipExtent[2] = mClipExtent[3] = 0;
124 
125  // derived variables
126  mXSize = mYSize = 0;
127  for ( int i = 0; i < 6; ++i )
128  mGeoTransform[i] = 0;
129 }
130 
131 void QgsAlignRaster::setClipExtent( double xmin, double ymin, double xmax, double ymax )
132 {
133  mClipExtent[0] = xmin;
134  mClipExtent[1] = ymin;
135  mClipExtent[2] = xmax;
136  mClipExtent[3] = ymax;
137 }
138 
140 {
141  setClipExtent( extent.xMinimum(), extent.yMinimum(),
142  extent.xMaximum(), extent.yMaximum() );
143 }
144 
146 {
147  return QgsRectangle( mClipExtent[0], mClipExtent[1],
148  mClipExtent[2], mClipExtent[3] );
149 }
150 
151 
152 bool QgsAlignRaster::setParametersFromRaster( const QString &filename, const QString &destWkt, QSizeF customCellSize, QPointF customGridOffset )
153 {
154  return setParametersFromRaster( RasterInfo( filename ), destWkt, customCellSize, customGridOffset );
155 }
156 
157 bool QgsAlignRaster::setParametersFromRaster( const RasterInfo &rasterInfo, const QString &customCRSWkt, QSizeF customCellSize, QPointF customGridOffset )
158 {
159  if ( customCRSWkt.isEmpty() || customCRSWkt == rasterInfo.crs() )
160  {
161  // use ref. layer to init input
162  mCrsWkt = rasterInfo.crs();
163 
164  if ( !customCellSize.isValid() )
165  {
166  mCellSizeX = rasterInfo.cellSize().width();
167  mCellSizeY = rasterInfo.cellSize().height();
168  }
169  else
170  {
171  mCellSizeX = customCellSize.width();
172  mCellSizeY = customCellSize.height();
173  }
174 
175  if ( customGridOffset.x() < 0 || customGridOffset.y() < 0 )
176  {
177  if ( !customCellSize.isValid() )
178  {
179  // using original raster's grid offset to be aligned with origin
180  mGridOffsetX = rasterInfo.gridOffset().x();
181  mGridOffsetY = rasterInfo.gridOffset().y();
182  }
183  else
184  {
185  // if using custom cell size: offset so that we are aligned
186  // with the original raster's origin point
187  mGridOffsetX = fmod_with_tolerance( rasterInfo.origin().x(), customCellSize.width() );
188  mGridOffsetY = fmod_with_tolerance( rasterInfo.origin().y(), customCellSize.height() );
189  }
190  }
191  else
192  {
193  mGridOffsetX = customGridOffset.x();
194  mGridOffsetY = customGridOffset.y();
195  }
196  }
197  else
198  {
199  QSizeF cs;
200  QPointF go;
201  if ( !suggestedWarpOutput( rasterInfo, customCRSWkt, &cs, &go ) )
202  {
203  mCrsWkt = QStringLiteral( "_error_" );
204  mCellSizeX = mCellSizeY = 0;
206  return false;
207  }
208 
209  mCrsWkt = customCRSWkt;
210 
211  if ( !customCellSize.isValid() )
212  {
213  mCellSizeX = cs.width();
214  mCellSizeY = cs.height();
215  }
216  else
217  {
218  mCellSizeX = customCellSize.width();
219  mCellSizeY = customCellSize.height();
220  }
221 
222  if ( customGridOffset.x() < 0 || customGridOffset.y() < 0 )
223  {
224  mGridOffsetX = go.x();
225  mGridOffsetY = go.y();
226  }
227  else
228  {
229  mGridOffsetX = customGridOffset.x();
230  mGridOffsetY = customGridOffset.x();
231  }
232  }
233  return true;
234 }
235 
236 
238 {
239  mErrorMessage.clear();
240 
241  if ( mCrsWkt == QLatin1String( "_error_" ) )
242  {
243  mErrorMessage = QObject::tr( "Unable to reproject." );
244  return false;
245  }
246 
247  if ( mCellSizeX == 0 || mCellSizeY == 0 )
248  {
249  mErrorMessage = QObject::tr( "Cell size must not be zero." );
250  return false;
251  }
252 
253  mXSize = mYSize = 0;
254  std::fill_n( mGeoTransform, 6, 0 );
255 
256  double finalExtent[4] = { 0, 0, 0, 0 };
257 
258  // for each raster: determine their extent in projected cfg
259  for ( int i = 0; i < mRasters.count(); ++i )
260  {
261  Item &r = mRasters[i];
262 
263  const RasterInfo info( r.inputFilename );
264 
265  QSizeF cs;
266  QgsRectangle extent;
267  if ( !suggestedWarpOutput( info, mCrsWkt, &cs, nullptr, &extent ) )
268  {
269  mErrorMessage = QString( "Failed to get suggested warp output.\n\n"
270  "File:\n%1\n\n"
271  "Source WKT:\n%2\n\nDestination WKT:\n%3" )
272  .arg( r.inputFilename,
273  info.mCrsWkt,
274  mCrsWkt );
275  return false;
276  }
277 
278  r.srcCellSizeInDestCRS = cs.width() * cs.height();
279 
280  if ( finalExtent[0] == 0 && finalExtent[1] == 0 && finalExtent[2] == 0 && finalExtent[3] == 0 )
281  {
282  // initialize with the first layer
283  finalExtent[0] = extent.xMinimum();
284  finalExtent[1] = extent.yMinimum();
285  finalExtent[2] = extent.xMaximum();
286  finalExtent[3] = extent.yMaximum();
287  }
288  else
289  {
290  // use intersection of rects
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();
295  }
296  }
297 
298  // count in extra clip extent (if present)
299  // 1. align requested rect to the grid - extend the rect if necessary
300  // 2. intersect with clip extent with final extent
301 
302  if ( !( mClipExtent[0] == 0 && mClipExtent[1] == 0 && mClipExtent[2] == 0 && mClipExtent[3] == 0 ) )
303  {
304  // extend clip extent to grid
305  const double clipX0 = floor_with_tolerance( ( mClipExtent[0] - mGridOffsetX ) / mCellSizeX ) * mCellSizeX + mGridOffsetX;
306  const double clipY0 = floor_with_tolerance( ( mClipExtent[1] - mGridOffsetY ) / mCellSizeY ) * mCellSizeY + mGridOffsetY;
307  const double clipX1 = ceil_with_tolerance( ( mClipExtent[2] - clipX0 ) / mCellSizeX ) * mCellSizeX + clipX0;
308  const double clipY1 = ceil_with_tolerance( ( mClipExtent[3] - clipY0 ) / mCellSizeY ) * mCellSizeY + clipY0;
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;
313  }
314 
315  // align to grid - shrink the rect if necessary
316  // output raster grid configuration (with no rotation/shear)
317  // ... and raster width/height
318 
319  const double originX = ceil_with_tolerance( ( finalExtent[0] - mGridOffsetX ) / mCellSizeX ) * mCellSizeX + mGridOffsetX;
320  const double originY = ceil_with_tolerance( ( finalExtent[1] - mGridOffsetY ) / mCellSizeY ) * mCellSizeY + mGridOffsetY;
321  const int xSize = floor_with_tolerance( ( finalExtent[2] - originX ) / mCellSizeX );
322  const int ySize = floor_with_tolerance( ( finalExtent[3] - originY ) / mCellSizeY );
323 
324  if ( xSize <= 0 || ySize <= 0 )
325  {
326  mErrorMessage = QObject::tr( "No common intersecting area." );
327  return false;
328  }
329 
330  mXSize = xSize;
331  mYSize = ySize;
332 
333  // build final geotransform...
334  mGeoTransform[0] = originX;
336  mGeoTransform[2] = 0;
337  mGeoTransform[3] = originY + ( mCellSizeY * ySize );
338  mGeoTransform[4] = 0;
340 
341  return true;
342 }
343 
344 
346 {
347  return QSize( mXSize, mYSize );
348 }
349 
351 {
352  return transform_to_extent( mGeoTransform, mXSize, mYSize );
353 }
354 
355 
357 {
358  mErrorMessage.clear();
359 
360  // consider extent of all layers and setup geotransform and output grid size
361  if ( !checkInputParameters() )
362  return false;
363 
364  //dump();
365 
366  const auto constMRasters = mRasters;
367  for ( const Item &r : constMRasters )
368  {
369  if ( !createAndWarp( r ) )
370  return false;
371  }
372  return true;
373 }
374 
375 
377 {
378  qDebug( "---ALIGN------------------" );
379  qDebug( "wkt %s", mCrsWkt.toLatin1().constData() );
380  qDebug( "w/h %d,%d", mXSize, mYSize );
381  qDebug( "transform" );
382  qDebug( "%6.2f %6.2f %6.2f", mGeoTransform[0], mGeoTransform[1], mGeoTransform[2] );
383  qDebug( "%6.2f %6.2f %6.2f", mGeoTransform[3], mGeoTransform[4], mGeoTransform[5] );
384 
385  const QgsRectangle e = transform_to_extent( mGeoTransform, mXSize, mYSize );
386  qDebug( "extent %s", e.toString().toLatin1().constData() );
387 }
388 
390 {
391  int bestIndex = -1;
392  double bestCellArea = INFINITY;
393  QSizeF cs;
394  int i = 0;
395 
396  // using WGS84 as a destination CRS... but maybe some projected CRS
397  // would be a better a choice to more accurately compute areas?
398  // (Why earth is not flat???)
399  const QgsCoordinateReferenceSystem destCRS( QStringLiteral( "EPSG:4326" ) );
400  const QString destWkt = destCRS.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED_GDAL );
401 
402  const auto constMRasters = mRasters;
403  for ( const Item &raster : constMRasters )
404  {
405  if ( !suggestedWarpOutput( RasterInfo( raster.inputFilename ), destWkt, &cs ) )
406  return false;
407 
408  const double cellArea = cs.width() * cs.height();
409  if ( cellArea < bestCellArea )
410  {
411  bestCellArea = cellArea;
412  bestIndex = i;
413  }
414  ++i;
415  }
416 
417  return bestIndex;
418 }
419 
420 
421 bool QgsAlignRaster::createAndWarp( const Item &raster )
422 {
423  GDALDriverH hDriver = GDALGetDriverByName( "GTiff" );
424  if ( !hDriver )
425  {
426  mErrorMessage = QStringLiteral( "GDALGetDriverByName(GTiff) failed." );
427  return false;
428  }
429 
430  // Open the source file.
431  const gdal::dataset_unique_ptr hSrcDS( GDALOpen( raster.inputFilename.toLocal8Bit().constData(), GA_ReadOnly ) );
432  if ( !hSrcDS )
433  {
434  mErrorMessage = QObject::tr( "Unable to open input file: %1" ).arg( raster.inputFilename );
435  return false;
436  }
437 
438  // Create output with same datatype as first input band.
439 
440  const int bandCount = GDALGetRasterCount( hSrcDS.get() );
441  const GDALDataType eDT = GDALGetRasterDataType( GDALGetRasterBand( hSrcDS.get(), 1 ) );
442 
443  // Create the output file.
444  const gdal::dataset_unique_ptr hDstDS( GDALCreate( hDriver, raster.outputFilename.toLocal8Bit().constData(), mXSize, mYSize,
445  bandCount, eDT, nullptr ) );
446  if ( !hDstDS )
447  {
448  mErrorMessage = QObject::tr( "Unable to create output file: %1" ).arg( raster.outputFilename );
449  return false;
450  }
451 
452  // Write out the projection definition.
453  GDALSetProjection( hDstDS.get(), mCrsWkt.toLatin1().constData() );
454  GDALSetGeoTransform( hDstDS.get(), mGeoTransform );
455 
456  // Copy the color table, if required.
457  GDALColorTableH hCT = GDALGetRasterColorTable( GDALGetRasterBand( hSrcDS.get(), 1 ) );
458  if ( hCT )
459  GDALSetRasterColorTable( GDALGetRasterBand( hDstDS.get(), 1 ), hCT );
460 
461  // -----------------------------------------------------------------------
462 
463  // Setup warp options.
464  gdal::warp_options_unique_ptr psWarpOptions( GDALCreateWarpOptions() );
465  psWarpOptions->hSrcDS = hSrcDS.get();
466  psWarpOptions->hDstDS = hDstDS.get();
467 
468  psWarpOptions->nBandCount = GDALGetRasterCount( hSrcDS.get() );
469  psWarpOptions->panSrcBands = ( int * ) CPLMalloc( sizeof( int ) * psWarpOptions->nBandCount );
470  psWarpOptions->panDstBands = ( int * ) CPLMalloc( sizeof( int ) * psWarpOptions->nBandCount );
471  for ( int i = 0; i < psWarpOptions->nBandCount; ++i )
472  {
473  psWarpOptions->panSrcBands[i] = i + 1;
474  psWarpOptions->panDstBands[i] = i + 1;
475  }
476 
477  psWarpOptions->eResampleAlg = static_cast< GDALResampleAlg >( raster.resampleMethod );
478 
479  // our progress function
480  psWarpOptions->pfnProgress = _progress;
481  psWarpOptions->pProgressArg = this;
482 
483  // Establish reprojection transformer.
484  psWarpOptions->pTransformerArg =
485  GDALCreateGenImgProjTransformer( hSrcDS.get(), GDALGetProjectionRef( hSrcDS.get() ),
486  hDstDS.get(), GDALGetProjectionRef( hDstDS.get() ),
487  FALSE, 0.0, 1 );
488  psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
489 
490  double rescaleArg[2];
491  if ( raster.rescaleValues )
492  {
493  rescaleArg[0] = raster.srcCellSizeInDestCRS; // source cell size
494  rescaleArg[1] = mCellSizeX * mCellSizeY; // destination cell size
495  psWarpOptions->pfnPreWarpChunkProcessor = rescalePreWarpChunkProcessor;
496  psWarpOptions->pfnPostWarpChunkProcessor = rescalePostWarpChunkProcessor;
497  psWarpOptions->pPreWarpProcessorArg = rescaleArg;
498  psWarpOptions->pPostWarpProcessorArg = rescaleArg;
499  // force use of float32 data type as that is what our pre/post-processor uses
500  psWarpOptions->eWorkingDataType = GDT_Float32;
501  }
502 
503  // Initialize and execute the warp operation.
504  GDALWarpOperation oOperation;
505  oOperation.Initialize( psWarpOptions.get() );
506  oOperation.ChunkAndWarpImage( 0, 0, mXSize, mYSize );
507 
508  GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
509  return true;
510 }
511 
512 bool QgsAlignRaster::suggestedWarpOutput( const QgsAlignRaster::RasterInfo &info, const QString &destWkt, QSizeF *cellSize, QPointF *gridOffset, QgsRectangle *rect )
513 {
514  // Create a transformer that maps from source pixel/line coordinates
515  // to destination georeferenced coordinates (not destination
516  // pixel line). We do that by omitting the destination dataset
517  // handle (setting it to nullptr).
518  void *hTransformArg = GDALCreateGenImgProjTransformer( info.mDataset.get(), info.mCrsWkt.toLatin1().constData(), nullptr, destWkt.toLatin1().constData(), FALSE, 0, 1 );
519  if ( !hTransformArg )
520  return false;
521 
522  // Get approximate output georeferenced bounds and resolution for file.
523  double adfDstGeoTransform[6];
524  double extents[4];
525  int nPixels = 0, nLines = 0;
526  CPLErr eErr;
527  eErr = GDALSuggestedWarpOutput2( info.mDataset.get(),
528  GDALGenImgProjTransform, hTransformArg,
529  adfDstGeoTransform, &nPixels, &nLines, extents, 0 );
530  GDALDestroyGenImgProjTransformer( hTransformArg );
531 
532  if ( eErr != CE_None )
533  return false;
534 
535  const QSizeF cs( std::fabs( adfDstGeoTransform[1] ), std::fabs( adfDstGeoTransform[5] ) );
536 
537  if ( rect )
538  *rect = QgsRectangle( extents[0], extents[1], extents[2], extents[3] );
539  if ( cellSize )
540  *cellSize = cs;
541  if ( gridOffset )
542  *gridOffset = QPointF( fmod_with_tolerance( adfDstGeoTransform[0], cs.width() ),
543  fmod_with_tolerance( adfDstGeoTransform[3], cs.height() ) );
544  return true;
545 }
546 
547 
548 //----------
549 
550 
551 QgsAlignRaster::RasterInfo::RasterInfo( const QString &layerpath )
552 {
553  mDataset.reset( GDALOpen( layerpath.toLocal8Bit().constData(), GA_ReadOnly ) );
554  if ( !mDataset )
555  return;
556 
557  mXSize = GDALGetRasterXSize( mDataset.get() );
558  mYSize = GDALGetRasterYSize( mDataset.get() );
559 
560  ( void ) GDALGetGeoTransform( mDataset.get(), mGeoTransform );
561 
562  // TODO: may be null or empty string
563  mCrsWkt = QString::fromLatin1( GDALGetProjectionRef( mDataset.get() ) );
564 
565  mBandCnt = GDALGetBandNumber( mDataset.get() );
566 }
567 
569 {
570  return QSizeF( std::fabs( mGeoTransform[1] ), std::fabs( mGeoTransform[5] ) );
571 }
572 
574 {
575  return QPointF( fmod_with_tolerance( mGeoTransform[0], cellSize().width() ),
576  fmod_with_tolerance( mGeoTransform[3], cellSize().height() ) );
577 }
578 
580 {
581  return transform_to_extent( mGeoTransform, mXSize, mYSize );
582 }
583 
585 {
586  return QPointF( mGeoTransform[0], mGeoTransform[3] );
587 }
588 
590 {
591  qDebug( "---RASTER INFO------------------" );
592  qDebug( "wkt %s", mCrsWkt.toLatin1().constData() );
593  qDebug( "w/h %d,%d", mXSize, mYSize );
594  qDebug( "cell x/y %f,%f", cellSize().width(), cellSize().width() );
595 
596  const QgsRectangle r = extent();
597  qDebug( "extent %s", r.toString().toLatin1().constData() );
598 
599  qDebug( "transform" );
600  qDebug( "%6.2f %6.2f %6.2f", mGeoTransform[0], mGeoTransform[1], mGeoTransform[2] );
601  qDebug( "%6.2f %6.2f %6.2f", mGeoTransform[3], mGeoTransform[4], mGeoTransform[5] );
602 }
603 
604 double QgsAlignRaster::RasterInfo::identify( double mx, double my )
605 {
606  GDALRasterBandH hBand = GDALGetRasterBand( mDataset.get(), 1 );
607 
608  // must not be rotated in order for this to work
609  const int px = int( ( mx - mGeoTransform[0] ) / mGeoTransform[1] );
610  const int py = int( ( my - mGeoTransform[3] ) / mGeoTransform[5] );
611 
612  float *pafScanline = ( float * ) CPLMalloc( sizeof( float ) );
613  const CPLErr err = GDALRasterIO( hBand, GF_Read, px, py, 1, 1,
614  pafScanline, 1, 1, GDT_Float32, 0, 0 );
615  const double value = err == CE_None ? pafScanline[0] : std::numeric_limits<double>::quiet_NaN();
616  CPLFree( pafScanline );
617 
618  return value;
619 }
620 
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).
@ WKT_PREFERRED_GDAL
Preferred format for conversion of CRS to WKT for use with the GDAL library.
QString toWkt(WktVariant variant=WKT1_GDAL, bool multiline=false, int indentationWidth=4) const
Returns a WKT representation of this CRS.
A rectangle specified with double values.
Definition: qgsrectangle.h:42
QString toString(int precision=16) const
Returns a string representation of form xmin,ymin : xmax,ymax Coordinates will be truncated to the sp...
double yMaximum() const SIP_HOLDGIL
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:193
double xMaximum() const SIP_HOLDGIL
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:183
double xMinimum() const SIP_HOLDGIL
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:188
double yMinimum() const SIP_HOLDGIL
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:198
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
Definition: qgsogrutils.h:139
std::unique_ptr< GDALWarpOptions, GDALWarpOptionsDeleter > warp_options_unique_ptr
Scoped GDAL warp options.
Definition: qgsogrutils.h:154
Definition of one raster layer for alignment.
bool rescaleValues
rescaling of values according to the change of pixel size
QgsAlignRaster::ResampleAlg resampleMethod
resampling method to be used
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.