QGIS API Documentation  2.17.0-Master (8784312)
qgsrasterprojector.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsrasterprojector.cpp - Raster projector
3  --------------------------------------
4  Date : Jan 16, 2011
5  Copyright : (C) 2005 by Radim Blazek
6  email : radim dot blazek at gmail dot com
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 #include <algorithm>
18 
19 #include "qgsrasterdataprovider.h"
20 #include "qgscrscache.h"
21 #include "qgslogger.h"
22 #include "qgsrasterprojector.h"
23 #include "qgscoordinatetransform.h"
24 
26  const QgsCoordinateReferenceSystem& theSrcCRS,
27  const QgsCoordinateReferenceSystem& theDestCRS,
28  int theSrcDatumTransform,
29  int theDestDatumTransform,
30  const QgsRectangle& theDestExtent,
31  int theDestRows, int theDestCols,
32  double theMaxSrcXRes, double theMaxSrcYRes,
33  const QgsRectangle& theExtent )
34  : QgsRasterInterface( nullptr )
35  , mSrcCRS( theSrcCRS )
36  , mDestCRS( theDestCRS )
37  , mSrcDatumTransform( theSrcDatumTransform )
38  , mDestDatumTransform( theDestDatumTransform )
39  , mDestExtent( theDestExtent )
40  , mExtent( theExtent )
41  , mDestRows( theDestRows ), mDestCols( theDestCols )
42  , pHelperTop( nullptr ), pHelperBottom( nullptr )
43  , mMaxSrcXRes( theMaxSrcXRes ), mMaxSrcYRes( theMaxSrcYRes )
44  , mPrecision( Approximate )
45  , mApproximate( true )
46 {
47  QgsDebugMsgLevel( "Entered", 4 );
48  QgsDebugMsgLevel( "theDestExtent = " + theDestExtent.toString(), 4 );
49 
50  calc();
51 }
52 
54  const QgsCoordinateReferenceSystem& theSrcCRS,
55  const QgsCoordinateReferenceSystem& theDestCRS,
56  const QgsRectangle& theDestExtent,
57  int theDestRows, int theDestCols,
58  double theMaxSrcXRes, double theMaxSrcYRes,
59  const QgsRectangle& theExtent )
60  : QgsRasterInterface( nullptr )
61  , mSrcCRS( theSrcCRS )
62  , mDestCRS( theDestCRS )
63  , mSrcDatumTransform( -1 )
64  , mDestDatumTransform( -1 )
65  , mDestExtent( theDestExtent )
66  , mExtent( theExtent )
67  , mDestRows( theDestRows ), mDestCols( theDestCols )
68  , pHelperTop( nullptr ), pHelperBottom( nullptr )
69  , mMaxSrcXRes( theMaxSrcXRes ), mMaxSrcYRes( theMaxSrcYRes )
70  , mPrecision( Approximate )
71  , mApproximate( false )
72 {
73  QgsDebugMsgLevel( "Entered", 4 );
74  QgsDebugMsgLevel( "theDestExtent = " + theDestExtent.toString(), 4 );
75 
76  calc();
77 }
78 
80  const QgsCoordinateReferenceSystem& theSrcCRS,
81  const QgsCoordinateReferenceSystem& theDestCRS,
82  double theMaxSrcXRes, double theMaxSrcYRes,
83  const QgsRectangle& theExtent )
84  : QgsRasterInterface( nullptr )
85  , mSrcCRS( theSrcCRS )
86  , mDestCRS( theDestCRS )
87  , mSrcDatumTransform( -1 )
88  , mDestDatumTransform( -1 )
89  , mExtent( theExtent )
90  , mDestRows( 0 )
91  , mDestCols( 0 )
92  , mDestXRes( 0.0 )
93  , mDestYRes( 0.0 )
94  , mSrcRows( 0 )
95  , mSrcCols( 0 )
96  , mSrcXRes( 0.0 )
97  , mSrcYRes( 0.0 )
98  , mDestRowsPerMatrixRow( 0.0 )
99  , mDestColsPerMatrixCol( 0.0 )
100  , pHelperTop( nullptr ), pHelperBottom( nullptr )
101  , mHelperTopRow( 0 )
102  , mCPCols( 0 )
103  , mCPRows( 0 )
104  , mSqrTolerance( 0.0 )
105  , mMaxSrcXRes( theMaxSrcXRes )
106  , mMaxSrcYRes( theMaxSrcYRes )
107  , mPrecision( Approximate )
108  , mApproximate( false )
109 {
110  QgsDebugMsgLevel( "Entered", 4 );
111 }
112 
114  : QgsRasterInterface( nullptr )
115  , mSrcDatumTransform( -1 )
116  , mDestDatumTransform( -1 )
117  , mDestRows( 0 )
118  , mDestCols( 0 )
119  , mDestXRes( 0.0 )
120  , mDestYRes( 0.0 )
121  , mSrcRows( 0 )
122  , mSrcCols( 0 )
123  , mSrcXRes( 0.0 )
124  , mSrcYRes( 0.0 )
125  , mDestRowsPerMatrixRow( 0.0 )
126  , mDestColsPerMatrixCol( 0.0 )
127  , pHelperTop( nullptr )
128  , pHelperBottom( nullptr )
129  , mHelperTopRow( 0 )
130  , mCPCols( 0 )
131  , mCPRows( 0 )
132  , mSqrTolerance( 0.0 )
133  , mMaxSrcXRes( 0 )
134  , mMaxSrcYRes( 0 )
135  , mPrecision( Approximate )
136  , mApproximate( false )
137 {
138  QgsDebugMsgLevel( "Entered", 4 );
139 }
140 
142  : QgsRasterInterface( nullptr )
143  , pHelperTop( nullptr )
144  , pHelperBottom( nullptr )
145  , mHelperTopRow( 0 )
146  , mCPCols( 0 )
147  , mCPRows( 0 )
148  , mSqrTolerance( 0 )
149  , mApproximate( false )
150 {
151  mSrcCRS = projector.mSrcCRS;
152  mDestCRS = projector.mDestCRS;
153  mSrcDatumTransform = projector.mSrcDatumTransform;
154  mDestDatumTransform = projector.mDestDatumTransform;
155  mMaxSrcXRes = projector.mMaxSrcXRes;
156  mMaxSrcYRes = projector.mMaxSrcYRes;
157  mExtent = projector.mExtent;
158  mDestRows = projector.mDestRows;
159  mDestCols = projector.mDestCols;
160  mDestXRes = projector.mDestXRes;
161  mDestYRes = projector.mDestYRes;
162  mSrcRows = projector.mSrcRows;
163  mSrcCols = projector.mSrcCols;
164  mSrcXRes = projector.mSrcXRes;
165  mSrcYRes = projector.mSrcYRes;
166  mDestRowsPerMatrixRow = projector.mDestRowsPerMatrixRow;
167  mDestColsPerMatrixCol = projector.mDestColsPerMatrixCol;
168  mPrecision = projector.mPrecision;
169 }
170 
172 {
173  if ( &projector != this )
174  {
175  mSrcCRS = projector.mSrcCRS;
176  mDestCRS = projector.mDestCRS;
177  mSrcDatumTransform = projector.mSrcDatumTransform;
178  mDestDatumTransform = projector.mDestDatumTransform;
179  mMaxSrcXRes = projector.mMaxSrcXRes;
180  mMaxSrcYRes = projector.mMaxSrcYRes;
181  mExtent = projector.mExtent;
182  mPrecision = projector.mPrecision;
183  }
184  return *this;
185 }
186 
188 {
189  QgsDebugMsgLevel( "Entered", 4 );
190  QgsRasterProjector * projector = new QgsRasterProjector( mSrcCRS, mDestCRS, mMaxSrcXRes, mMaxSrcYRes, mExtent );
191  projector->mSrcDatumTransform = mSrcDatumTransform;
192  projector->mDestDatumTransform = mDestDatumTransform;
193  projector->mPrecision = mPrecision;
194  return projector;
195 }
196 
198 {
199  delete[] pHelperTop;
200  delete[] pHelperBottom;
201 }
202 
204 {
205  if ( mInput ) return mInput->bandCount();
206 
207  return 0;
208 }
209 
211 {
212  if ( mInput ) return mInput->dataType( bandNo );
213 
214  return QGis::UnknownDataType;
215 }
216 
217 void QgsRasterProjector::setCRS( const QgsCoordinateReferenceSystem & theSrcCRS, const QgsCoordinateReferenceSystem & theDestCRS, int srcDatumTransform, int destDatumTransform )
218 {
219  mSrcCRS = theSrcCRS;
220  mDestCRS = theDestCRS;
221  mSrcDatumTransform = srcDatumTransform;
222  mDestDatumTransform = destDatumTransform;
223 }
224 
225 void QgsRasterProjector::calc()
226 {
227  QgsDebugMsgLevel( "Entered", 4 );
228  mCPMatrix.clear();
229  mCPLegalMatrix.clear();
230  delete[] pHelperTop;
231  pHelperTop = nullptr;
232  delete[] pHelperBottom;
233  pHelperBottom = nullptr;
234 
235  // Get max source resolution and extent if possible
236  mMaxSrcXRes = 0;
237  mMaxSrcYRes = 0;
238  if ( mInput )
239  {
240  QgsRasterDataProvider *provider = dynamic_cast<QgsRasterDataProvider*>( mInput->srcInput() );
241  if ( provider )
242  {
243  if ( provider->capabilities() & QgsRasterDataProvider::Size )
244  {
245  mMaxSrcXRes = provider->extent().width() / provider->xSize();
246  mMaxSrcYRes = provider->extent().height() / provider->ySize();
247  }
248  // Get source extent
249  if ( mExtent.isEmpty() )
250  {
251  mExtent = provider->extent();
252  }
253  }
254  }
255 
256  mDestXRes = mDestExtent.width() / ( mDestCols );
257  mDestYRes = mDestExtent.height() / ( mDestRows );
258 
259  // Calculate tolerance
260  // TODO: Think it over better
261  // Note: we are checking on matrix each even point, that means that the real error
262  // in that moment is approximately half size
263  double myDestRes = mDestXRes < mDestYRes ? mDestXRes : mDestYRes;
264  mSqrTolerance = myDestRes * myDestRes;
265 
266  const QgsCoordinateTransform* inverseCt = QgsCoordinateTransformCache::instance()->transform( mDestCRS.authid(), mSrcCRS.authid(), mDestDatumTransform, mSrcDatumTransform );
267 
268  if ( mPrecision == Approximate )
269  {
270  mApproximate = true;
271  }
272  else
273  {
274  mApproximate = false;
275  }
276 
277  // Always try to calculate mCPMatrix, it is used in calcSrcExtent() for both Approximate and Exact
278  // Initialize the matrix by corners and middle points
279  mCPCols = mCPRows = 3;
280  for ( int i = 0; i < mCPRows; i++ )
281  {
282  QList<QgsPoint> myRow;
283  myRow.append( QgsPoint() );
284  myRow.append( QgsPoint() );
285  myRow.append( QgsPoint() );
286  mCPMatrix.insert( i, myRow );
287  // And the legal points
288  QList<bool> myLegalRow;
289  myLegalRow.append( bool( false ) );
290  myLegalRow.append( bool( false ) );
291  myLegalRow.append( bool( false ) );
292  mCPLegalMatrix.insert( i, myLegalRow );
293  }
294  for ( int i = 0; i < mCPRows; i++ )
295  {
296  calcRow( i, inverseCt );
297  }
298 
299  while ( true )
300  {
301  bool myColsOK = checkCols( inverseCt );
302  if ( !myColsOK )
303  {
304  insertRows( inverseCt );
305  }
306  bool myRowsOK = checkRows( inverseCt );
307  if ( !myRowsOK )
308  {
309  insertCols( inverseCt );
310  }
311  if ( myColsOK && myRowsOK )
312  {
313  QgsDebugMsgLevel( "CP matrix within tolerance", 4 );
314  break;
315  }
316  // What is the maximum reasonable size of transformatio matrix?
317  // TODO: consider better when to break - ratio
318  if ( mCPRows * mCPCols > 0.25 * mDestRows * mDestCols )
319  //if ( mCPRows * mCPCols > mDestRows * mDestCols )
320  {
321  QgsDebugMsgLevel( "Too large CP matrix", 4 );
322  mApproximate = false;
323  break;
324  }
325  }
326  QgsDebugMsgLevel( QString( "CPMatrix size: mCPRows = %1 mCPCols = %2" ).arg( mCPRows ).arg( mCPCols ), 4 );
327  mDestRowsPerMatrixRow = static_cast< float >( mDestRows ) / ( mCPRows - 1 );
328  mDestColsPerMatrixCol = static_cast< float >( mDestCols ) / ( mCPCols - 1 );
329 
330  QgsDebugMsgLevel( "CPMatrix:", 5 );
331  QgsDebugMsgLevel( cpToString(), 5 );
332 
333  // init helper points
334  pHelperTop = new QgsPoint[mDestCols];
335  pHelperBottom = new QgsPoint[mDestCols];
336  calcHelper( 0, pHelperTop );
337  calcHelper( 1, pHelperBottom );
338  mHelperTopRow = 0;
339 
340  // Calculate source dimensions
341  calcSrcExtent();
342  calcSrcRowsCols();
343  mSrcYRes = mSrcExtent.height() / mSrcRows;
344  mSrcXRes = mSrcExtent.width() / mSrcCols;
345 }
346 
347 void QgsRasterProjector::calcSrcExtent()
348 {
349  /* Run around the mCPMatrix and find source extent */
350  // Attention, source limits are not necessarily on destination edges, e.g.
351  // for destination EPSG:32661 Polar Stereographic and source EPSG:4326,
352  // the maximum y may be in the middle of destination extent
353  // TODO: How to find extent exactly and quickly?
354  // For now, we run through all matrix
355  // mCPMatrix is used for both Approximate and Exact because QgsCoordinateTransform::transformBoundingBox()
356  // is not precise enough, see #13665
357  QgsPoint myPoint = mCPMatrix[0][0];
358  mSrcExtent = QgsRectangle( myPoint.x(), myPoint.y(), myPoint.x(), myPoint.y() );
359  for ( int i = 0; i < mCPRows; i++ )
360  {
361  for ( int j = 0; j < mCPCols ; j++ )
362  {
363  myPoint = mCPMatrix[i][j];
364  if ( mCPLegalMatrix[i][j] )
365  {
366  mSrcExtent.combineExtentWith( myPoint.x(), myPoint.y() );
367  }
368  }
369  }
370  // Expand a bit to avoid possible approx coords falling out because of representation error?
371 
372  // Combine with maximum source extent
373  mSrcExtent = mSrcExtent.intersect( &mExtent );
374 
375  // If mMaxSrcXRes, mMaxSrcYRes are defined (fixed src resolution)
376  // align extent to src resolution to avoid jumping of reprojected pixels
377  // when shifting resampled grid.
378  // Important especially if we are over mMaxSrcXRes, mMaxSrcYRes limits
379  // Note however, that preceding filters (like resampler) may read data
380  // on different resolution.
381 
382  QgsDebugMsgLevel( "mSrcExtent = " + mSrcExtent.toString(), 4 );
383  QgsDebugMsgLevel( "mExtent = " + mExtent.toString(), 4 );
384  if ( !mExtent.isEmpty() )
385  {
386  if ( mMaxSrcXRes > 0 )
387  {
388  // with floor/ceil it should work correctly also for mSrcExtent.xMinimum() < mExtent.xMinimum()
389  double col = floor(( mSrcExtent.xMinimum() - mExtent.xMinimum() ) / mMaxSrcXRes );
390  double x = mExtent.xMinimum() + col * mMaxSrcXRes;
391  mSrcExtent.setXMinimum( x );
392 
393  col = ceil(( mSrcExtent.xMaximum() - mExtent.xMinimum() ) / mMaxSrcXRes );
394  x = mExtent.xMinimum() + col * mMaxSrcXRes;
395  mSrcExtent.setXMaximum( x );
396  }
397  if ( mMaxSrcYRes > 0 )
398  {
399  double row = floor(( mExtent.yMaximum() - mSrcExtent.yMaximum() ) / mMaxSrcYRes );
400  double y = mExtent.yMaximum() - row * mMaxSrcYRes;
401  mSrcExtent.setYMaximum( y );
402 
403  row = ceil(( mExtent.yMaximum() - mSrcExtent.yMinimum() ) / mMaxSrcYRes );
404  y = mExtent.yMaximum() - row * mMaxSrcYRes;
405  mSrcExtent.setYMinimum( y );
406  }
407  }
408  QgsDebugMsgLevel( "mSrcExtent = " + mSrcExtent.toString(), 4 );
409 }
410 
411 QString QgsRasterProjector::cpToString()
412 {
413  QString myString;
414  for ( int i = 0; i < mCPRows; i++ )
415  {
416  if ( i > 0 )
417  myString += '\n';
418  for ( int j = 0; j < mCPCols; j++ )
419  {
420  if ( j > 0 )
421  myString += " ";
422  QgsPoint myPoint = mCPMatrix[i][j];
423  if ( mCPLegalMatrix[i][j] )
424  {
425  myString += myPoint.toString();
426  }
427  else
428  {
429  myString += "(-,-)";
430  }
431  }
432  }
433  return myString;
434 }
435 
436 void QgsRasterProjector::calcSrcRowsCols()
437 {
438  // Wee need to calculate minimum cell size in the source
439  // TODO: Think it over better, what is the right source resolution?
440  // Taking distances between cell centers projected to source along source
441  // axis would result in very high resolution
442  // TODO: different resolution for rows and cols ?
443 
444  double myMinSize = std::numeric_limits<double>::max();
445 
446  if ( mApproximate )
447  {
448  // For now, we take cell sizes projected to source but not to source axes
449  double myDestColsPerMatrixCell = static_cast< double >( mDestCols ) / mCPCols;
450  double myDestRowsPerMatrixCell = static_cast< double >( mDestRows ) / mCPRows;
451  QgsDebugMsgLevel( QString( "myDestColsPerMatrixCell = %1 myDestRowsPerMatrixCell = %2" ).arg( myDestColsPerMatrixCell ).arg( myDestRowsPerMatrixCell ), 4 );
452  for ( int i = 0; i < mCPRows - 1; i++ )
453  {
454  for ( int j = 0; j < mCPCols - 1; j++ )
455  {
456  QgsPoint myPointA = mCPMatrix[i][j];
457  QgsPoint myPointB = mCPMatrix[i][j+1];
458  QgsPoint myPointC = mCPMatrix[i+1][j];
459  if ( mCPLegalMatrix[i][j] && mCPLegalMatrix[i][j+1] && mCPLegalMatrix[i+1][j] )
460  {
461  double mySize = sqrt( myPointA.sqrDist( myPointB ) ) / myDestColsPerMatrixCell;
462  if ( mySize < myMinSize )
463  myMinSize = mySize;
464 
465  mySize = sqrt( myPointA.sqrDist( myPointC ) ) / myDestRowsPerMatrixCell;
466  if ( mySize < myMinSize )
467  myMinSize = mySize;
468  }
469  }
470  }
471  }
472  else
473  {
474  // take highest from corners, points in in the middle of corners and center (3 x 3 )
475  const QgsCoordinateTransform* inverseCt = QgsCoordinateTransformCache::instance()->transform( mDestCRS.authid(), mSrcCRS.authid(), mDestDatumTransform, mSrcDatumTransform );
476  //double
477  QgsRectangle srcExtent;
478  int srcXSize, srcYSize;
479  if ( extentSize( inverseCt, mDestExtent, mDestCols, mDestRows, srcExtent, srcXSize, srcYSize ) )
480  {
481  double srcXRes = srcExtent.width() / srcXSize;
482  double srcYRes = srcExtent.height() / srcYSize;
483  myMinSize = std::min( srcXRes, srcYRes );
484  }
485  else
486  {
487  QgsDebugMsg( "Cannot get src extent/size" );
488  }
489  }
490 
491  // Make it a bit higher resolution
492  // TODO: find the best coefficient, attention, increasing resolution for WMS
493  // is changing WMS content
494  myMinSize *= 0.75;
495 
496  QgsDebugMsgLevel( QString( "mMaxSrcXRes = %1 mMaxSrcYRes = %2" ).arg( mMaxSrcXRes ).arg( mMaxSrcYRes ), 4 );
497  // mMaxSrcXRes, mMaxSrcYRes may be 0 - no limit (WMS)
498  double myMinXSize = mMaxSrcXRes > myMinSize ? mMaxSrcXRes : myMinSize;
499  double myMinYSize = mMaxSrcYRes > myMinSize ? mMaxSrcYRes : myMinSize;
500  QgsDebugMsgLevel( QString( "myMinXSize = %1 myMinYSize = %2" ).arg( myMinXSize ).arg( myMinYSize ), 4 );
501  QgsDebugMsgLevel( QString( "mSrcExtent.width = %1 mSrcExtent.height = %2" ).arg( mSrcExtent.width() ).arg( mSrcExtent.height() ), 4 );
502 
503  // we have to round to keep alignment set in calcSrcExtent
504  mSrcRows = static_cast< int >( qRound( mSrcExtent.height() / myMinYSize ) );
505  mSrcCols = static_cast< int >( qRound( mSrcExtent.width() / myMinXSize ) );
506 
507  QgsDebugMsgLevel( QString( "mSrcRows = %1 mSrcCols = %2" ).arg( mSrcRows ).arg( mSrcCols ), 4 );
508 }
509 
510 
511 inline void QgsRasterProjector::destPointOnCPMatrix( int theRow, int theCol, double *theX, double *theY )
512 {
513  *theX = mDestExtent.xMinimum() + theCol * mDestExtent.width() / ( mCPCols - 1 );
514  *theY = mDestExtent.yMaximum() - theRow * mDestExtent.height() / ( mCPRows - 1 );
515 }
516 
517 inline int QgsRasterProjector::matrixRow( int theDestRow )
518 {
519  return static_cast< int >( floor(( theDestRow + 0.5 ) / mDestRowsPerMatrixRow ) );
520 }
521 inline int QgsRasterProjector::matrixCol( int theDestCol )
522 {
523  return static_cast< int >( floor(( theDestCol + 0.5 ) / mDestColsPerMatrixCol ) );
524 }
525 
526 QgsPoint QgsRasterProjector::srcPoint( int theDestRow, int theCol )
527 {
528  Q_UNUSED( theDestRow );
529  Q_UNUSED( theCol );
530  return QgsPoint();
531 }
532 
533 void QgsRasterProjector::calcHelper( int theMatrixRow, QgsPoint *thePoints )
534 {
535  // TODO?: should we also precalc dest cell center coordinates for x and y?
536  for ( int myDestCol = 0; myDestCol < mDestCols; myDestCol++ )
537  {
538  double myDestX = mDestExtent.xMinimum() + ( myDestCol + 0.5 ) * mDestXRes;
539 
540  int myMatrixCol = matrixCol( myDestCol );
541 
542  double myDestXMin, myDestYMin, myDestXMax, myDestYMax;
543 
544  destPointOnCPMatrix( theMatrixRow, myMatrixCol, &myDestXMin, &myDestYMin );
545  destPointOnCPMatrix( theMatrixRow, myMatrixCol + 1, &myDestXMax, &myDestYMax );
546 
547  double xfrac = ( myDestX - myDestXMin ) / ( myDestXMax - myDestXMin );
548 
549  QgsPoint &mySrcPoint0 = mCPMatrix[theMatrixRow][myMatrixCol];
550  QgsPoint &mySrcPoint1 = mCPMatrix[theMatrixRow][myMatrixCol+1];
551  double s = mySrcPoint0.x() + ( mySrcPoint1.x() - mySrcPoint0.x() ) * xfrac;
552  double t = mySrcPoint0.y() + ( mySrcPoint1.y() - mySrcPoint0.y() ) * xfrac;
553 
554  thePoints[myDestCol].setX( s );
555  thePoints[myDestCol].setY( t );
556  }
557 }
558 void QgsRasterProjector::nextHelper()
559 {
560  // We just switch pHelperTop and pHelperBottom, memory is not lost
561  QgsPoint *tmp;
562  tmp = pHelperTop;
563  pHelperTop = pHelperBottom;
564  pHelperBottom = tmp;
565  calcHelper( mHelperTopRow + 2, pHelperBottom );
566  mHelperTopRow++;
567 }
568 
569 bool QgsRasterProjector::srcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform* ct )
570 {
571  if ( mApproximate )
572  {
573  return approximateSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol );
574  }
575  else
576  {
577  return preciseSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol, ct );
578  }
579 }
580 
581 bool QgsRasterProjector::preciseSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform* ct )
582 {
583 #ifdef QGISDEBUG
584  QgsDebugMsgLevel( QString( "theDestRow = %1" ).arg( theDestRow ), 5 );
585  QgsDebugMsgLevel( QString( "theDestRow = %1 mDestExtent.yMaximum() = %2 mDestYRes = %3" ).arg( theDestRow ).arg( mDestExtent.yMaximum() ).arg( mDestYRes ), 5 );
586 #endif
587 
588  // Get coordinate of center of destination cell
589  double x = mDestExtent.xMinimum() + ( theDestCol + 0.5 ) * mDestXRes;
590  double y = mDestExtent.yMaximum() - ( theDestRow + 0.5 ) * mDestYRes;
591  double z = 0;
592 
593 #ifdef QGISDEBUG
594  QgsDebugMsgLevel( QString( "x = %1 y = %2" ).arg( x ).arg( y ), 5 );
595 #endif
596 
597  if ( ct )
598  {
599  ct->transformInPlace( x, y, z );
600  }
601 
602 #ifdef QGISDEBUG
603  QgsDebugMsgLevel( QString( "x = %1 y = %2" ).arg( x ).arg( y ), 5 );
604 #endif
605 
606  if ( !mExtent.contains( QgsPoint( x, y ) ) )
607  {
608  return false;
609  }
610  // Get source row col
611  *theSrcRow = static_cast< int >( floor(( mSrcExtent.yMaximum() - y ) / mSrcYRes ) );
612  *theSrcCol = static_cast< int >( floor(( x - mSrcExtent.xMinimum() ) / mSrcXRes ) );
613 #ifdef QGISDEBUG
614  QgsDebugMsgLevel( QString( "mSrcExtent.yMinimum() = %1 mSrcExtent.yMaximum() = %2 mSrcYRes = %3" ).arg( mSrcExtent.yMinimum() ).arg( mSrcExtent.yMaximum() ).arg( mSrcYRes ), 5 );
615  QgsDebugMsgLevel( QString( "theSrcRow = %1 theSrcCol = %2" ).arg( *theSrcRow ).arg( *theSrcCol ), 5 );
616 #endif
617 
618  // With epsg 32661 (Polar Stereographic) it was happening that *theSrcCol == mSrcCols
619  // For now silently correct limits to avoid crashes
620  // TODO: review
621  // should not happen
622  if ( *theSrcRow >= mSrcRows ) return false;
623  if ( *theSrcRow < 0 ) return false;
624  if ( *theSrcCol >= mSrcCols ) return false;
625  if ( *theSrcCol < 0 ) return false;
626 
627  return true;
628 }
629 
630 bool QgsRasterProjector::approximateSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol )
631 {
632  int myMatrixRow = matrixRow( theDestRow );
633  int myMatrixCol = matrixCol( theDestCol );
634 
635  if ( myMatrixRow > mHelperTopRow )
636  {
637  // TODO: make it more robust (for random, not sequential reading)
638  nextHelper();
639  }
640 
641  double myDestY = mDestExtent.yMaximum() - ( theDestRow + 0.5 ) * mDestYRes;
642 
643  // See the schema in javax.media.jai.WarpGrid doc (but up side down)
644  // TODO: use some kind of cache of values which can be reused
645  double myDestXMin, myDestYMin, myDestXMax, myDestYMax;
646 
647  destPointOnCPMatrix( myMatrixRow + 1, myMatrixCol, &myDestXMin, &myDestYMin );
648  destPointOnCPMatrix( myMatrixRow, myMatrixCol + 1, &myDestXMax, &myDestYMax );
649 
650  double yfrac = ( myDestY - myDestYMin ) / ( myDestYMax - myDestYMin );
651 
652  QgsPoint &myTop = pHelperTop[theDestCol];
653  QgsPoint &myBot = pHelperBottom[theDestCol];
654 
655  // Warning: this is very SLOW compared to the following code!:
656  //double mySrcX = myBot.x() + (myTop.x() - myBot.x()) * yfrac;
657  //double mySrcY = myBot.y() + (myTop.y() - myBot.y()) * yfrac;
658 
659  double tx = myTop.x();
660  double ty = myTop.y();
661  double bx = myBot.x();
662  double by = myBot.y();
663  double mySrcX = bx + ( tx - bx ) * yfrac;
664  double mySrcY = by + ( ty - by ) * yfrac;
665 
666  if ( !mExtent.contains( QgsPoint( mySrcX, mySrcY ) ) )
667  {
668  return false;
669  }
670 
671  // TODO: check again cell selection (coor is in the middle)
672 
673  *theSrcRow = static_cast< int >( floor(( mSrcExtent.yMaximum() - mySrcY ) / mSrcYRes ) );
674  *theSrcCol = static_cast< int >( floor(( mySrcX - mSrcExtent.xMinimum() ) / mSrcXRes ) );
675 
676  // For now silently correct limits to avoid crashes
677  // TODO: review
678  // should not happen
679  if ( *theSrcRow >= mSrcRows ) return false;
680  if ( *theSrcRow < 0 ) return false;
681  if ( *theSrcCol >= mSrcCols ) return false;
682  if ( *theSrcCol < 0 ) return false;
683 
684  return true;
685 }
686 
687 void QgsRasterProjector::insertRows( const QgsCoordinateTransform* ct )
688 {
689  for ( int r = 0; r < mCPRows - 1; r++ )
690  {
691  QList<QgsPoint> myRow;
692  QList<bool> myLegalRow;
693  myRow.reserve( mCPCols );
694  myLegalRow.reserve( mCPCols );
695  for ( int c = 0; c < mCPCols; ++c )
696  {
697  myRow.append( QgsPoint() );
698  myLegalRow.append( false );
699  }
700  QgsDebugMsgLevel( QString( "insert new row at %1" ).arg( 1 + r*2 ), 3 );
701  mCPMatrix.insert( 1 + r*2, myRow );
702  mCPLegalMatrix.insert( 1 + r*2, myLegalRow );
703  }
704  mCPRows += mCPRows - 1;
705  for ( int r = 1; r < mCPRows - 1; r += 2 )
706  {
707  calcRow( r, ct );
708  }
709 }
710 
711 void QgsRasterProjector::insertCols( const QgsCoordinateTransform* ct )
712 {
713  for ( int r = 0; r < mCPRows; r++ )
714  {
715  for ( int c = 0; c < mCPCols - 1; c++ )
716  {
717  mCPMatrix[r].insert( 1 + c*2, QgsPoint() );
718  mCPLegalMatrix[r].insert( 1 + c*2, false );
719  }
720  }
721  mCPCols += mCPCols - 1;
722  for ( int c = 1; c < mCPCols - 1; c += 2 )
723  {
724  calcCol( c, ct );
725  }
726 
727 }
728 
729 void QgsRasterProjector::calcCP( int theRow, int theCol, const QgsCoordinateTransform* ct )
730 {
731  double myDestX, myDestY;
732  destPointOnCPMatrix( theRow, theCol, &myDestX, &myDestY );
733  QgsPoint myDestPoint( myDestX, myDestY );
734  try
735  {
736  if ( ct )
737  {
738  mCPMatrix[theRow][theCol] = ct->transform( myDestPoint );
739  mCPLegalMatrix[theRow][theCol] = true;
740  }
741  else
742  {
743  mCPLegalMatrix[theRow][theCol] = false;
744  }
745  }
746  catch ( QgsCsException &e )
747  {
748  Q_UNUSED( e );
749  // Caught an error in transform
750  mCPLegalMatrix[theRow][theCol] = false;
751  }
752 }
753 
754 bool QgsRasterProjector::calcRow( int theRow, const QgsCoordinateTransform* ct )
755 {
756  QgsDebugMsgLevel( QString( "theRow = %1" ).arg( theRow ), 3 );
757  for ( int i = 0; i < mCPCols; i++ )
758  {
759  calcCP( theRow, i, ct );
760  }
761 
762  return true;
763 }
764 
765 bool QgsRasterProjector::calcCol( int theCol, const QgsCoordinateTransform* ct )
766 {
767  QgsDebugMsgLevel( QString( "theCol = %1" ).arg( theCol ), 3 );
768  for ( int i = 0; i < mCPRows; i++ )
769  {
770  calcCP( i, theCol, ct );
771  }
772 
773  return true;
774 }
775 
776 bool QgsRasterProjector::checkCols( const QgsCoordinateTransform* ct )
777 {
778  if ( !ct )
779  {
780  return false;
781  }
782 
783  for ( int c = 0; c < mCPCols; c++ )
784  {
785  for ( int r = 1; r < mCPRows - 1; r += 2 )
786  {
787  double myDestX, myDestY;
788  destPointOnCPMatrix( r, c, &myDestX, &myDestY );
789  QgsPoint myDestPoint( myDestX, myDestY );
790 
791  QgsPoint mySrcPoint1 = mCPMatrix[r-1][c];
792  QgsPoint mySrcPoint2 = mCPMatrix[r][c];
793  QgsPoint mySrcPoint3 = mCPMatrix[r+1][c];
794 
795  QgsPoint mySrcApprox(( mySrcPoint1.x() + mySrcPoint3.x() ) / 2, ( mySrcPoint1.y() + mySrcPoint3.y() ) / 2 );
796  if ( !mCPLegalMatrix[r-1][c] || !mCPLegalMatrix[r][c] || !mCPLegalMatrix[r+1][c] )
797  {
798  // There was an error earlier in transform, just abort
799  return false;
800  }
801  try
802  {
803  QgsPoint myDestApprox = ct->transform( mySrcApprox, QgsCoordinateTransform::ReverseTransform );
804  double mySqrDist = myDestApprox.sqrDist( myDestPoint );
805  if ( mySqrDist > mSqrTolerance )
806  {
807  return false;
808  }
809  }
810  catch ( QgsCsException &e )
811  {
812  Q_UNUSED( e );
813  // Caught an error in transform
814  return false;
815  }
816  }
817  }
818  return true;
819 }
820 
821 bool QgsRasterProjector::checkRows( const QgsCoordinateTransform* ct )
822 {
823  if ( !ct )
824  {
825  return false;
826  }
827 
828  for ( int r = 0; r < mCPRows; r++ )
829  {
830  for ( int c = 1; c < mCPCols - 1; c += 2 )
831  {
832  double myDestX, myDestY;
833  destPointOnCPMatrix( r, c, &myDestX, &myDestY );
834 
835  QgsPoint myDestPoint( myDestX, myDestY );
836  QgsPoint mySrcPoint1 = mCPMatrix[r][c-1];
837  QgsPoint mySrcPoint2 = mCPMatrix[r][c];
838  QgsPoint mySrcPoint3 = mCPMatrix[r][c+1];
839 
840  QgsPoint mySrcApprox(( mySrcPoint1.x() + mySrcPoint3.x() ) / 2, ( mySrcPoint1.y() + mySrcPoint3.y() ) / 2 );
841  if ( !mCPLegalMatrix[r][c-1] || !mCPLegalMatrix[r][c] || !mCPLegalMatrix[r][c+1] )
842  {
843  // There was an error earlier in transform, just abort
844  return false;
845  }
846  try
847  {
848  QgsPoint myDestApprox = ct->transform( mySrcApprox, QgsCoordinateTransform::ReverseTransform );
849  double mySqrDist = myDestApprox.sqrDist( myDestPoint );
850  if ( mySqrDist > mSqrTolerance )
851  {
852  return false;
853  }
854  }
855  catch ( QgsCsException &e )
856  {
857  Q_UNUSED( e );
858  // Caught an error in transform
859  return false;
860  }
861  }
862  }
863  return true;
864 }
865 
867 {
868  switch ( precision )
869  {
870  case Approximate:
871  return tr( "Approximate" );
872  case Exact:
873  return tr( "Exact" );
874  }
875  return "Unknown";
876 }
877 
878 QgsRasterBlock * QgsRasterProjector::block( int bandNo, QgsRectangle const & extent, int width, int height )
879 {
880  QgsDebugMsgLevel( QString( "extent:\n%1" ).arg( extent.toString() ), 4 );
881  QgsDebugMsgLevel( QString( "width = %1 height = %2" ).arg( width ).arg( height ), 4 );
882  if ( !mInput )
883  {
884  QgsDebugMsgLevel( "Input not set", 4 );
885  return new QgsRasterBlock();
886  }
887 
888  if ( ! mSrcCRS.isValid() || ! mDestCRS.isValid() || mSrcCRS == mDestCRS )
889  {
890  QgsDebugMsgLevel( "No projection necessary", 4 );
891  return mInput->block( bandNo, extent, width, height );
892  }
893 
894  mDestExtent = extent;
895  mDestRows = height;
896  mDestCols = width;
897  calc();
898 
899  QgsDebugMsgLevel( QString( "srcExtent:\n%1" ).arg( srcExtent().toString() ), 4 );
900  QgsDebugMsgLevel( QString( "srcCols = %1 srcRows = %2" ).arg( srcCols() ).arg( srcRows() ), 4 );
901 
902  // If we zoom out too much, projector srcRows / srcCols maybe 0, which can cause problems in providers
903  if ( srcRows() <= 0 || srcCols() <= 0 )
904  {
905  QgsDebugMsgLevel( "Zero srcRows or srcCols", 4 );
906  return new QgsRasterBlock();
907  }
908 
909  QgsRasterBlock *inputBlock = mInput->block( bandNo, srcExtent(), srcCols(), srcRows() );
910  if ( !inputBlock || inputBlock->isEmpty() )
911  {
912  QgsDebugMsg( "No raster data!" );
913  delete inputBlock;
914  return new QgsRasterBlock();
915  }
916 
917  qgssize pixelSize = QgsRasterBlock::typeSize( mInput->dataType( bandNo ) );
918 
919  QgsRasterBlock *outputBlock;
920  if ( inputBlock->hasNoDataValue() )
921  {
922  outputBlock = new QgsRasterBlock( inputBlock->dataType(), width, height, inputBlock->noDataValue() );
923  }
924  else
925  {
926  outputBlock = new QgsRasterBlock( inputBlock->dataType(), width, height );
927  }
928  if ( !outputBlock->isValid() )
929  {
930  QgsDebugMsg( "Cannot create block" );
931  delete inputBlock;
932  return outputBlock;
933  }
934 
935  // set output to no data, it should be fast
936  outputBlock->setIsNoData();
937 
938  // No data: because isNoData()/setIsNoData() is slow with respect to simple memcpy,
939  // we use if only if necessary:
940  // 1) no data value exists (numerical) -> memcpy, not necessary isNoData()/setIsNoData()
941  // 2) no data value does not exist but it may contain no data (numerical no data bitmap)
942  // -> must use isNoData()/setIsNoData()
943  // 3) no data are not used (no no data value, no no data bitmap) -> simple memcpy
944  // 4) image - simple memcpy
945 
946  // To copy no data values stored in bitmaps we have to use isNoData()/setIsNoData(),
947  // we cannot fill output block with no data because we use memcpy for data, not setValue().
948  bool doNoData = !QgsRasterBlock::typeIsNumeric( inputBlock->dataType() ) && inputBlock->hasNoData() && !inputBlock->hasNoDataValue();
949 
950  const QgsCoordinateTransform* inverseCt = nullptr;
951  if ( !mApproximate )
952  {
953  inverseCt = QgsCoordinateTransformCache::instance()->transform( mDestCRS.authid(), mSrcCRS.authid(), mDestDatumTransform, mSrcDatumTransform );
954  }
955 
956  outputBlock->setIsNoData();
957 
958  int srcRow, srcCol;
959  for ( int i = 0; i < height; ++i )
960  {
961  for ( int j = 0; j < width; ++j )
962  {
963  bool inside = srcRowCol( i, j, &srcRow, &srcCol, inverseCt );
964  if ( !inside ) continue; // we have everything set to no data
965 
966  qgssize srcIndex = static_cast< qgssize >( srcRow ) * mSrcCols + srcCol;
967  QgsDebugMsgLevel( QString( "row = %1 col = %2 srcRow = %3 srcCol = %4" ).arg( i ).arg( j ).arg( srcRow ).arg( srcCol ), 5 );
968 
969  // isNoData() may be slow so we check doNoData first
970  if ( doNoData && inputBlock->isNoData( srcRow, srcCol ) )
971  {
972  outputBlock->setIsNoData( i, j );
973  continue;
974  }
975 
976  qgssize destIndex = static_cast< qgssize >( i ) * width + j;
977  char *srcBits = inputBlock->bits( srcIndex );
978  char *destBits = outputBlock->bits( destIndex );
979  if ( !srcBits )
980  {
981  QgsDebugMsg( QString( "Cannot get input block data: row = %1 col = %2" ).arg( i ).arg( j ) );
982  continue;
983  }
984  if ( !destBits )
985  {
986  QgsDebugMsg( QString( "Cannot set output block data: srcRow = %1 srcCol = %2" ).arg( srcRow ).arg( srcCol ) );
987  continue;
988  }
989  memcpy( destBits, srcBits, pixelSize );
990  outputBlock->setIsData( i, j );
991  }
992  }
993 
994  delete inputBlock;
995 
996  return outputBlock;
997 }
998 
999 bool QgsRasterProjector::destExtentSize( const QgsRectangle& theSrcExtent, int theSrcXSize, int theSrcYSize,
1000  QgsRectangle& theDestExtent, int& theDestXSize, int& theDestYSize )
1001 {
1002  if ( theSrcExtent.isEmpty() || theSrcXSize <= 0 || theSrcYSize <= 0 )
1003  {
1004  return false;
1005  }
1006  const QgsCoordinateTransform* ct = QgsCoordinateTransformCache::instance()->transform( mSrcCRS.authid(), mDestCRS.authid(), mSrcDatumTransform, mDestDatumTransform );
1007 
1008  return extentSize( ct, theSrcExtent, theSrcXSize, theSrcYSize, theDestExtent, theDestXSize, theDestYSize );
1009 }
1010 
1012  const QgsRectangle& theSrcExtent, int theSrcXSize, int theSrcYSize,
1013  QgsRectangle& theDestExtent, int& theDestXSize, int& theDestYSize )
1014 {
1015  if ( theSrcExtent.isEmpty() || theSrcXSize <= 0 || theSrcYSize <= 0 )
1016  {
1017  return false;
1018  }
1019 
1020  theDestExtent = ct->transformBoundingBox( theSrcExtent );
1021 
1022  // We reproject pixel rectangle from 9 points matrix of source extent, of course, it gives
1023  // bigger xRes,yRes than reprojected edges (envelope)
1024  double srcXStep = theSrcExtent.width() / 3;
1025  double srcYStep = theSrcExtent.height() / 3;
1026  double srcXRes = theSrcExtent.width() / theSrcXSize;
1027  double srcYRes = theSrcExtent.height() / theSrcYSize;
1028  double destXRes = std::numeric_limits<double>::max();
1029  double destYRes = std::numeric_limits<double>::max();
1030 
1031  for ( int i = 0; i < 3; i++ )
1032  {
1033  double x = theSrcExtent.xMinimum() + i * srcXStep;
1034  for ( int j = 0; j < 3; j++ )
1035  {
1036  double y = theSrcExtent.yMinimum() + j * srcYStep;
1037  QgsRectangle srcRectangle( x - srcXRes / 2, y - srcYRes / 2, x + srcXRes / 2, y + srcYRes / 2 );
1038  QgsRectangle destRectangle = ct->transformBoundingBox( srcRectangle );
1039  if ( destRectangle.width() > 0 )
1040  {
1041  destXRes = std::min( destXRes, destRectangle.width() );
1042  }
1043  if ( destRectangle.height() > 0 )
1044  {
1045  destYRes = std::min( destYRes, destRectangle.height() );
1046  }
1047  }
1048  }
1049  theDestXSize = std::max( 1, static_cast< int >( theDestExtent.width() / destYRes ) );
1050  theDestYSize = std::max( 1, static_cast< int >( theDestExtent.height() / destYRes ) );
1051 
1052  return true;
1053 }
1054 
virtual int bandCount() const =0
Get number of bands.
void clear()
A rectangle specified with double values.
Definition: qgsrectangle.h:35
Precision precision() const
bool isEmpty() const
test if rectangle is empty.
void setCRS(const QgsCoordinateReferenceSystem &theSrcCRS, const QgsCoordinateReferenceSystem &theDestCRS, int srcDatumTransform=-1, int destDatumTransform=-1)
set source and destination CRS
Unknown or unspecified type.
Definition: qgis.h:135
void setXMaximum(double x)
Set the maximum x value.
Definition: qgsrectangle.h:172
Approximate (default), fast but possibly inaccurate.
double yMaximum() const
Get the y maximum value (top side of rectangle)
Definition: qgsrectangle.h:197
static bool typeIsNumeric(QGis::DataType type)
Returns true if data type is numeric.
#define QgsDebugMsg(str)
Definition: qgslogger.h:33
void reserve(int alloc)
bool contains(const QgsRectangle &rect) const
return true when rectangle contains other rectangle
virtual const QgsRasterInterface * srcInput() const
Get source / raw input, the first in pipe, usually provider.
double noDataValue() const
Return no data value.
QGis::DataType dataType() const
Returns data type.
QgsPoint transform(const QgsPoint &p, TransformDirection direction=ForwardTransform) const
Transform the point from Source Coordinate System to Destination Coordinate System If the direction i...
static QgsCoordinateTransformCache * instance()
Definition: qgscrscache.cpp:22
double sqrDist(double x, double y) const
Returns the squared distance between this point a specified x, y coordinate.
Definition: qgspoint.cpp:353
bool isNoData(int row, int column)
Check if value at position is no data.
double x() const
Get the x value of the point.
Definition: qgspoint.h:185
virtual int ySize() const
bool setIsNoData(int row, int column)
Set no data on pixel.
double ANALYSIS_EXPORT max(double x, double y)
Returns the maximum of two doubles or the first argument if both are equal.
const QgsCoordinateTransform * transform(const QString &srcAuthId, const QString &destAuthId, int srcDatumTransform=-1, int destDatumTransform=-1)
Returns coordinate transformation.
Definition: qgscrscache.cpp:41
void transformInPlace(double &x, double &y, double &z, TransformDirection direction=ForwardTransform) const
void append(const T &value)
~QgsRasterProjector()
The destructor.
bool hasNoData() const
Returns true if the block may contain no data.
Raster data container.
double yMinimum() const
Get the y minimum value (bottom side of rectangle)
Definition: qgsrectangle.h:202
QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height) override
Read block of data using given extent and size.
double xMaximum() const
Get the x maximum value (right side of rectangle)
Definition: qgsrectangle.h:187
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:34
void setYMinimum(double y)
Set the minimum y value.
Definition: qgsrectangle.h:177
QString toString() const
String representation of the point (x,y)
Definition: qgspoint.cpp:134
static int typeSize(int dataType)
virtual QGis::DataType dataType(int bandNo) const =0
Returns data type for the band specified by number.
Base class for processing filters like renderers, reprojector, resampler etc.
A class to represent a point.
Definition: qgspoint.h:117
unsigned long long qgssize
Qgssize is used instead of size_t, because size_t is stdlib type, unknown by SIP, and it would be har...
Definition: qgis.h:500
void setX(double x)
Sets the x value of the point.
Definition: qgspoint.h:162
virtual int capabilities() const
Returns a bitmask containing the supported capabilities.
void setY(double y)
Sets the y value of the point.
Definition: qgspoint.h:170
void combineExtentWith(const QgsRectangle &rect)
expand the rectangle so that covers both the original rectangle and the given rectangle ...
virtual QgsRectangle extent() override=0
Get the extent of the data source.
virtual QgsRectangle extent()
Get the extent of the interface.
bool isValid() const
Returns whether this CRS is correctly initialized and usable.
bool hasNoDataValue() const
True if the block has no data value.
char * bits(int row, int column)
Get pointer to data.
QgsRasterProjector * clone() const override
Clone itself, create deep copy.
Precision
Precision defines if each pixel is reprojected or approximate reprojection based on an approximation ...
virtual int xSize() const
Get raster size.
void insert(int i, const T &value)
void setYMaximum(double y)
Set the maximum y value.
Definition: qgsrectangle.h:182
QgsRectangle intersect(const QgsRectangle *rect) const
return the intersection with the given rectangle
Class for storing a coordinate reference system (CRS)
DataType
Raster data types.
Definition: qgis.h:133
virtual QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height)=0
Read block of data using given extent and size.
QString authid() const
Returns the authority identifier for the CRS, which includes both the authority (eg EPSG) and the CRS...
Class for doing transforms between two map coordinate systems.
static bool extentSize(const QgsCoordinateTransform *ct, const QgsRectangle &theSrcExtent, int theSrcXSize, int theSrcYSize, QgsRectangle &theDestExtent, int &theDestXSize, int &theDestYSize)
Calculate destination extent and size from source extent and size.
double y() const
Get the y value of the point.
Definition: qgspoint.h:193
Exact, precise but slow.
Custom exception class for Coordinate Reference System related exceptions.
QGis::DataType dataType(int bandNo) const override
Returns data type for the band specified by number.
double ANALYSIS_EXPORT min(double x, double y)
Returns the minimum of two doubles or the first argument if both are equal.
double width() const
Width of the rectangle.
Definition: qgsrectangle.h:207
bool destExtentSize(const QgsRectangle &theSrcExtent, int theSrcXSize, int theSrcYSize, QgsRectangle &theDestExtent, int &theDestXSize, int &theDestYSize)
Calculate destination extent and size from source extent and size.
QgsRasterInterface * mInput
QgsRasterProjector & operator=(const QgsRasterProjector &projector)
QString toString(bool automaticPrecision=false) const
returns string representation of form xmin,ymin xmax,ymax
QString arg(qlonglong a, int fieldWidth, int base, const QChar &fillChar) const
double xMinimum() const
Get the x minimum value (left side of rectangle)
Definition: qgsrectangle.h:192
static QString precisionLabel(Precision precision)
void setXMinimum(double x)
Set the minimum x value.
Definition: qgsrectangle.h:167
QgsRectangle transformBoundingBox(const QgsRectangle &theRect, TransformDirection direction=ForwardTransform, const bool handle180Crossover=false) const
Transform a QgsRectangle to the dest Coordinate system If the direction is ForwardTransform then coor...
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:212
bool isEmpty() const
Returns true if block is empty, i.e.
int bandCount() const override
Get number of bands.
Base class for raster data providers.