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