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