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