Skip to content
Permalink
Browse files

align extent according to src resolution/origin

git-svn-id: http://svn.osgeo.org/qgis/trunk@15533 c8812cc2-4d05-0410-92ff-de0c093fc19c
  • Loading branch information
rblazek
rblazek committed Mar 18, 2011
1 parent bf2f284 commit 44cc8a73a439f6fb2fb3a155c2d06ca6b8a7a421
Showing with 62 additions and 24 deletions.
  1. +19 −19 src/core/qgsrasterdataprovider.cpp
  2. +38 −4 src/core/qgsrasterprojector.cpp
  3. +5 −1 src/core/qgsrasterprojector.h
@@ -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" );
@@ -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 ) );
}
@@ -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;

0 comments on commit 44cc8a7

Please sign in to comment.
You can’t perform that action at this time.