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