128 changes: 93 additions & 35 deletions src/core/raster/qgsrasterprojector.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,21 +16,25 @@
***************************************************************************/

#include "qgsrasterdataprovider.h"
#include "qgscrscache.h"
#include "qgslogger.h"
#include "qgsrasterprojector.h"
#include "qgscoordinatetransform.h"

QgsRasterProjector::QgsRasterProjector(
QgsCoordinateReferenceSystem theSrcCRS,
QgsCoordinateReferenceSystem theDestCRS,
int theSrcDatumTransform,
int theDestDatumTransform,
QgsRectangle theDestExtent,
int theDestRows, int theDestCols,
double theMaxSrcXRes, double theMaxSrcYRes,
QgsRectangle theExtent )
: QgsRasterInterface( 0 )
, mSrcCRS( theSrcCRS )
, mDestCRS( theDestCRS )
, mCoordinateTransform( theDestCRS, theSrcCRS )
, mSrcDatumTransform( theSrcDatumTransform )
, mDestDatumTransform( theDestDatumTransform )
, mDestExtent( theDestExtent )
, mExtent( theExtent )
, mDestRows( theDestRows ), mDestCols( theDestCols )
Expand All @@ -46,22 +50,46 @@ QgsRasterProjector::QgsRasterProjector(
QgsRasterProjector::QgsRasterProjector(
QgsCoordinateReferenceSystem theSrcCRS,
QgsCoordinateReferenceSystem theDestCRS,
QgsRectangle theDestExtent,
int theDestRows, int theDestCols,
double theMaxSrcXRes, double theMaxSrcYRes,
QgsRectangle theExtent )
: QgsRasterInterface( 0 )
, mSrcCRS( theSrcCRS )
, mDestCRS( theDestCRS )
, mCoordinateTransform( theDestCRS, theSrcCRS )
, mSrcDatumTransform( -1 )
, mDestDatumTransform( -1 )
, mDestExtent( theDestExtent )
, mExtent( theExtent )
, mDestRows( theDestRows ), mDestCols( theDestCols )
, pHelperTop( 0 ), pHelperBottom( 0 )
, mMaxSrcXRes( theMaxSrcXRes ), mMaxSrcYRes( theMaxSrcYRes )
{
QgsDebugMsg( "Entered" );
QgsDebugMsg( "theDestExtent = " + theDestExtent.toString() );

calc();
}

QgsRasterProjector::QgsRasterProjector()
QgsRasterProjector::QgsRasterProjector(
QgsCoordinateReferenceSystem theSrcCRS,
QgsCoordinateReferenceSystem theDestCRS,
double theMaxSrcXRes, double theMaxSrcYRes,
QgsRectangle theExtent )
: QgsRasterInterface( 0 )
, mSrcCRS( theSrcCRS )
, mDestCRS( theDestCRS )
, mSrcDatumTransform( -1 )
, mDestDatumTransform( -1 )
, mExtent( theExtent )
, pHelperTop( 0 ), pHelperBottom( 0 )
, mMaxSrcXRes( theMaxSrcXRes ), mMaxSrcYRes( theMaxSrcYRes )
{
QgsDebugMsg( "Entered" );
}

QgsRasterProjector::QgsRasterProjector()
: QgsRasterInterface( 0 ), mSrcDatumTransform( -1 ), mDestDatumTransform( -1 ) , pHelperTop( 0 ), pHelperBottom( 0 )
{
QgsDebugMsg( "Entered" );
}
Expand All @@ -71,11 +99,11 @@ QgsRasterProjector::QgsRasterProjector( const QgsRasterProjector &projector )
{
mSrcCRS = projector.mSrcCRS;
mDestCRS = projector.mDestCRS;
mSrcDatumTransform = projector.mSrcDatumTransform;
mDestDatumTransform = projector.mDestDatumTransform;
mMaxSrcXRes = projector.mMaxSrcXRes;
mMaxSrcYRes = projector.mMaxSrcYRes;
mExtent = projector.mExtent;
mCoordinateTransform.setSourceCrs( mSrcCRS );
mCoordinateTransform.setDestCRS( mDestCRS );
}

QgsRasterProjector & QgsRasterProjector::operator=( const QgsRasterProjector & projector )
Expand All @@ -84,11 +112,11 @@ QgsRasterProjector & QgsRasterProjector::operator=( const QgsRasterProjector & p
{
mSrcCRS = projector.mSrcCRS;
mDestCRS = projector.mDestCRS;
mSrcDatumTransform = projector.mSrcDatumTransform;
mDestDatumTransform = projector.mDestDatumTransform;
mMaxSrcXRes = projector.mMaxSrcXRes;
mMaxSrcYRes = projector.mMaxSrcYRes;
mExtent = projector.mExtent;
mCoordinateTransform.setSourceCrs( mSrcCRS );
mCoordinateTransform.setDestCRS( mDestCRS );
}
return *this;
}
Expand All @@ -97,6 +125,8 @@ QgsRasterInterface * QgsRasterProjector::clone() const
{
QgsDebugMsg( "Entered" );
QgsRasterProjector * projector = new QgsRasterProjector( mSrcCRS, mDestCRS, mMaxSrcXRes, mMaxSrcYRes, mExtent );
projector->mSrcDatumTransform = mSrcDatumTransform;
projector->mDestDatumTransform = mDestDatumTransform;
return projector;
}

Expand All @@ -120,12 +150,12 @@ QGis::DataType QgsRasterProjector::dataType( int bandNo ) const
return QGis::UnknownDataType;
}

void QgsRasterProjector::setCRS( const QgsCoordinateReferenceSystem & theSrcCRS, const QgsCoordinateReferenceSystem & theDestCRS )
void QgsRasterProjector::setCRS( const QgsCoordinateReferenceSystem & theSrcCRS, const QgsCoordinateReferenceSystem & theDestCRS, int srcDatumTransform, int destDatumTransform )
{
mSrcCRS = theSrcCRS;
mDestCRS = theDestCRS;
mCoordinateTransform.setSourceCrs( theDestCRS );
mCoordinateTransform.setDestCRS( theSrcCRS );
mSrcDatumTransform = srcDatumTransform;
mDestDatumTransform = destDatumTransform;
}

void QgsRasterProjector::calc()
Expand Down Expand Up @@ -166,6 +196,8 @@ void QgsRasterProjector::calc()
double myDestRes = mDestXRes < mDestYRes ? mDestXRes : mDestYRes;
mSqrTolerance = myDestRes * myDestRes;

const QgsCoordinateTransform* ct = QgsCoordinateTransformCache::instance()->transform( mDestCRS.authid(), mSrcCRS.authid(), mDestDatumTransform, mSrcDatumTransform );

// Initialize the matrix by corners and middle points
mCPCols = mCPRows = 3;
for ( int i = 0; i < mCPRows; i++ )
Expand All @@ -184,20 +216,20 @@ void QgsRasterProjector::calc()
}
for ( int i = 0; i < mCPRows; i++ )
{
calcRow( i );
calcRow( i, ct );
}

while ( true )
{
bool myColsOK = checkCols();
bool myColsOK = checkCols( ct );
if ( !myColsOK )
{
insertRows();
insertRows( ct );
}
bool myRowsOK = checkRows();
bool myRowsOK = checkRows( ct );
if ( !myRowsOK )
{
insertCols();
insertCols( ct );
}
if ( myColsOK && myRowsOK )
{
Expand Down Expand Up @@ -432,19 +464,19 @@ void QgsRasterProjector::nextHelper()
mHelperTopRow++;
}

void QgsRasterProjector::srcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol )
void QgsRasterProjector::srcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform* ct )
{
if ( mApproximate )
{
approximateSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol );
}
else
{
preciseSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol );
preciseSrcRowCol( theDestRow, theDestCol, theSrcRow, theSrcCol, ct );
}
}

void QgsRasterProjector::preciseSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol )
void QgsRasterProjector::preciseSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform* ct )
{
#ifdef QGISDEBUG
QgsDebugMsgLevel( QString( "theDestRow = %1" ).arg( theDestRow ), 5 );
Expand All @@ -460,7 +492,10 @@ void QgsRasterProjector::preciseSrcRowCol( int theDestRow, int theDestCol, int *
QgsDebugMsgLevel( QString( "x = %1 y = %2" ).arg( x ).arg( y ), 5 );
#endif

mCoordinateTransform.transformInPlace( x, y, z );
if ( ct )
{
ct->transformInPlace( x, y, z );
}

#ifdef QGISDEBUG
QgsDebugMsgLevel( QString( "x = %1 y = %2" ).arg( x ).arg( y ), 5 );
Expand Down Expand Up @@ -545,7 +580,7 @@ void QgsRasterProjector::approximateSrcRowCol( int theDestRow, int theDestCol, i
Q_ASSERT( *theSrcCol < mSrcCols );
}

void QgsRasterProjector::insertRows()
void QgsRasterProjector::insertRows( const QgsCoordinateTransform* ct )
{
for ( int r = 0; r < mCPRows - 1; r++ )
{
Expand All @@ -563,11 +598,11 @@ void QgsRasterProjector::insertRows()
mCPRows += mCPRows - 1;
for ( int r = 1; r < mCPRows - 1; r += 2 )
{
calcRow( r );
calcRow( r, ct );
}
}

void QgsRasterProjector::insertCols()
void QgsRasterProjector::insertCols( const QgsCoordinateTransform* ct )
{
for ( int r = 0; r < mCPRows; r++ )
{
Expand All @@ -582,20 +617,27 @@ void QgsRasterProjector::insertCols()
mCPCols += mCPCols - 1;
for ( int c = 1; c < mCPCols - 1; c += 2 )
{
calcCol( c );
calcCol( c, ct );
}

}

void QgsRasterProjector::calcCP( int theRow, int theCol )
void QgsRasterProjector::calcCP( int theRow, int theCol, const QgsCoordinateTransform* ct )
{
double myDestX, myDestY;
destPointOnCPMatrix( theRow, theCol, &myDestX, &myDestY );
QgsPoint myDestPoint( myDestX, myDestY );
try
{
mCPMatrix[theRow][theCol] = mCoordinateTransform.transform( myDestPoint );
mCPLegalMatrix[theRow][theCol] = true;
if ( ct )
{
mCPMatrix[theRow][theCol] = ct->transform( myDestPoint );
mCPLegalMatrix[theRow][theCol] = true;
}
else
{
mCPLegalMatrix[theRow][theCol] = false;
}
}
catch ( QgsCsException &e )
{
Expand All @@ -605,30 +647,35 @@ void QgsRasterProjector::calcCP( int theRow, int theCol )
}
}

bool QgsRasterProjector::calcRow( int theRow )
bool QgsRasterProjector::calcRow( int theRow, const QgsCoordinateTransform* ct )
{
QgsDebugMsgLevel( QString( "theRow = %1" ).arg( theRow ), 3 );
for ( int i = 0; i < mCPCols; i++ )
{
calcCP( theRow, i );
calcCP( theRow, i, ct );
}

return true;
}

bool QgsRasterProjector::calcCol( int theCol )
bool QgsRasterProjector::calcCol( int theCol, const QgsCoordinateTransform* ct )
{
QgsDebugMsgLevel( QString( "theCol = %1" ).arg( theCol ), 3 );
for ( int i = 0; i < mCPRows; i++ )
{
calcCP( i, theCol );
calcCP( i, theCol, ct );
}

return true;
}

bool QgsRasterProjector::checkCols()
bool QgsRasterProjector::checkCols( const QgsCoordinateTransform* ct )
{
if ( !ct )
{
return false;
}

for ( int c = 0; c < mCPCols; c++ )
{
for ( int r = 1; r < mCPRows - 1; r += 2 )
Expand All @@ -649,7 +696,7 @@ bool QgsRasterProjector::checkCols()
}
try
{
QgsPoint myDestApprox = mCoordinateTransform.transform( mySrcApprox, QgsCoordinateTransform::ReverseTransform );
QgsPoint myDestApprox = ct->transform( mySrcApprox, QgsCoordinateTransform::ReverseTransform );
double mySqrDist = myDestApprox.sqrDist( myDestPoint );
if ( mySqrDist > mSqrTolerance )
{
Expand All @@ -667,8 +714,13 @@ bool QgsRasterProjector::checkCols()
return true;
}

bool QgsRasterProjector::checkRows()
bool QgsRasterProjector::checkRows( const QgsCoordinateTransform* ct )
{
if ( !ct )
{
return false;
}

for ( int r = 0; r < mCPRows; r++ )
{
for ( int c = 1; c < mCPCols - 1; c += 2 )
Expand All @@ -689,7 +741,7 @@ bool QgsRasterProjector::checkRows()
}
try
{
QgsPoint myDestApprox = mCoordinateTransform.transform( mySrcApprox, QgsCoordinateTransform::ReverseTransform );
QgsPoint myDestApprox = ct->transform( mySrcApprox, QgsCoordinateTransform::ReverseTransform );
double mySqrDist = myDestApprox.sqrDist( myDestPoint );
if ( mySqrDist > mSqrTolerance )
{
Expand Down Expand Up @@ -775,12 +827,18 @@ QgsRasterBlock * QgsRasterProjector::block( int bandNo, QgsRectangle const & ex
// we cannot fill output block with no data because we use memcpy for data, not setValue().
bool doNoData = inputBlock->hasNoData() && !inputBlock->hasNoDataValue();

const QgsCoordinateTransform* ct = 0;
if ( !mApproximate )
{
ct = QgsCoordinateTransformCache::instance()->transform( mDestCRS.authid(), mSrcCRS.authid(), mDestDatumTransform, mSrcDatumTransform );
}

int srcRow, srcCol;
for ( int i = 0; i < height; ++i )
{
for ( int j = 0; j < width; ++j )
{
srcRowCol( i, j, &srcRow, &srcCol );
srcRowCol( i, j, &srcRow, &srcCol, ct );
qgssize srcIndex = ( qgssize )srcRow * mSrcCols + srcCol;
QgsDebugMsgLevel( QString( "row = %1 col = %2 srcRow = %3 srcCol = %4" ).arg( i ).arg( j ).arg( srcRow ).arg( srcCol ), 5 );

Expand Down
40 changes: 28 additions & 12 deletions src/core/raster/qgsrasterprojector.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,18 @@ class CORE_EXPORT QgsRasterProjector : public QgsRasterInterface
* it calculates grid of points in source CRS for target CRS + extent
* which are used to calculate affine transformation matrices.
*/

QgsRasterProjector(
QgsCoordinateReferenceSystem theSrcCRS,
QgsCoordinateReferenceSystem theDestCRS,
int theSrcDatumTransform,
int theDestDatumTransform,
QgsRectangle theDestExtent,
int theDestRows, int theDestCols,
double theMaxSrcXRes, double theMaxSrcYRes,
QgsRectangle theExtent
);

QgsRasterProjector(
QgsCoordinateReferenceSystem theSrcCRS,
QgsCoordinateReferenceSystem theDestCRS,
Expand Down Expand Up @@ -74,7 +86,8 @@ class CORE_EXPORT QgsRasterProjector : public QgsRasterInterface
QGis::DataType dataType( int bandNo ) const;

/** \brief set source and destination CRS */
void setCRS( const QgsCoordinateReferenceSystem & theSrcCRS, const QgsCoordinateReferenceSystem & theDestCRS );
void setCRS( const QgsCoordinateReferenceSystem & theSrcCRS, const QgsCoordinateReferenceSystem & theDestCRS,
int srcDatumTransform = -1, int destDatumTransform = -1 );

/** \brief Get source CRS */
QgsCoordinateReferenceSystem srcCrs() const { return mSrcCRS; }
Expand All @@ -101,7 +114,7 @@ class CORE_EXPORT QgsRasterProjector : public QgsRasterInterface
void setSrcCols( int theCols ) { mSrcCols = theCols; mSrcYRes = mSrcExtent.width() / mSrcCols; }

/** \brief Get source row and column indexes for current source extent and resolution */
void srcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol );
void srcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform* ct );

int dstRows() const { return mDestRows; }
int dstCols() const { return mDestCols; }
Expand All @@ -117,7 +130,7 @@ class CORE_EXPORT QgsRasterProjector : public QgsRasterInterface
QgsPoint srcPoint( int theRow, int theCol );

/** \brief Get precise source row and column indexes for current source extent and resolution */
inline void preciseSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol );
inline void preciseSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol, const QgsCoordinateTransform* ct );

/** \brief Get approximate source row and column indexes for current source extent and resolution */
inline void approximateSrcRowCol( int theDestRow, int theDestCol, int *theSrcRow, int *theSrcCol );
Expand All @@ -126,19 +139,19 @@ class CORE_EXPORT QgsRasterProjector : public QgsRasterInterface
void calc();

/** \brief insert rows to matrix */
void insertRows();
void insertRows( const QgsCoordinateTransform* ct );

/** \brief insert columns to matrix */
void insertCols();
void insertCols( const QgsCoordinateTransform* ct );

/* calculate single control point in current matrix */
void calcCP( int theRow, int theCol );
void calcCP( int theRow, int theCol, const QgsCoordinateTransform* ct );

/** \brief calculate matrix row */
bool calcRow( int theRow );
bool calcRow( int theRow, const QgsCoordinateTransform* ct );

/** \brief calculate matrix column */
bool calcCol( int theCol );
bool calcCol( int theCol, const QgsCoordinateTransform* ct );

/** \brief calculate source extent */
void calcSrcExtent();
Expand All @@ -148,11 +161,11 @@ class CORE_EXPORT QgsRasterProjector : public QgsRasterInterface

/** \brief check error along columns
* returns true if within threshold */
bool checkCols();
bool checkCols( const QgsCoordinateTransform* ct );

/** \brief check error along rows
* returns true if within threshold */
bool checkRows();
bool checkRows( const QgsCoordinateTransform* ct );

/** Calculate array of src helper points */
void calcHelper( int theMatrixRow, QgsPoint *thePoints );
Expand All @@ -169,8 +182,11 @@ class CORE_EXPORT QgsRasterProjector : public QgsRasterInterface
/** Destination CRS */
QgsCoordinateReferenceSystem mDestCRS;

/** Reverse coordinate transform (from destination to source) */
QgsCoordinateTransform mCoordinateTransform;
/** Source datum transformation id (or -1 if none) */
int mSrcDatumTransform;

/** Destination datum transformation id (or -1 if none) */
int mDestDatumTransform;

/** Destination extent */
QgsRectangle mDestExtent;
Expand Down
3 changes: 3 additions & 0 deletions src/core/raster/qgsrasterviewport.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,9 @@ struct QgsRasterViewPort

/** \brief Target coordinate system */
QgsCoordinateReferenceSystem mDestCRS;

int mSrcDatumTransform;
int mDestDatumTransform;
};

#endif //QGSRASTERVIEWPORT_H