Skip to content

Commit

Permalink
When saving/loading GCP files, use source layer coordinates for source
Browse files Browse the repository at this point in the history 
points instead of always using source pixel row/columns

The "always use pixels" approach does not work well when a source
image already has some georeferencing information attached.

Also add tests for this, and fix loading of "enabled" status for
points when loading points from a file
  • Loading branch information
@nyalldawson
nyalldawson committed Feb 9, 2022
1 parent fa424a4 commit baf0534
Show file tree
Hide file tree
Showing 5 changed files with 197 additions and 71 deletions.
91 changes: 91 additions & 0 deletions src/app/georeferencer/qgsgcplist.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include "qgsproject.h"

#include "qgsgcplist.h"
#include <QDir>

void QgsGCPList::createGCPVectors( QVector<QgsPointXY> &sourcePoints, QVector<QgsPointXY> &destinationPoints, const QgsCoordinateReferenceSystem &targetCrs, const QgsCoordinateTransformContext &context ) const
{
Expand Down Expand Up @@ -72,3 +73,93 @@ QList<QgsGcpPoint> QgsGCPList::asPoints() const
}
return res;
}

bool QgsGCPList::saveGcps( const QString &filePath, const QgsCoordinateReferenceSystem &targetCrs, const QgsCoordinateTransformContext &context, QString &error ) const
{
error.clear();

QFile pointFile( filePath );
if ( pointFile.open( QIODevice::WriteOnly | QIODevice::Truncate ) )
{
QTextStream points( &pointFile );
points << QStringLiteral( "#CRS: %1" ).arg( targetCrs.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED ) ) << endl;
points << "mapX,mapY,sourceX,sourceY,enable,dX,dY,residual" << endl;
for ( QgsGeorefDataPoint *pt : *this )
{
const QgsPointXY transformedDestinationPoint = pt->transformedDestinationPoint( targetCrs, context );
points << QStringLiteral( "%1,%2,%3,%4,%5,%6,%7,%8" )
.arg( qgsDoubleToString( transformedDestinationPoint.x() ),
qgsDoubleToString( transformedDestinationPoint.y() ),
qgsDoubleToString( pt->sourcePoint().x() ),
qgsDoubleToString( pt->sourcePoint().y() ) )
.arg( pt->isEnabled() )
.arg( qgsDoubleToString( pt->residual().x() ),
qgsDoubleToString( pt->residual().y() ),
qgsDoubleToString( std::sqrt( pt->residual().x() * pt->residual().x() + pt->residual().y() * pt->residual().y() ) ) )
<< endl;
}
return true;
}
else
{
error = QObject::tr( "Could not write to GCP points file %1." ).arg( QDir::toNativeSeparators( filePath ) );
return false;
}
}

QList<QgsGcpPoint> QgsGCPList::loadGcps( const QString &filePath, const QgsCoordinateReferenceSystem &defaultDestinationCrs, QString &error )
{
error.clear();
QFile pointFile( filePath );
if ( !pointFile.open( QIODevice::ReadOnly ) )
{
error = QObject::tr( "Could not open GCP points file %1." ).arg( QDir::toNativeSeparators( filePath ) );
return {};
}

QTextStream points( &pointFile );
int lineNumber = 0;
QString line = points.readLine();
lineNumber++;

int i = 0;
QgsCoordinateReferenceSystem destinationCrs;
if ( line.contains( QLatin1String( "#CRS: " ) ) )
{
destinationCrs = QgsCoordinateReferenceSystem( line.remove( QStringLiteral( "#CRS: " ) ) );
line = points.readLine();
lineNumber++;
}
else
destinationCrs = defaultDestinationCrs;

QList<QgsGcpPoint> res;
while ( !points.atEnd() )
{
line = points.readLine();
lineNumber++;
QStringList ls;
if ( line.contains( ',' ) ) // in previous format "\t" is delimiter of points in new - ","
ls = line.split( ',' ); // points from new georeferencer
else
ls = line.split( '\t' ); // points from prev georeferencer

if ( ls.count() < 4 )
{
error = QObject::tr( "Malformed content at line %1" ).arg( lineNumber );
return {};
}

const QgsPointXY destinationPoint( ls.at( 0 ).toDouble(), ls.at( 1 ).toDouble() ); // map x,y
const QgsPointXY sourcePoint( ls.at( 2 ).toDouble(), ls.at( 3 ).toDouble() ); // source x,y
bool enable = true;
if ( ls.count() >= 5 )
{
enable = ls.at( 4 ).toInt();
}
res.append( QgsGcpPoint( sourcePoint, destinationPoint, destinationCrs, enable ) );

++i;
}
return res;
}
27 changes: 27 additions & 0 deletions src/app/georeferencer/qgsgcplist.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,33 @@ class APP_EXPORT QgsGCPList : public QList<QgsGeorefDataPoint * >
*/
QList< QgsGcpPoint > asPoints() const;

/**
* Saves the GCPs to a text file.
*
* \param filePath destination file path
* \param targetCrs target CRS for destination points
* \param context transform context
* \param error will be set to a descriptive error message if saving fails
*
* \returns TRUE on success
*/
bool saveGcps( const QString &filePath,
const QgsCoordinateReferenceSystem &targetCrs, const QgsCoordinateTransformContext &context,
QString &error ) const;

/**
* Loads GCPs from a text file.
*
* \param filePath source file path
* \param defaultDestinationCrs default CRS to use for destination points if no destination CRS information is present in text file.
* \param error will be set to a descriptive error message if loading fails
*
* \returns TRUE on success
*/
static QList< QgsGcpPoint > loadGcps( const QString &filePath,
const QgsCoordinateReferenceSystem &defaultDestinationCrs,
QString &error );

};

#endif
85 changes: 15 additions & 70 deletions src/app/georeferencer/qgsgeorefmainwindow.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -278,7 +278,8 @@ void QgsGeoreferencerMainWindow::openRaster( const QString &fileName )

// load previously added points
mGCPpointsFileName = mRasterFileName + ".points";
( void )loadGCPs();
QString error;
( void )loadGCPs( error );

if ( mLayer )
mCanvas->setExtent( mLayer->extent() );
Expand Down Expand Up @@ -663,9 +664,10 @@ void QgsGeoreferencerMainWindow::loadGCPsDialog()
if ( mGCPpointsFileName.isEmpty() )
return;

if ( !loadGCPs() )
QString error;
if ( !loadGCPs( error ) )
{
mMessageBar->pushMessage( tr( "Load GCP Points" ), tr( "Invalid GCP file. File could not be read." ), Qgis::MessageLevel::Critical );
mMessageBar->pushMessage( tr( "Load GCP Points" ), error, Qgis::MessageLevel::Critical );
}
else
{
Expand Down Expand Up @@ -1238,52 +1240,18 @@ void QgsGeoreferencerMainWindow::writeSettings()
}

// GCP points
bool QgsGeoreferencerMainWindow::loadGCPs( /*bool verbose*/ )
bool QgsGeoreferencerMainWindow::loadGCPs( QString &error )
{
QFile pointFile( mGCPpointsFileName );
if ( !pointFile.open( QIODevice::ReadOnly ) )
{
const QList< QgsGcpPoint > points = QgsGCPList::loadGcps( mGCPpointsFileName,
QgsProject::instance()->crs(),
error );
if ( !error.isEmpty() )
return false;
}

clearGCPData();

QTextStream points( &pointFile );
QString line = points.readLine();
int i = 0;
QgsCoordinateReferenceSystem destinationCrs;
if ( line.contains( QLatin1String( "#CRS: " ) ) )
{
destinationCrs = QgsCoordinateReferenceSystem( line.remove( QStringLiteral( "#CRS: " ) ) );
line = points.readLine();
}
else
destinationCrs = QgsProject::instance()->crs();

while ( !points.atEnd() )
for ( const QgsGcpPoint &point : points )
{
line = points.readLine();
QStringList ls;
if ( line.contains( ',' ) ) // in previous format "\t" is delimiter of points in new - ","
ls = line.split( ',' ); // points from new georeferencer
else
ls = line.split( '\t' ); // points from prev georeferencer

if ( ls.count() < 4 )
return false;

const QgsPointXY destinationCoordinate( ls.at( 0 ).toDouble(), ls.at( 1 ).toDouble() ); // map x,y
const QgsPointXY pixelCoords( ls.at( 2 ).toDouble(), ls.at( 3 ).toDouble() ); // pixel x,y
const QgsPointXY sourceLayerCoordinate = mGeorefTransform.toSourceCoordinate( pixelCoords );
if ( ls.count() == 5 )
{
bool enable = ls.at( 4 ).toInt();
addPoint( sourceLayerCoordinate, destinationCoordinate, destinationCrs, enable, false );
}
else
addPoint( sourceLayerCoordinate, destinationCoordinate, destinationCrs, true, false );

++i;
addPoint( point.sourcePoint(), point.destinationPoint(), point.destinationPointCrs(), point.isEnabled(), false );
}

mSavedPoints = mPoints.asPoints();
Expand All @@ -1301,38 +1269,15 @@ bool QgsGeoreferencerMainWindow::loadGCPs( /*bool verbose*/ )

void QgsGeoreferencerMainWindow::saveGCPs()
{
QFile pointFile( mGCPpointsFileName );
if ( pointFile.open( QIODevice::WriteOnly | QIODevice::Truncate ) )
QString error;
if ( mPoints.saveGcps( mGCPpointsFileName, mProjection, QgsProject::instance()->transformContext(), error ) )
{
QTextStream points( &pointFile );
points << QStringLiteral( "#CRS: %1" ).arg( mProjection.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED ) ) << endl;
points << "mapX,mapY,pixelX,pixelY,enable,dX,dY,residual" << endl;
for ( QgsGeorefDataPoint *pt : std::as_const( mPoints ) )
{
const QgsPointXY sourcePixel = mGeorefTransform.toSourcePixel( pt->sourcePoint() );

const QgsPointXY transformedDestinationPoint = pt->transformedDestinationPoint( mProjection, QgsProject::instance()->transformContext() );
points << QStringLiteral( "%1,%2,%3,%4,%5,%6,%7,%8" )
.arg( qgsDoubleToString( transformedDestinationPoint.x() ),
qgsDoubleToString( transformedDestinationPoint.y() ),
qgsDoubleToString( sourcePixel.x() ),
qgsDoubleToString( sourcePixel.y() ) )
.arg( pt->isEnabled() )
.arg( qgsDoubleToString( pt->residual().x() ),
qgsDoubleToString( pt->residual().y() ),
qgsDoubleToString( std::sqrt( pt->residual().x() * pt->residual().x() + pt->residual().y() * pt->residual().y() ) ) )
<< endl;
}

mSavedPoints = mPoints.asPoints();
}
else
{
mMessageBar->pushMessage( tr( "Write Error" ), tr( "Could not write to GCP points file %1." ).arg( mGCPpointsFileName ), Qgis::MessageLevel::Critical );
return;
mMessageBar->pushMessage( tr( "Write Error" ), error, Qgis::MessageLevel::Critical );
}

// showMessageInLog(tr("GCP points saved in"), mGCPpointsFileName);
}

QgsGeoreferencerMainWindow::SaveGCPs QgsGeoreferencerMainWindow::checkNeedGCPSave()
Expand Down
2 changes: 1 addition & 1 deletion src/app/georeferencer/qgsgeorefmainwindow.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,7 +164,7 @@ class QgsGeoreferencerMainWindow : public QMainWindow, private Ui::QgsGeorefPlug
void writeSettings();

// gcp points
bool loadGCPs( /*bool verbose = true*/ );
bool loadGCPs( QString &error );
void saveGCPs();
QgsGeoreferencerMainWindow::SaveGCPs checkNeedGCPSave();

Expand Down
63 changes: 63 additions & 0 deletions tests/src/app/testqgsgeoreferencer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@ class TestQgsGeoreferencer : public QObject
void testGcpPoint();
void testGeorefDataPoint();
void testGcpList();
void testSaveLoadGcps();
void testTransformImageNoGeoference();
void testTransformImageWithExistingGeoreference();
void testRasterChangeCoords();
Expand Down Expand Up @@ -220,6 +221,68 @@ void TestQgsGeoreferencer::testGcpList()

}

void TestQgsGeoreferencer::testSaveLoadGcps()
{
// test saving and loading GCPs
QgsGCPList list;
QgsMapCanvas c1;
QgsMapCanvas c2;
list.append( new QgsGeorefDataPoint( &c1, &c2,
QgsPointXY( 111, 222 ), QgsPointXY( -30, 40 ), QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4326" ) ),
true ) );
list.append( new QgsGeorefDataPoint( &c1, &c2,
QgsPointXY( 11, 22 ), QgsPointXY( 16697923, -3503549 ), QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:3857" ) ),
true ) );
// disabled!
list.append( new QgsGeorefDataPoint( &c1, &c2,
QgsPointXY( 33, 44 ), QgsPointXY( 100, 200 ), QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:3857" ) ),
false ) );

QTemporaryDir dir;
QVERIFY( dir.isValid() );
const QString tempFilename = dir.filePath( QStringLiteral( "test.points" ) );

QString error;

QVERIFY( list.saveGcps( tempFilename, QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:3857" ) ),
QgsProject::instance()->transformContext(), error ) );
QVERIFY( error.isEmpty() );


// load
QList<QgsGcpPoint> res = QgsGCPList::loadGcps( QStringLiteral( "not real" ),
QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:3111" ) ),
error );
QVERIFY( !error.isEmpty() );

res = QgsGCPList::loadGcps( tempFilename,
QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:3111" ) ),
error );
QVERIFY( error.isEmpty() );
QCOMPARE( res.size(), 3 );

QCOMPARE( res.at( 0 ).sourcePoint().x(), 111 );
QCOMPARE( res.at( 0 ).sourcePoint().y(), 222 );
QGSCOMPARENEAR( res.at( 0 ).destinationPoint().x(), -3339584, 10000 );
QGSCOMPARENEAR( res.at( 0 ).destinationPoint().y(), 4865942, 10000 );
QVERIFY( res.at( 0 ).isEnabled() );
QCOMPARE( res.at( 0 ).destinationPointCrs().authid(), QStringLiteral( "EPSG:3857" ) );

QCOMPARE( res.at( 1 ).sourcePoint().x(), 11 );
QCOMPARE( res.at( 1 ).sourcePoint().y(), 22 );
QCOMPARE( res.at( 1 ).destinationPoint().x(), 16697923 );
QCOMPARE( res.at( 1 ).destinationPoint().y(), -3503549 );
QVERIFY( res.at( 1 ).isEnabled() );
QCOMPARE( res.at( 1 ).destinationPointCrs().authid(), QStringLiteral( "EPSG:3857" ) );

QCOMPARE( res.at( 2 ).sourcePoint().x(), 33 );
QCOMPARE( res.at( 2 ).sourcePoint().y(), 44 );
QCOMPARE( res.at( 2 ).destinationPoint().x(), 100 );
QCOMPARE( res.at( 2 ).destinationPoint().y(), 200 );
QVERIFY( !res.at( 2 ).isEnabled() );
QCOMPARE( res.at( 2 ).destinationPointCrs().authid(), QStringLiteral( "EPSG:3857" ) );
}

void TestQgsGeoreferencer::testTransformImageNoGeoference()
{
QgsGeorefTransform transform( QgsGcpTransformerInterface::TransformMethod::Linear );
Expand Down

0 comments on commit baf0534

Please sign in to comment.