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