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