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