Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Support reprojection in GUI, added unit tests for reprojection
  • Loading branch information
wonder-sk committed Jun 29, 2015
1 parent d7a9493 commit b452bc0
Show file tree
Hide file tree
Showing 6 changed files with 166 additions and 45 deletions.
4 changes: 2 additions & 2 deletions python/analysis/raster/qgsalignraster.sip
Expand Up @@ -124,9 +124,9 @@ class QgsAlignRaster
QgsRectangle clipExtent() const;

//! Set destination CRS, cell size and grid offset from a raster file
void setParametersFromRaster( const RasterInfo& rasterInfo );
bool setParametersFromRaster( const RasterInfo& rasterInfo, const QString& destWkt = QString() );
//! Set destination CRS, cell size and grid offset from a raster file
void setParametersFromRaster( const QString& filename );
bool setParametersFromRaster( const QString& filename, const QString& destWkt = QString() );

//! Run the alignment process
//! @return true on success
Expand Down
118 changes: 77 additions & 41 deletions src/analysis/raster/qgsalignraster.cpp
Expand Up @@ -145,19 +145,36 @@ QgsRectangle QgsAlignRaster::clipExtent() const
}


void QgsAlignRaster::setParametersFromRaster( const QString& filename )
bool QgsAlignRaster::setParametersFromRaster( const QString& filename, const QString& destWkt )
{
setParametersFromRaster( RasterInfo( filename ) );
return setParametersFromRaster( RasterInfo( filename ), destWkt );
}

void QgsAlignRaster::setParametersFromRaster( const RasterInfo& rasterInfo )
bool QgsAlignRaster::setParametersFromRaster( const RasterInfo& rasterInfo, const QString& destWkt )
{
// use ref. layer to init input
mCrsWkt = rasterInfo.crs();
mCellSizeX = rasterInfo.cellSize().width();
mCellSizeY = rasterInfo.cellSize().height();
mGridOffsetX = rasterInfo.gridOffset().x();
mGridOffsetY = rasterInfo.gridOffset().y();
if ( destWkt.isEmpty() || destWkt.toAscii() == rasterInfo.crs() )
{
// use ref. layer to init input
mCrsWkt = rasterInfo.crs();
mCellSizeX = rasterInfo.cellSize().width();
mCellSizeY = rasterInfo.cellSize().height();
mGridOffsetX = rasterInfo.gridOffset().x();
mGridOffsetY = rasterInfo.gridOffset().y();
}
else
{
QSizeF cs;
QPointF go;
if ( !suggestedWarpOutput( rasterInfo, destWkt.toAscii(), &cs, &go ) )
return false;

mCrsWkt = destWkt.toAscii();
mCellSizeX = cs.width();
mCellSizeY = cs.height();
mGridOffsetX = go.x();
mGridOffsetY = go.y();
}
return true;
}


Expand All @@ -172,52 +189,36 @@ bool QgsAlignRaster::determineTransformAndSize()

RasterInfo info( r.inputFilename );

// Create a transformer that maps from source pixel/line coordinates
// to destination georeferenced coordinates (not destination
// pixel line). We do that by omitting the destination dataset
// handle (setting it to NULL).
void* hTransformArg = GDALCreateGenImgProjTransformer( info.mDataset, info.mCrsWkt.constData(), NULL, mCrsWkt.constData(), FALSE, 0, 1 );
if ( !hTransformArg )
QSizeF cs;
QgsRectangle extent;
if ( !suggestedWarpOutput( info, mCrsWkt, &cs, 0, &extent ) )
{
mErrorMessage = QString( "GDALCreateGenImgProjTransformer failed.\n\n"
"Source WKT:\n%1\n\nDestination WKT:\n%2" )
mErrorMessage = QString( "Failed to get suggested warp output.\n\n"
"File:\n%1\n\n"
"Source WKT:\n%2\n\nDestination WKT:\n%3" )
.arg( r.inputFilename )
.arg( QString::fromAscii( info.mCrsWkt ) )
.arg( QString::fromAscii( mCrsWkt ) );
return false;
}

// Get approximate output georeferenced bounds and resolution for file.
double adfDstGeoTransform[6];
double extents[4];
int nPixels = 0, nLines = 0;
CPLErr eErr;
eErr = GDALSuggestedWarpOutput2( info.mDataset,
GDALGenImgProjTransform, hTransformArg,
adfDstGeoTransform, &nPixels, &nLines, extents, 0 );
GDALDestroyGenImgProjTransformer( hTransformArg );
r.srcCellSizeInDestCRS = fabs( adfDstGeoTransform[1] * adfDstGeoTransform[5] );

if ( eErr != CE_None )
{
mErrorMessage = QString( "GDALSuggestedWarpOutput2 failed.\n\n" + r.inputFilename );
return false;
}
r.srcCellSizeInDestCRS = cs.width() * cs.height();

if ( finalExtent[0] == 0 && finalExtent[1] == 0 && finalExtent[2] == 0 && finalExtent[3] == 0 )
{
// initialize with the first layer
finalExtent[0] = extents[0];
finalExtent[1] = extents[1];
finalExtent[2] = extents[2];
finalExtent[3] = extents[3];
finalExtent[0] = extent.xMinimum();
finalExtent[1] = extent.yMinimum();
finalExtent[2] = extent.xMaximum();
finalExtent[3] = extent.yMaximum();
}
else
{
// use intersection of rects
if ( extents[0] > finalExtent[0] ) finalExtent[0] = extents[0];
if ( extents[1] > finalExtent[1] ) finalExtent[1] = extents[1];
if ( extents[2] < finalExtent[2] ) finalExtent[2] = extents[2];
if ( extents[3] < finalExtent[3] ) finalExtent[3] = extents[3];
if ( extent.xMinimum() > finalExtent[0] ) finalExtent[0] = extent.xMinimum();
if ( extent.yMinimum() > finalExtent[1] ) finalExtent[1] = extent.yMinimum();
if ( extent.xMaximum() < finalExtent[2] ) finalExtent[2] = extent.xMaximum();
if ( extent.yMaximum() < finalExtent[3] ) finalExtent[3] = extent.yMaximum();
}
}

Expand Down Expand Up @@ -393,6 +394,41 @@ bool QgsAlignRaster::createAndWarp( const Item& raster )
return true;
}

bool QgsAlignRaster::suggestedWarpOutput( const QgsAlignRaster::RasterInfo& info, const QByteArray& destWkt, QSizeF* cellSize, QPointF* gridOffset, QgsRectangle* rect )
{
// Create a transformer that maps from source pixel/line coordinates
// to destination georeferenced coordinates (not destination
// pixel line). We do that by omitting the destination dataset
// handle (setting it to NULL).
void* hTransformArg = GDALCreateGenImgProjTransformer( info.mDataset, info.mCrsWkt.constData(), NULL, destWkt.constData(), FALSE, 0, 1 );
if ( !hTransformArg )
return false;

// Get approximate output georeferenced bounds and resolution for file.
double adfDstGeoTransform[6];
double extents[4];
int nPixels = 0, nLines = 0;
CPLErr eErr;
eErr = GDALSuggestedWarpOutput2( info.mDataset,
GDALGenImgProjTransform, hTransformArg,
adfDstGeoTransform, &nPixels, &nLines, extents, 0 );
GDALDestroyGenImgProjTransformer( hTransformArg );

if ( eErr != CE_None )
return false;

QSizeF cs( qAbs( adfDstGeoTransform[1] ), qAbs( adfDstGeoTransform[5] ) );

if ( rect )
*rect = QgsRectangle( extents[0], extents[1], extents[2], extents[3] );
if ( cellSize )
*cellSize = cs;
if ( gridOffset )
*gridOffset = QPointF( fmod_with_tolerance( adfDstGeoTransform[0], cs.width() ),
fmod_with_tolerance( adfDstGeoTransform[3], cs.height() ) );
return true;
}


//----------

Expand Down
9 changes: 7 additions & 2 deletions src/analysis/raster/qgsalignraster.h
Expand Up @@ -32,6 +32,8 @@ typedef void* GDALDatasetH;
* - coordinate reference system
* - cell size and raster size
* - offset of the raster grid
*
* @note added in QGIS 2.12
*/
class ANALYSIS_EXPORT QgsAlignRaster
{
Expand Down Expand Up @@ -170,9 +172,9 @@ class ANALYSIS_EXPORT QgsAlignRaster
QgsRectangle clipExtent() const;

//! Set destination CRS, cell size and grid offset from a raster file
void setParametersFromRaster( const RasterInfo& rasterInfo );
bool setParametersFromRaster( const RasterInfo& rasterInfo, const QString& destWkt = QString() );
//! Set destination CRS, cell size and grid offset from a raster file
void setParametersFromRaster( const QString& filename );
bool setParametersFromRaster( const QString& filename, const QString& destWkt = QString() );

//! Run the alignment process
//! @return true on success
Expand All @@ -193,6 +195,9 @@ class ANALYSIS_EXPORT QgsAlignRaster
//! Internal function for processing of one raster (1. create output, 2. do the alignment)
bool createAndWarp( const Item& raster );

//! Determine suggested output of raster warp to a different CRS. Returns true on success
static bool suggestedWarpOutput( const RasterInfo& info, const QByteArray& destWkt, QSizeF* cellSize = 0, QPointF* gridOffset = 0, QgsRectangle* rect = 0 );

protected:

// set by the client
Expand Down
26 changes: 26 additions & 0 deletions src/app/qgsalignrasterdialog.cpp
Expand Up @@ -72,6 +72,7 @@ QgsAlignRasterDialog::QgsAlignRasterDialog( QWidget *parent )
connect( mBtnEdit, SIGNAL( clicked( bool ) ), this, SLOT( editLayer() ) );

connect( mCboReferenceLayer, SIGNAL( currentIndexChanged( int ) ), this, SLOT( updateConfigFromReferenceLayer() ) );
connect( mCrsSelector, SIGNAL( crsChanged( QgsCoordinateReferenceSystem ) ), this, SLOT( destinationCrsChanged() ) );

mClipExtentGroupBox->setChecked( false );

Expand Down Expand Up @@ -184,6 +185,31 @@ void QgsAlignRasterDialog::updateConfigFromReferenceLayer()
mClipExtentGroupBox->setOriginalExtent( mAlign->clipExtent(), destCRS );
}

void QgsAlignRasterDialog::destinationCrsChanged()
{
if ( mCrsSelector->crs().toWkt() == mAlign->destinationCRS() )
return;

int index = mCboReferenceLayer->currentIndex();
if ( index < 0 )
return;

if ( !mAlign->setParametersFromRaster( mAlign->rasters().at( index ).inputFilename, mCrsSelector->crs().toWkt() ) )
{
QMessageBox::warning( this, tr( "Align Rasters" ), tr( "Cannot reproject reference layer to the chosen destination CRS.\n\nPlease select a different CRS" ) );
return;
}

QSizeF cellSize = mAlign->cellSize();
mSpinCellSizeX->setValue( cellSize.width() );
mSpinCellSizeY->setValue( cellSize.height() );

QPointF gridOffset = mAlign->gridOffset();
mSpinGridOffsetX->setValue( gridOffset.x() );
mSpinGridOffsetY->setValue( gridOffset.y() );
}


void QgsAlignRasterDialog::runAlign()
{
setEnabled( false );
Expand Down
2 changes: 2 additions & 0 deletions src/app/qgsalignrasterdialog.h
Expand Up @@ -26,6 +26,8 @@ class QgsAlignRasterDialog : public QDialog, private Ui::QgsAlignRasterDialog

void runAlign();

void destinationCrsChanged();

protected:
void populateLayersView();

Expand Down
52 changes: 52 additions & 0 deletions tests/src/analysis/testqgsalignraster.cpp
Expand Up @@ -16,6 +16,8 @@
#include <QtTest/QtTest>

#include "qgsalignraster.h"
#include "qgsapplication.h"
#include "qgscoordinatereferencesystem.h"
#include "qgsrectangle.h"

#include <QDir>
Expand All @@ -42,6 +44,8 @@ class TestAlignRaster : public QObject
GDALAllRegister();

SRC_FILE = QString( TEST_DATA_DIR ) + QDir::separator() + "float1-16.tif";

QgsApplication::init(); // needed for CRS database
}

void testRasterInfo()
Expand Down Expand Up @@ -210,6 +214,54 @@ class TestAlignRaster : public QObject
QCOMPARE( out.identify( 106.2, -6.4 ), 14. ); // = (1+2+5+6)
}

void testReprojectToOtherCRS()
{
QString tmpFile( _tempFile( "reproject-utm-47n" ) );

// reproject from WGS84 to UTM zone 47N
// (the true UTM zone for this raster is 48N, but here it is
// more obvious the different shape of the resulting raster)
QgsCoordinateReferenceSystem destCRS( "EPSG:32647" );
QVERIFY( destCRS.isValid() );

QgsAlignRaster align;
QgsAlignRaster::List rasters;
rasters << QgsAlignRaster::Item( SRC_FILE, tmpFile );
align.setRasters( rasters );
align.setParametersFromRaster( SRC_FILE, destCRS.toWkt() );
bool res = align.run();
QVERIFY( res );

QgsAlignRaster::RasterInfo out( tmpFile );
QVERIFY( out.isValid() );
QgsCoordinateReferenceSystem outCRS( QString::fromAscii( out.crs() ) );
QCOMPARE( outCRS, destCRS );
QCOMPARE( out.rasterSize(), QSize( 4, 4 ) );
// let's stick to integers to keep the test more robust
QCOMPARE( qRound( out.cellSize().width() ), 22293 ); // ~ 22293.256065
QCOMPARE( qRound( out.cellSize().height() ), 22293 ); // ~ 22293.256065
QCOMPARE( qRound( out.gridOffset().x() ), 4327 ); // ~ 4327.168434
QCOMPARE( qRound( out.gridOffset().y() ), 637 ); // ~ 637.007990
QCOMPARE( out.identify( 1308405, -746611 ), 10. );
}

void testInvalidReprojection()
{
QString tmpFile( _tempFile( "reproject-invalid" ) );

// reprojection to British National Grid with raster in Jakarta area clearly cannot work
QgsCoordinateReferenceSystem destCRS( "EPSG:27700" );
QVERIFY( destCRS.isValid() );

QgsAlignRaster align;
QgsAlignRaster::List rasters;
rasters << QgsAlignRaster::Item( SRC_FILE, tmpFile );
align.setRasters( rasters );
align.setParametersFromRaster( SRC_FILE, destCRS.toWkt() );
bool res = align.run();
QVERIFY( !res );
}

};

QTEST_MAIN( TestAlignRaster )
Expand Down

0 comments on commit b452bc0

Please sign in to comment.