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