Skip to content

Commit b452bc0

Browse files
wonder-skwonder-sk
committedJun 29, 2015
Support reprojection in GUI, added unit tests for reprojection
·
ltr-3_40archive/master_2
1 parent d7a9493 commit b452bc0

File tree

6 files changed

+166
-45
lines changed

6 files changed

+166
-45
lines changed
 

‎python/analysis/raster/qgsalignraster.sip

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -124,9 +124,9 @@ class QgsAlignRaster
124124
QgsRectangle clipExtent() const;
125125

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

131131
//! Run the alignment process
132132
//! @return true on success

‎src/analysis/raster/qgsalignraster.cpp

Lines changed: 77 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -145,19 +145,36 @@ QgsRectangle QgsAlignRaster::clipExtent() const
145145
}
146146

147147

148-
void QgsAlignRaster::setParametersFromRaster( const QString& filename )
148+
bool QgsAlignRaster::setParametersFromRaster( const QString& filename, const QString& destWkt )
149149
{
150-
setParametersFromRaster( RasterInfo( filename ) );
150+
return setParametersFromRaster( RasterInfo( filename ), destWkt );
151151
}
152152

153-
void QgsAlignRaster::setParametersFromRaster( const RasterInfo& rasterInfo )
153+
bool QgsAlignRaster::setParametersFromRaster( const RasterInfo& rasterInfo, const QString& destWkt )
154154
{
155-
// use ref. layer to init input
156-
mCrsWkt = rasterInfo.crs();
157-
mCellSizeX = rasterInfo.cellSize().width();
158-
mCellSizeY = rasterInfo.cellSize().height();
159-
mGridOffsetX = rasterInfo.gridOffset().x();
160-
mGridOffsetY = rasterInfo.gridOffset().y();
155+
if ( destWkt.isEmpty() || destWkt.toAscii() == rasterInfo.crs() )
156+
{
157+
// use ref. layer to init input
158+
mCrsWkt = rasterInfo.crs();
159+
mCellSizeX = rasterInfo.cellSize().width();
160+
mCellSizeY = rasterInfo.cellSize().height();
161+
mGridOffsetX = rasterInfo.gridOffset().x();
162+
mGridOffsetY = rasterInfo.gridOffset().y();
163+
}
164+
else
165+
{
166+
QSizeF cs;
167+
QPointF go;
168+
if ( !suggestedWarpOutput( rasterInfo, destWkt.toAscii(), &cs, &go ) )
169+
return false;
170+
171+
mCrsWkt = destWkt.toAscii();
172+
mCellSizeX = cs.width();
173+
mCellSizeY = cs.height();
174+
mGridOffsetX = go.x();
175+
mGridOffsetY = go.y();
176+
}
177+
return true;
161178
}
162179

163180

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

173190
RasterInfo info( r.inputFilename );
174191

175-
// Create a transformer that maps from source pixel/line coordinates
176-
// to destination georeferenced coordinates (not destination
177-
// pixel line). We do that by omitting the destination dataset
178-
// handle (setting it to NULL).
179-
void* hTransformArg = GDALCreateGenImgProjTransformer( info.mDataset, info.mCrsWkt.constData(), NULL, mCrsWkt.constData(), FALSE, 0, 1 );
180-
if ( !hTransformArg )
192+
QSizeF cs;
193+
QgsRectangle extent;
194+
if ( !suggestedWarpOutput( info, mCrsWkt, &cs, 0, &extent ) )
181195
{
182-
mErrorMessage = QString( "GDALCreateGenImgProjTransformer failed.\n\n"
183-
"Source WKT:\n%1\n\nDestination WKT:\n%2" )
196+
mErrorMessage = QString( "Failed to get suggested warp output.\n\n"
197+
"File:\n%1\n\n"
198+
"Source WKT:\n%2\n\nDestination WKT:\n%3" )
199+
.arg( r.inputFilename )
184200
.arg( QString::fromAscii( info.mCrsWkt ) )
185201
.arg( QString::fromAscii( mCrsWkt ) );
186202
return false;
187203
}
188204

189-
// Get approximate output georeferenced bounds and resolution for file.
190-
double adfDstGeoTransform[6];
191-
double extents[4];
192-
int nPixels = 0, nLines = 0;
193-
CPLErr eErr;
194-
eErr = GDALSuggestedWarpOutput2( info.mDataset,
195-
GDALGenImgProjTransform, hTransformArg,
196-
adfDstGeoTransform, &nPixels, &nLines, extents, 0 );
197-
GDALDestroyGenImgProjTransformer( hTransformArg );
198-
r.srcCellSizeInDestCRS = fabs( adfDstGeoTransform[1] * adfDstGeoTransform[5] );
199-
200-
if ( eErr != CE_None )
201-
{
202-
mErrorMessage = QString( "GDALSuggestedWarpOutput2 failed.\n\n" + r.inputFilename );
203-
return false;
204-
}
205+
r.srcCellSizeInDestCRS = cs.width() * cs.height();
205206

206207
if ( finalExtent[0] == 0 && finalExtent[1] == 0 && finalExtent[2] == 0 && finalExtent[3] == 0 )
207208
{
208209
// initialize with the first layer
209-
finalExtent[0] = extents[0];
210-
finalExtent[1] = extents[1];
211-
finalExtent[2] = extents[2];
212-
finalExtent[3] = extents[3];
210+
finalExtent[0] = extent.xMinimum();
211+
finalExtent[1] = extent.yMinimum();
212+
finalExtent[2] = extent.xMaximum();
213+
finalExtent[3] = extent.yMaximum();
213214
}
214215
else
215216
{
216217
// use intersection of rects
217-
if ( extents[0] > finalExtent[0] ) finalExtent[0] = extents[0];
218-
if ( extents[1] > finalExtent[1] ) finalExtent[1] = extents[1];
219-
if ( extents[2] < finalExtent[2] ) finalExtent[2] = extents[2];
220-
if ( extents[3] < finalExtent[3] ) finalExtent[3] = extents[3];
218+
if ( extent.xMinimum() > finalExtent[0] ) finalExtent[0] = extent.xMinimum();
219+
if ( extent.yMinimum() > finalExtent[1] ) finalExtent[1] = extent.yMinimum();
220+
if ( extent.xMaximum() < finalExtent[2] ) finalExtent[2] = extent.xMaximum();
221+
if ( extent.yMaximum() < finalExtent[3] ) finalExtent[3] = extent.yMaximum();
221222
}
222223
}
223224

@@ -393,6 +394,41 @@ bool QgsAlignRaster::createAndWarp( const Item& raster )
393394
return true;
394395
}
395396

397+
bool QgsAlignRaster::suggestedWarpOutput( const QgsAlignRaster::RasterInfo& info, const QByteArray& destWkt, QSizeF* cellSize, QPointF* gridOffset, QgsRectangle* rect )
398+
{
399+
// Create a transformer that maps from source pixel/line coordinates
400+
// to destination georeferenced coordinates (not destination
401+
// pixel line). We do that by omitting the destination dataset
402+
// handle (setting it to NULL).
403+
void* hTransformArg = GDALCreateGenImgProjTransformer( info.mDataset, info.mCrsWkt.constData(), NULL, destWkt.constData(), FALSE, 0, 1 );
404+
if ( !hTransformArg )
405+
return false;
406+
407+
// Get approximate output georeferenced bounds and resolution for file.
408+
double adfDstGeoTransform[6];
409+
double extents[4];
410+
int nPixels = 0, nLines = 0;
411+
CPLErr eErr;
412+
eErr = GDALSuggestedWarpOutput2( info.mDataset,
413+
GDALGenImgProjTransform, hTransformArg,
414+
adfDstGeoTransform, &nPixels, &nLines, extents, 0 );
415+
GDALDestroyGenImgProjTransformer( hTransformArg );
416+
417+
if ( eErr != CE_None )
418+
return false;
419+
420+
QSizeF cs( qAbs( adfDstGeoTransform[1] ), qAbs( adfDstGeoTransform[5] ) );
421+
422+
if ( rect )
423+
*rect = QgsRectangle( extents[0], extents[1], extents[2], extents[3] );
424+
if ( cellSize )
425+
*cellSize = cs;
426+
if ( gridOffset )
427+
*gridOffset = QPointF( fmod_with_tolerance( adfDstGeoTransform[0], cs.width() ),
428+
fmod_with_tolerance( adfDstGeoTransform[3], cs.height() ) );
429+
return true;
430+
}
431+
396432

397433
//----------
398434

‎src/analysis/raster/qgsalignraster.h

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,8 @@ typedef void* GDALDatasetH;
3232
* - coordinate reference system
3333
* - cell size and raster size
3434
* - offset of the raster grid
35+
*
36+
* @note added in QGIS 2.12
3537
*/
3638
class ANALYSIS_EXPORT QgsAlignRaster
3739
{
@@ -170,9 +172,9 @@ class ANALYSIS_EXPORT QgsAlignRaster
170172
QgsRectangle clipExtent() const;
171173

172174
//! Set destination CRS, cell size and grid offset from a raster file
173-
void setParametersFromRaster( const RasterInfo& rasterInfo );
175+
bool setParametersFromRaster( const RasterInfo& rasterInfo, const QString& destWkt = QString() );
174176
//! Set destination CRS, cell size and grid offset from a raster file
175-
void setParametersFromRaster( const QString& filename );
177+
bool setParametersFromRaster( const QString& filename, const QString& destWkt = QString() );
176178

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

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

198203
// set by the client

‎src/app/qgsalignrasterdialog.cpp

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -72,6 +72,7 @@ QgsAlignRasterDialog::QgsAlignRasterDialog( QWidget *parent )
7272
connect( mBtnEdit, SIGNAL( clicked( bool ) ), this, SLOT( editLayer() ) );
7373

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

7677
mClipExtentGroupBox->setChecked( false );
7778

@@ -184,6 +185,31 @@ void QgsAlignRasterDialog::updateConfigFromReferenceLayer()
184185
mClipExtentGroupBox->setOriginalExtent( mAlign->clipExtent(), destCRS );
185186
}
186187

188+
void QgsAlignRasterDialog::destinationCrsChanged()
189+
{
190+
if ( mCrsSelector->crs().toWkt() == mAlign->destinationCRS() )
191+
return;
192+
193+
int index = mCboReferenceLayer->currentIndex();
194+
if ( index < 0 )
195+
return;
196+
197+
if ( !mAlign->setParametersFromRaster( mAlign->rasters().at( index ).inputFilename, mCrsSelector->crs().toWkt() ) )
198+
{
199+
QMessageBox::warning( this, tr( "Align Rasters" ), tr( "Cannot reproject reference layer to the chosen destination CRS.\n\nPlease select a different CRS" ) );
200+
return;
201+
}
202+
203+
QSizeF cellSize = mAlign->cellSize();
204+
mSpinCellSizeX->setValue( cellSize.width() );
205+
mSpinCellSizeY->setValue( cellSize.height() );
206+
207+
QPointF gridOffset = mAlign->gridOffset();
208+
mSpinGridOffsetX->setValue( gridOffset.x() );
209+
mSpinGridOffsetY->setValue( gridOffset.y() );
210+
}
211+
212+
187213
void QgsAlignRasterDialog::runAlign()
188214
{
189215
setEnabled( false );

‎src/app/qgsalignrasterdialog.h

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,8 @@ class QgsAlignRasterDialog : public QDialog, private Ui::QgsAlignRasterDialog
2626

2727
void runAlign();
2828

29+
void destinationCrsChanged();
30+
2931
protected:
3032
void populateLayersView();
3133

‎tests/src/analysis/testqgsalignraster.cpp

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,8 @@
1616
#include <QtTest/QtTest>
1717

1818
#include "qgsalignraster.h"
19+
#include "qgsapplication.h"
20+
#include "qgscoordinatereferencesystem.h"
1921
#include "qgsrectangle.h"
2022

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

4446
SRC_FILE = QString( TEST_DATA_DIR ) + QDir::separator() + "float1-16.tif";
47+
48+
QgsApplication::init(); // needed for CRS database
4549
}
4650

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

217+
void testReprojectToOtherCRS()
218+
{
219+
QString tmpFile( _tempFile( "reproject-utm-47n" ) );
220+
221+
// reproject from WGS84 to UTM zone 47N
222+
// (the true UTM zone for this raster is 48N, but here it is
223+
// more obvious the different shape of the resulting raster)
224+
QgsCoordinateReferenceSystem destCRS( "EPSG:32647" );
225+
QVERIFY( destCRS.isValid() );
226+
227+
QgsAlignRaster align;
228+
QgsAlignRaster::List rasters;
229+
rasters << QgsAlignRaster::Item( SRC_FILE, tmpFile );
230+
align.setRasters( rasters );
231+
align.setParametersFromRaster( SRC_FILE, destCRS.toWkt() );
232+
bool res = align.run();
233+
QVERIFY( res );
234+
235+
QgsAlignRaster::RasterInfo out( tmpFile );
236+
QVERIFY( out.isValid() );
237+
QgsCoordinateReferenceSystem outCRS( QString::fromAscii( out.crs() ) );
238+
QCOMPARE( outCRS, destCRS );
239+
QCOMPARE( out.rasterSize(), QSize( 4, 4 ) );
240+
// let's stick to integers to keep the test more robust
241+
QCOMPARE( qRound( out.cellSize().width() ), 22293 ); // ~ 22293.256065
242+
QCOMPARE( qRound( out.cellSize().height() ), 22293 ); // ~ 22293.256065
243+
QCOMPARE( qRound( out.gridOffset().x() ), 4327 ); // ~ 4327.168434
244+
QCOMPARE( qRound( out.gridOffset().y() ), 637 ); // ~ 637.007990
245+
QCOMPARE( out.identify( 1308405, -746611 ), 10. );
246+
}
247+
248+
void testInvalidReprojection()
249+
{
250+
QString tmpFile( _tempFile( "reproject-invalid" ) );
251+
252+
// reprojection to British National Grid with raster in Jakarta area clearly cannot work
253+
QgsCoordinateReferenceSystem destCRS( "EPSG:27700" );
254+
QVERIFY( destCRS.isValid() );
255+
256+
QgsAlignRaster align;
257+
QgsAlignRaster::List rasters;
258+
rasters << QgsAlignRaster::Item( SRC_FILE, tmpFile );
259+
align.setRasters( rasters );
260+
align.setParametersFromRaster( SRC_FILE, destCRS.toWkt() );
261+
bool res = align.run();
262+
QVERIFY( !res );
263+
}
264+
213265
};
214266

215267
QTEST_MAIN( TestAlignRaster )

0 commit comments

Comments
 (0)
Please sign in to comment.