QGIS API Documentation  2.7.0-Master
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules 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;
102  mSrcDatumTransform = projector.mSrcDatumTransform;
103  mDestDatumTransform = projector.mDestDatumTransform;
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;
115  mSrcDatumTransform = projector.mSrcDatumTransform;
116  mDestDatumTransform = projector.mDestDatumTransform;
117  mMaxSrcXRes = projector.mMaxSrcXRes;
118  mMaxSrcYRes = projector.mMaxSrcYRes;
119  mExtent = projector.mExtent;
120  }
121  return *this;
122 }
123 
125 {
126  QgsDebugMsg( "Entered" );
127  QgsRasterProjector * projector = new QgsRasterProjector( mSrcCRS, mDestCRS, mMaxSrcXRes, mMaxSrcYRes, mExtent );
128  projector->mSrcDatumTransform = mSrcDatumTransform;
129  projector->mDestDatumTransform = mDestDatumTransform;
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 
161 void QgsRasterProjector::calc()
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 
189  mDestXRes = mDestExtent.width() / ( mDestCols );
190  mDestYRes = mDestExtent.height() / ( mDestRows );
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 
199  const QgsCoordinateTransform* ct = QgsCoordinateTransformCache::instance()->transform( mDestCRS.authid(), mSrcCRS.authid(), mDestDatumTransform, mSrcDatumTransform );
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 );
254  QgsDebugMsgLevel( cpToString(), 5 );
255 
256  // Calculate source dimensions
257  calcSrcExtent();
258  calcSrcRowsCols();
259  mSrcYRes = mSrcExtent.height() / mSrcRows;
260  mSrcXRes = mSrcExtent.width() / mSrcCols;
261 
262  // init helper points
263  pHelperTop = new QgsPoint[mDestCols];
264  pHelperBottom = new QgsPoint[mDestCols];
265  calcHelper( 0, pHelperTop );
266  calcHelper( 1, pHelperBottom );
267  mHelperTopRow = 0;
268 }
269 
270 void QgsRasterProjector::calcSrcExtent()
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
294  mSrcExtent = mSrcExtent.intersect( &mExtent );
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 
332 QString QgsRasterProjector::cpToString()
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 
357 void QgsRasterProjector::calcSrcRowsCols()
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 }
459 void QgsRasterProjector::nextHelper()
460 {
461  // We just switch pHelperTop and pHelperBottom, memory is not lost
462  QgsPoint *tmp;
463  tmp = pHelperTop;
464  pHelperTop = pHelperBottom;
465  pHelperBottom = tmp;
466  calcHelper( mHelperTopRow + 2, pHelperBottom );
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 
588 void QgsRasterProjector::insertRows( const QgsCoordinateTransform* ct )
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 
610 void QgsRasterProjector::insertCols( const QgsCoordinateTransform* ct )
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 
655 bool QgsRasterProjector::calcRow( int theRow, const QgsCoordinateTransform* ct )
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 
666 bool QgsRasterProjector::calcCol( int theCol, const QgsCoordinateTransform* ct )
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 
677 bool QgsRasterProjector::checkCols( const QgsCoordinateTransform* ct )
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 
722 bool QgsRasterProjector::checkRows( const QgsCoordinateTransform* ct )
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 
783  mDestExtent = extent;
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  {
842  ct = QgsCoordinateTransformCache::instance()->transform( mDestCRS.authid(), mSrcCRS.authid(), mDestDatumTransform, mSrcDatumTransform );
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.
A rectangle specified with double values.
Definition: qgsrectangle.h:35
QgsRasterInterface * clone() const
Clone itself, create deep copy.
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:163
double yMaximum() const
Get the y maximum value (top side of rectangle)
Definition: qgsrectangle.h:188
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.
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:193
double xMaximum() const
Get the x maximum value (right side of rectangle)
Definition: qgsrectangle.h:178
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:34
virtual QgsRectangle extent()=0
Get the extent of the data source.
void setYMinimum(double y)
Set the minimum y value.
Definition: qgsrectangle.h:168
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.
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:420
void setX(double x)
Definition: qgspoint.h:103
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:111
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.
QGis::DataType dataType(int bandNo) const
Returns data type for the band specified by number.
void setYMaximum(double y)
Set the maximum y value.
Definition: qgsrectangle.h:173
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.
double width() const
Width of the rectangle.
Definition: qgsrectangle.h:198
QgsRasterInterface * mInput
QgsRasterBlock * block(int bandNo, const QgsRectangle &extent, int width, int height)
Read block of data using given extent and size.
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:183
void setXMinimum(double x)
Set the minimum x value.
Definition: qgsrectangle.h:158
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:203
bool isEmpty() const
Returns true if block is empty, i.e.
Base class for raster data providers.