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