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