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