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