QGIS API Documentation  2.99.0-Master (9fdd060)
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
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:264
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:131
void setY(double y)
Sets the y value of the point.
Definition: qgspointxy.h:114
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:105
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:119
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.
QString toString() const
String representation of the point (x,y)
Definition: qgspointxy.cpp:45
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:109
Exact, precise but slow.
Custom exception class for Coordinate Reference System related exceptions.
Definition: qgsexception.h:66
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:138
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.