Skip to content
Permalink
Browse files

[FEATURE] selection of datum transformation for otf-reprojection

  • Loading branch information
mhugent committed Nov 8, 2013
2 parents b21adc1 + 6aa4d95 commit 3a14b5384977897de5ac35b74e89d1ae821d212e
@@ -245,6 +245,8 @@ class QgsMapRenderer : QObject
//! Added in QGIS v1.4
void setLabelingEngine( QgsLabelingEngineInterface* iface /Transfer/ );

const QgsCoordinateTransform* transformation( const QgsMapLayer *layer ) const;

signals:

void drawingProgress( int current, int total );
@@ -28,6 +28,29 @@ class QgsVectorFileWriter
ErrInvalidLayer, // added in 2.0
};

//added in 2.0
enum SymbologyExport
{
NoSymbology = 0, //export only data
FeatureSymbology, //Keeps the number of features and export symbology per feature
SymbolLayerSymbology //Exports one feature per symbol layer (considering symbol levels)
};

static WriterError writeAsVectorFormat( QgsVectorLayer* layer,
const QString& fileName,
const QString& fileEncoding,
const QgsCoordinateTransform* ct,
const QString& driverName = "ESRI Shapefile",
bool onlySelected = false,
QString *errorMessage = 0,
const QStringList &datasourceOptions = QStringList(), // added in 1.6
const QStringList &layerOptions = QStringList(), // added in 1.6
bool skipAttributeCreation = false, // added in 1.6
QString *newFilename = 0, // added in 1.9
SymbologyExport symbologyExport = NoSymbology, //added in 2.0
double symbologyScale = 1.0 // added in 2.0
);

/** Write contents of vector layer to an (OGR supported) vector formt
@note: this method was added in version 1.5*/
static WriterError writeAsVectorFormat( QgsVectorLayer* layer,
@@ -40,7 +63,9 @@ class QgsVectorFileWriter
const QStringList &datasourceOptions = QStringList(), // added in 1.6
const QStringList &layerOptions = QStringList(), // added in 1.6
bool skipAttributeCreation = false, // added in 1.6
QString *newFilename = 0 // added in 1.9
QString *newFilename = 0, // added in 1.9
SymbologyExport symbologyExport = NoSymbology, //added in 2.0
double symbologyScale = 1.0 // added in 2.0
);

/** create shapefile and initialize it */
BIN +37 KB (100%) resources/srs.db
Binary file not shown.
@@ -118,6 +118,7 @@
#include "qgscustomization.h"
#include "qgscustomprojectiondialog.h"
#include "qgsdatasourceuri.h"
#include "qgsdatumtransformdialog.h"
#include "qgsdecorationcopyright.h"
#include "qgsdecorationnortharrow.h"
#include "qgsdecorationscalebar.h"
@@ -4611,10 +4612,13 @@ void QgisApp::saveAsVectorFileGeneral( bool saveOnlySelection, QgsVectorLayer* v
QString format = dialog->format();
QStringList datasourceOptions = dialog->datasourceOptions();

QgsCoordinateTransform* ct = 0;

switch ( dialog->crs() )
{
case -2: // Project CRS
destCRS = mMapCanvas->mapRenderer()->destinationCrs();

break;
case -1: // Layer CRS
destCRS = vlayer->crs();
@@ -4625,14 +4629,39 @@ void QgisApp::saveAsVectorFileGeneral( bool saveOnlySelection, QgsVectorLayer* v
break;
}

if ( destCRS.isValid() && destCRS != vlayer->crs() )
{
ct = new QgsCoordinateTransform( vlayer->crs(), destCRS );

//ask user about datum transformation
QList< QList< int > > dt = QgsCoordinateTransform::datumTransformations( vlayer->crs(), destCRS );
if ( dt.size() > 1 )
{
QgsDatumTransformDialog d( vlayer->name(), dt );
if ( d.exec() == QDialog::Accepted )
{
QList< int > sdt = d.selectedDatumTransform();
if ( sdt.size() > 0 )
{
ct->setSourceDatumTransform( sdt.at( 0 ) );
}
if ( sdt.size() > 1 )
{
ct->setDestinationDatumTransform( sdt.at( 1 ) );
}
ct->initialise();
}
}
}

// ok if the file existed it should be deleted now so we can continue...
QApplication::setOverrideCursor( Qt::WaitCursor );

QgsVectorFileWriter::WriterError error;
QString errorMessage;
QString newFilename;
error = QgsVectorFileWriter::writeAsVectorFormat(
vlayer, vectorFilename, encoding, &destCRS, format,
vlayer, vectorFilename, encoding, ct, format,
saveOnlySelection,
&errorMessage,
datasourceOptions, dialog->layerOptions(),
@@ -4641,6 +4670,8 @@ void QgisApp::saveAsVectorFileGeneral( bool saveOnlySelection, QgsVectorLayer* v
( QgsVectorFileWriter::SymbologyExport )( dialog->symbologyExport() ),
dialog->scaleDenominator() );

delete ct;

QApplication::restoreOverrideCursor();

if ( error == QgsVectorFileWriter::NoError )
@@ -1677,10 +1677,13 @@ bool QgsCoordinateReferenceSystem::loadIDs( QHash<int, QString> &wkts )

int QgsCoordinateReferenceSystem::syncDb()
{
QString dbFilePath = QgsApplication::srsDbFilePath();
syncDatumTransform( dbFilePath );

int inserted = 0, updated = 0, deleted = 0, errors = 0;

sqlite3 *database;
if ( sqlite3_open( QgsApplication::srsDbFilePath().toUtf8().constData(), &database ) != SQLITE_OK )
if ( sqlite3_open( dbFilePath.toUtf8().constData(), &database ) != SQLITE_OK )
{
qCritical( "Could not open database: %s [%s]\n", QgsApplication::srsDbFilePath().toLocal8Bit().constData(), sqlite3_errmsg( database ) );
return -1;
@@ -1933,3 +1936,125 @@ int QgsCoordinateReferenceSystem::syncDb()
else
return updated + inserted;
}

bool QgsCoordinateReferenceSystem::syncDatumTransform( const QString& dbPath )
{
QString filename = CPLFindFile( "gdal", "datum_shift.csv" );

QFile f( filename );
if ( !f.open( QIODevice::ReadOnly ) )
{
return false;
}

sqlite3* db;
int openResult = sqlite3_open( dbPath.toUtf8().constData(), &db );
if ( openResult != SQLITE_OK )
{
return false;
}

if ( sqlite3_exec( db, "BEGIN TRANSACTION", 0, 0, 0 ) != SQLITE_OK )
{
qCritical( "Could not begin transaction: %s [%s]\n", QgsApplication::srsDbFilePath().toLocal8Bit().constData(), sqlite3_errmsg( db ) );
return false;
}


QTextStream textStream( &f );
textStream.readLine();

QString line, coord_op, source_crs, target_crs, coord_op_method,
p1, p2, p3, p4, p5, p6, p7;

while ( !textStream.atEnd() )
{
line = textStream.readLine();
QStringList csList = line.split( "," );
int csSize = csList.size();
if ( csSize < 22 )
{
continue;
}

coord_op = csList[1];
source_crs = csList[2];
target_crs = csList[3];
coord_op_method = csList[csSize - 9];
p1 = csList[csSize - 8];
p1 = p1.isEmpty() ? "NULL" : p1;
p2 = csList[csSize - 7];
p2 = p2.isEmpty() ? "NULL" : p2;
p3 = csList[csSize - 6];
p3 = p3.isEmpty() ? "NULL" : p3;
p4 = csList[csSize - 5];
p4 = p4.isEmpty() ? "NULL" : p4;
p5 = csList[csSize - 4];
p5 = p5.isEmpty() ? "NULL" : p5;
p6 = csList[csSize - 3];
p6 = p6.isEmpty() ? "NULL" : p6;
p7 = csList[csSize - 2];
p7 = p7.isEmpty() ? "NULL" : p7;

//entry already in db?
sqlite3_stmt* stmt;
QString cOpCode;
QString sql = QString( "SELECT coord_op_code FROM tbl_datum_transform WHERE coord_op_code=%1" ).arg( coord_op );
int prepareRes = sqlite3_prepare( db, sql.toAscii(), sql.size(), &stmt, NULL );
if ( prepareRes != SQLITE_OK )
{
continue;
}

if ( sqlite3_step( stmt ) == SQLITE_ROW )
{
cOpCode = ( const char * ) sqlite3_column_text( stmt, 0 );
}
sqlite3_finalize( stmt );

if ( !cOpCode.isEmpty() )
{
//already in database, do update
QgsDebugMsg( "Trying datum transform update" );
sql = QString( "UPDATE tbl_datum_transform SET source_crs = %2, target_crs = %3, coord_op_method = %4, p1 = %5, p2 = %6, p3 = %7, p4 = %8, p5 = %9, p6 = %10, p7 = %11 WHERE coord_op = %1" )
.arg( coord_op ).arg( source_crs ).arg( target_crs ).arg( coord_op_method ).arg( p1 ).arg( p2 ).arg( p3 ).arg( p4 ).arg( p5 ).arg( p6 ).arg( p7 );
}
{
//not yet in database, do insert
QgsDebugMsg( "Trying datum transform insert" );
sql = QString( "INSERT INTO tbl_datum_transform ( coord_op_code, source_crs_code, target_crs_code, coord_op_method_code, p1, p2, p3, p4, p5, p6, p7 ) VALUES ( %1, %2, %3, %4, %5, %6, %7, %8, %9, %10, %11 )" )
.arg( coord_op ).arg( source_crs ).arg( target_crs ).arg( coord_op_method ).arg( p1 ).arg( p2 ).arg( p3 ).arg( p4 ).arg( p5 ).arg( p6 ).arg( p7 );

}

if ( sqlite3_exec( db, sql.toUtf8(), 0, 0, 0 ) != SQLITE_OK )
{
QgsDebugMsg( "Error" );
}
}

if ( sqlite3_exec( db, "COMMIT", 0, 0, 0 ) != SQLITE_OK )
{
qCritical( "Could not commit transaction: %s [%s]\n", QgsApplication::srsDbFilePath().toLocal8Bit().constData(), sqlite3_errmsg( db ) );
return false;
}

sqlite3_close( db );
return true;
}

QString QgsCoordinateReferenceSystem::geographicCRSAuthId() const
{
if ( geographicFlag() )
{
return mAuthId;
}
else if ( mCRS )
{
return OSRGetAuthorityName( mCRS, "GEOGCS" ) + QString( ":" ) + OSRGetAuthorityCode( mCRS, "GEOGCS" );
}
else
{
return "";
}
}
@@ -354,6 +354,9 @@ class CORE_EXPORT QgsCoordinateReferenceSystem
*/
bool saveAsUserCRS( QString name );

/**Returns auth id of related geographic CRS*/
QString geographicCRSAuthId() const;

// Mutators -----------------------------------
// We don't want to expose these to the public api since they wont create
// a fully valid crs. Programmers should use the createFrom* methods rather
@@ -465,6 +468,7 @@ class CORE_EXPORT QgsCoordinateReferenceSystem

static bool loadIDs( QHash<int, QString> &wkts );
static bool loadWkts( QHash<int, QString> &wkts, const char *filename );
static bool syncDatumTransform( const QString& dbPath );

//!Whether this is a coordinate system has inverted axis
mutable int mAxisInverted;

12 comments on commit 3a14b53

@palmerj

This comment has been minimized.

Copy link
Contributor

@palmerj palmerj replied Nov 8, 2013

Can we please add the NZGD49<->NZGD2000 (~WGS84) datum grid shift transformation as well? I think this SQL should do it:

INSERT INTO tbl_datum_transform (
coord_op_code,
source_crs_code,
target_crs_code,
coord_op_method_code,
p1
)
VALUES (
100002,
4272,
4326,
9615,
'nzgd2kgrid0005.gsb'
);

For more information bout this transformation see here: http://www.linz.govt.nz/geodetic/conversion-coordinates/geodetic-datum-conversion/nzgd1949-nzgd2000/distortion-grid

@gioman

This comment has been minimized.

Copy link
Contributor

@gioman gioman replied Nov 9, 2013

I have also grids and transformations for Portugal to suggest to be added. What is the right way to do it?

@palmerj

This comment has been minimized.

Copy link
Contributor

@palmerj palmerj replied Nov 9, 2013

This is is really great functionality and something I've been waiting a long time for - many thanks Marco!!

A few notes from my testing so far:

  • It would also be good to improve the datum transformation method dialogue to only have the different methods listed once. Currently for the EPSG:4272 <-> EPSG:4167 it displays the 3 different methods 4 times:

datum_transform

  • Would be nice to have the dialogue stating the transformation method as a column. i.e: 3 parameter, 7 parameter, Grid shift
  • Somewhere on the dialogue some text stating the datums the shift applies to. i.e in my example "NZGD49 / New Zealand Map Grid" to "NZGD2000"
  • Would be really nice to have an option to set the default transformation method between datums in the QGIS application. Maybe this dialogue would give you a check box to set it as the default. If the user wants to change it latter then they would need to go the QGIS settings to delete or change it.
@mhugent

This comment has been minimized.

Copy link
Contributor Author

@mhugent mhugent replied Nov 11, 2013

@palmerj: Thanks, the entry for NZGD49<->NZGD2000 is in srs.db now.
Remembering default transformations sounds reasonable. I'm trying to add that in the options dialog.

@mhugent

This comment has been minimized.

Copy link
Contributor Author

@mhugent mhugent replied Nov 11, 2013

@gioman: just send me the defails for the ntv2-Transformation that should go to srs.db (similar to the lines Jeremy sent).

@palmerj

This comment has been minimized.

Copy link
Contributor

@palmerj palmerj replied Nov 11, 2013

@mhugent thanks for the dialogue improvement to remember the default transformations.

Any chance to stop the duplicates showing in the dialogue as per my screenshot in #commitcomment-4556980. I think the bug is in the QgsCoordinateTransform::datumTransformations method.

@mhugent

This comment has been minimized.

Copy link
Contributor Author

@mhugent mhugent replied Nov 11, 2013

@palmerj
Yes, sounds like there might be something wrong in that method. Though I have difficulties here to reproduce the duplicated entries. Also with the NZ coordinate systems, I get only three entries plus the ntv2 one. Another possibility would be that something went wrong with the datum transform sync and produced duplicated entry in tbl_datum_transrom

@dakcarto

This comment has been minimized.

Copy link
Member

@dakcarto dakcarto replied Nov 12, 2013

@mhugent thank you for this work! however, on Mac when building WITH_TESTS I get ~2,000 errors when trying to do:

QgsCoordinateReferenceSystem::syncDatumTransform( const QString& dbPath )
3a14b53#diff-8ddbdd08af8ef2dce21153a7c23a4b55R2032

@etiennesky

This comment has been minimized.

Copy link
Contributor

@etiennesky etiennesky replied Nov 12, 2013

same here in linux debug build

Synchronizing CRS database with GDAL/PROJ definitions.
rc/core/qgscoordinatereferencesystem.cpp: 2018: (syncDatumTransform) Trying datum transform update
rc/core/qgscoordinatereferencesystem.cpp: 2024: (syncDatumTransform) Trying datum transform insert
rc/core/qgscoordinatereferencesystem.cpp: 2032: (syncDatumTransform) Error
rc/core/qgscoordinatereferencesystem.cpp: 2018: (syncDatumTransform) Trying datum transform update
rc/core/qgscoordinatereferencesystem.cpp: 2024: (syncDatumTransform) Trying datum transform insert
rc/core/qgscoordinatereferencesystem.cpp: 2032: (syncDatumTransform) Error
...

@mhugent

This comment has been minimized.

Copy link
Contributor Author

@mhugent mhugent replied Nov 12, 2013

Oh, what a stupid bug... Please try now, it should be fixed.
This bug might also be the cause of the duplicated entries mentioned by Jeremy. No idea why that worked on my system until now.

@dakcarto

This comment has been minimized.

Copy link
Member

@dakcarto dakcarto replied Nov 13, 2013

@mhugent, I was still having an issue (minus the debug output on insert):

rc/core/qgscoordinatereferencesystem.cpp: 2018: (syncDatumTransform) Trying datum transform update
rc/core/qgscoordinatereferencesystem.cpp: 2033: (syncDatumTransform) Error
...

Fixed by appending missing '_code' to UPDATE column names

After fix, still got 700+ debug lines during compile (these show up as errors in Qt Creator), when building CMAKE_BUILD_TYPE=RelWithDebInfo:
rc/core/qgscoordinatereferencesystem.cpp: 2018: (syncDatumTransform) Trying datum transform update

So maybe have those SQL commands not produce debug output unless the message level is higher? e.g.:
QgsDebugMsgLevel( 'Trying datum transform update', 4 );

@palmerj

This comment has been minimized.

Copy link
Contributor

@palmerj palmerj replied Nov 13, 2013

@mhugent I also get the error as reported by @dakcarto however my duplicate problem has now been fixed. Thanks!

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