align extent according to src resolution/origin

git-svn-id: http://svn.osgeo.org/qgis/trunk/qgis@15533 c8812cc2-4d05-0410-92ff-de0c093fc19c

Author: rblazek

src/core/qgsrasterdataprovider.cpp

@@ -47,7 +47,7 @@ void QgsRasterDataProvider::readBlock( int bandNo, QgsRectangle  const & viewExt
     mMaxSrcYRes = extent().height() / ySize();
   }
 
-  QgsRasterProjector myProjector( theSrcCRS, theDestCRS, viewExtent, height, width, mMaxSrcXRes, mMaxSrcYRes );
+  QgsRasterProjector myProjector( theSrcCRS, theDestCRS, viewExtent, height, width, mMaxSrcXRes, mMaxSrcYRes, mExtent );
 
   QgsDebugMsg( QString( "create projector time  (ms): %1" ).arg( time.elapsed() ) );
 
@@ -182,10 +182,10 @@ QString QgsRasterDataProvider::lastErrorFormat()
 
 QByteArray QgsRasterDataProvider::noValueBytes( int theBandNo )
 {
-  int type = dataType(theBandNo);
-  int size = dataTypeSize(theBandNo) / 8;
+  int type = dataType( theBandNo );
+  int size = dataTypeSize( theBandNo ) / 8;
   QByteArray ba;
-  ba.resize(size);
+  ba.resize( size );
   char * data = ba.data();
   double noval = mNoDataValue[theBandNo-1];
   unsigned char uc;
@@ -195,35 +195,35 @@ QByteArray QgsRasterDataProvider::noValueBytes( int theBandNo )
   int i;
   float f;
   double d;
-  switch (type) 
+  switch ( type )
   {
     case QgsRasterDataProvider::Byte:
-      uc = (unsigned char)noval;
-      memcpy ( data, &uc, size);
+      uc = ( unsigned char )noval;
+      memcpy( data, &uc, size );
       break;
     case QgsRasterDataProvider::UInt16:
-      us = (unsigned short)noval;
-      memcpy ( data, &us, size);
+      us = ( unsigned short )noval;
+      memcpy( data, &us, size );
       break;
     case QgsRasterDataProvider::Int16:
-      s = (short)noval;
-      memcpy ( data, &s, size);
+      s = ( short )noval;
+      memcpy( data, &s, size );
       break;
     case QgsRasterDataProvider::UInt32:
-      ui = (unsigned int)noval;
-      memcpy ( data, &ui, size);
+      ui = ( unsigned int )noval;
+      memcpy( data, &ui, size );
       break;
     case QgsRasterDataProvider::Int32:
-      i = (int)noval;
-      memcpy ( data, &i, size);
+      i = ( int )noval;
+      memcpy( data, &i, size );
       break;
     case QgsRasterDataProvider::Float32:
-      f = (float)noval;
-      memcpy ( data, &f, size);
+      f = ( float )noval;
+      memcpy( data, &f, size );
       break;
     case QgsRasterDataProvider::Float64:
-      d = (double)noval;
-      memcpy ( data, &d, size);
+      d = ( double )noval;
+      memcpy( data, &d, size );
       break;
     default:
       QgsLogger::warning( "GDAL data type is not supported" );

src/core/qgsrasterprojector.cpp

@@ -27,12 +27,14 @@ QgsRasterProjector::QgsRasterProjector(
   QgsCoordinateReferenceSystem theDestCRS,
   QgsRectangle theDestExtent,
   int theDestRows, int theDestCols,
-  double theMaxSrcXRes, double theMaxSrcYRes )
+  double theMaxSrcXRes, double theMaxSrcYRes,
+  QgsRectangle theExtent )
     : mSrcCRS( theSrcCRS )
     , mDestCRS( theDestCRS )
     , mDestExtent( theDestExtent )
     , mDestRows( theDestRows ), mDestCols( theDestCols )
     , mMaxSrcXRes( theMaxSrcXRes ), mMaxSrcYRes( theMaxSrcYRes )
+    , mExtent( theExtent )
 {
   QgsDebugMsg( "Entered" );
   QgsDebugMsg( "theDestExtent = " + theDestExtent.toString() );
@@ -135,7 +137,37 @@ void QgsRasterProjector::calcSrcExtent()
     myPoint = mCPMatrix[mCPRows-1][i];
     mSrcExtent.combineExtentWith( myPoint.x(), myPoint.y() );
   }
-  // Expand a bit to avoid possible approx coords falling out because of representation error
+  // Expand a bit to avoid possible approx coords falling out because of representation error?
+
+  // If mMaxSrcXRes, mMaxSrcYRes are defined (fixed src resolution)
+  // align extent to src resolution to avoid jumping reprojected pixels
+  // because of shifting resampled grid
+  // Important especially if we are over mMaxSrcXRes, mMaxSrcYRes limits
+
+  QgsDebugMsg( "mSrcExtent = " + mSrcExtent.toString() );
+
+  if ( mMaxSrcXRes > 0 )
+  {
+    // with floor/ceil it should work correctly also for mSrcExtent.xMinimum() < mExtent.xMinimum()
+    double col = floor(( mSrcExtent.xMinimum() - mExtent.xMinimum() ) / mMaxSrcXRes );
+    double x = mExtent.xMinimum() + col * mMaxSrcXRes;
+    mSrcExtent.setXMinimum( x );
+
+    col = ceil(( mSrcExtent.xMaximum() - mExtent.xMinimum() ) / mMaxSrcXRes );
+    x = mExtent.xMinimum() + col * mMaxSrcXRes;
+    mSrcExtent.setXMaximum( x );
+  }
+  if ( mMaxSrcYRes > 0 )
+  {
+    double row = floor(( mExtent.yMaximum() - mSrcExtent.yMaximum() ) / mMaxSrcYRes );
+    double y = mExtent.yMaximum() - row * mMaxSrcYRes;
+    mSrcExtent.setYMaximum( y );
+
+    row = ceil(( mExtent.yMaximum() - mSrcExtent.yMinimum() ) / mMaxSrcYRes );
+    y = mExtent.yMaximum() - row * mMaxSrcYRes;
+    mSrcExtent.setYMinimum( y );
+  }
+
 
   QgsDebugMsg( "mSrcExtent = " + mSrcExtent.toString() );
 }
@@ -197,8 +229,10 @@ void QgsRasterProjector::calcSrcRowsCols()
   double myMinYSize = mMaxSrcYRes > myMinSize ? mMaxSrcYRes : myMinSize;
   QgsDebugMsg( QString( "myMinXSize = %1 myMinYSize = %2" ).arg( myMinXSize ).arg( myMinYSize ) );
   QgsDebugMsg( QString( "mSrcExtent.width = %1 mSrcExtent.height = %2" ).arg( mSrcExtent.width() ).arg( mSrcExtent.height() ) );
-  mSrcRows = ( int ) ceil( mSrcExtent.height() / myMinYSize );
-  mSrcCols = ( int ) ceil( mSrcExtent.width() / myMinXSize );
+
+  // we have to round to keep alignment set in calcSrcExtent
+  mSrcRows = ( int ) qRound( mSrcExtent.height() / myMinYSize );
+  mSrcCols = ( int ) qRound( mSrcExtent.width() / myMinXSize );
 
   QgsDebugMsg( QString( "mSrcRows = %1 mSrcCols = %2" ).arg( mSrcRows ).arg( mSrcCols ) );
 }

src/core/qgsrasterprojector.h

@@ -51,7 +51,8 @@ class QgsRasterProjector
       QgsCoordinateReferenceSystem theDestCRS,
       QgsRectangle theDestExtent,
       int theDestRows, int theDestCols,
-      double theMaxSrcXRes, double theMaxSrcYRes
+      double theMaxSrcXRes, double theMaxSrcYRes,
+      QgsRectangle theExtent
     );
 
     /** \brief The destructor */
@@ -136,6 +137,9 @@ class QgsRasterProjector
     /** Source extent */
     QgsRectangle mSrcExtent;
 
+    /** Source raster extent */
+    QgsRectangle mExtent;
+
     /** Number of destination rows */
     int mDestRows;