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