QGIS API Documentation  2.99.0-Master (314842d)
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<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 }
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  QgsPoint 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  QgsPoint 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  QgsPoint myPointA = mCPMatrix[i][j];
326  QgsPoint myPointB = mCPMatrix[i][j + 1];
327  QgsPoint 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, QgsPoint *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  QgsPoint &mySrcPoint0 = mCPMatrix[matrixRow][myMatrixCol];
411  QgsPoint &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  QgsPoint *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( QgsPoint( 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  QgsPoint &myTop = pHelperTop[destCol];
515  QgsPoint &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( QgsPoint( 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<QgsPoint> 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( QgsPoint() );
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, QgsPoint() );
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  QgsPoint 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  QgsPoint myDestPoint( myDestX, myDestY );
652 
653  QgsPoint mySrcPoint1 = mCPMatrix[r - 1][c];
654  QgsPoint mySrcPoint2 = mCPMatrix[r][c];
655  QgsPoint mySrcPoint3 = mCPMatrix[r + 1][c];
656 
657  QgsPoint 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  QgsPoint 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  QgsPoint myDestPoint( myDestX, myDestY );
698  QgsPoint mySrcPoint1 = mCPMatrix[r][c - 1];
699  QgsPoint mySrcPoint2 = mCPMatrix[r][c];
700  QgsPoint mySrcPoint3 = mCPMatrix[r][c + 1];
701 
702  QgsPoint 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  QgsPoint 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:36
double y
Definition: qgspoint.h:42
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:36
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: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
void setCrs(const QgsCoordinateReferenceSystem &srcCRS, const QgsCoordinateReferenceSystem &destCRS, int srcDatumTransform=-1, int destDatumTransform=-1)
set source and destination CRS
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:215
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:63
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.
A class to represent a point.
Definition: qgspoint.h:37
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:45
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:349
void setX(double x)
Sets the x value of the point.
Definition: qgspoint.h:95
void setY(double y)
Sets the y value of the point.
Definition: qgspoint.h:103
double yMinimum() const
Get the y minimum value (bottom side of rectangle)
Definition: qgsrectangle.h:210
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.
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:62
Transform from destination to source CRS.
This class represents a coordinate reference system (CRS).
Class for doing transforms between two map coordinate systems.
double xMinimum() const
Get the x minimum value (left side of rectangle)
Definition: qgsrectangle.h:200
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:264
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
Height of the rectangle.
Definition: qgsrectangle.h:220
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:41