QGIS API Documentation  3.16.0-Hannover (43b64b13f3)
qgsdatumtransform.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsdatumtransform.cpp
3  ------------------------
4  begin : Dec 2017
5  copyright : (C) 2017 Nyall Dawson
6  email : nyall dot dawson at gmail dot com
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 #include "qgsdatumtransform.h"
19 #include "qgsapplication.h"
20 #include "qgssqliteutils.h"
21 #include <sqlite3.h>
22 
23 #if PROJ_VERSION_MAJOR>=6
24 #include "qgsprojutils.h"
25 #include <proj.h>
26 #endif
27 
28 QList<QgsDatumTransform::TransformDetails> QgsDatumTransform::operations( const QgsCoordinateReferenceSystem &source, const QgsCoordinateReferenceSystem &destination, bool includeSuperseded )
29 {
30  QList< QgsDatumTransform::TransformDetails > res;
31 #if PROJ_VERSION_MAJOR<6
32  Q_UNUSED( source )
33  Q_UNUSED( destination )
34  Q_UNUSED( includeSuperseded )
35 #else
36  if ( !source.projObject() || !destination.projObject() )
37  return res;
38 
39  PJ_CONTEXT *pjContext = QgsProjContext::get();
40 
41  PJ_OPERATION_FACTORY_CONTEXT *operationContext = proj_create_operation_factory_context( pjContext, nullptr );
42 
43  // We want to return ALL grids, not just those available for use
44  proj_operation_factory_context_set_grid_availability_use( pjContext, operationContext, PROJ_GRID_AVAILABILITY_IGNORED );
45 
46  // See https://lists.osgeo.org/pipermail/proj/2019-May/008604.html
47  proj_operation_factory_context_set_spatial_criterion( pjContext, operationContext, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION );
48 
49  if ( includeSuperseded )
50  proj_operation_factory_context_set_discard_superseded( pjContext, operationContext, false );
51 
52  if ( PJ_OBJ_LIST *ops = proj_create_operations( pjContext, source.projObject(), destination.projObject(), operationContext ) )
53  {
54  int count = proj_list_get_count( ops );
55  for ( int i = 0; i < count; ++i )
56  {
57  QgsProjUtils::proj_pj_unique_ptr op( proj_list_get( pjContext, ops, i ) );
58  if ( !op )
59  continue;
60 
61  QgsDatumTransform::TransformDetails details = transformDetailsFromPj( op.get() );
62  if ( !details.proj.isEmpty() )
63  res.push_back( details );
64 
65  }
66  proj_list_destroy( ops );
67  }
68  proj_operation_factory_context_destroy( operationContext );
69 #endif
70  return res;
71 }
72 
73 QList< QgsDatumTransform::TransformPair > QgsDatumTransform::datumTransformations( const QgsCoordinateReferenceSystem &srcCRS, const QgsCoordinateReferenceSystem &destCRS )
74 {
75  QList< QgsDatumTransform::TransformPair > transformations;
76 
77  QString srcGeoId = srcCRS.geographicCrsAuthId();
78  QString destGeoId = destCRS.geographicCrsAuthId();
79 
80  if ( srcGeoId.isEmpty() || destGeoId.isEmpty() )
81  {
82  return transformations;
83  }
84 
85  QStringList srcSplit = srcGeoId.split( ':' );
86  QStringList destSplit = destGeoId.split( ':' );
87 
88  if ( srcSplit.size() < 2 || destSplit.size() < 2 )
89  {
90  return transformations;
91  }
92 
93  int srcAuthCode = srcSplit.at( 1 ).toInt();
94  int destAuthCode = destSplit.at( 1 ).toInt();
95 
96  if ( srcAuthCode == destAuthCode )
97  {
98  return transformations; //crs have the same datum
99  }
100 
101  QList<int> directTransforms;
102  searchDatumTransform( QStringLiteral( "SELECT coord_op_code FROM tbl_datum_transform WHERE source_crs_code=%1 AND target_crs_code=%2 ORDER BY deprecated ASC,preferred DESC" ).arg( srcAuthCode ).arg( destAuthCode ),
103  directTransforms );
104  QList<int> reverseDirectTransforms;
105  searchDatumTransform( QStringLiteral( "SELECT coord_op_code FROM tbl_datum_transform WHERE source_crs_code = %1 AND target_crs_code=%2 ORDER BY deprecated ASC,preferred DESC" ).arg( destAuthCode ).arg( srcAuthCode ),
106  reverseDirectTransforms );
107  QList<int> srcToWgs84;
108  searchDatumTransform( QStringLiteral( "SELECT coord_op_code FROM tbl_datum_transform WHERE (source_crs_code=%1 AND target_crs_code=%2) OR (source_crs_code=%2 AND target_crs_code=%1) ORDER BY deprecated ASC,preferred DESC" ).arg( srcAuthCode ).arg( 4326 ),
109  srcToWgs84 );
110  QList<int> destToWgs84;
111  searchDatumTransform( QStringLiteral( "SELECT coord_op_code FROM tbl_datum_transform WHERE (source_crs_code=%1 AND target_crs_code=%2) OR (source_crs_code=%2 AND target_crs_code=%1) ORDER BY deprecated ASC,preferred DESC" ).arg( destAuthCode ).arg( 4326 ),
112  destToWgs84 );
113 
114  //add direct datum transformations
115  for ( int transform : qgis::as_const( directTransforms ) )
116  {
117  transformations.push_back( QgsDatumTransform::TransformPair( transform, -1 ) );
118  }
119 
120  //add direct datum transformations
121  for ( int transform : qgis::as_const( reverseDirectTransforms ) )
122  {
123  transformations.push_back( QgsDatumTransform::TransformPair( -1, transform ) );
124  }
125 
126  for ( int srcTransform : qgis::as_const( srcToWgs84 ) )
127  {
128  for ( int destTransform : qgis::as_const( destToWgs84 ) )
129  {
130  transformations.push_back( QgsDatumTransform::TransformPair( srcTransform, destTransform ) );
131  }
132  }
133 
134  return transformations;
135 }
136 
137 void QgsDatumTransform::searchDatumTransform( const QString &sql, QList< int > &transforms )
138 {
140  int openResult = database.open_v2( QgsApplication::srsDatabaseFilePath(), SQLITE_OPEN_READONLY, nullptr );
141  if ( openResult != SQLITE_OK )
142  {
143  return;
144  }
145 
147  int prepareRes;
148  statement = database.prepare( sql, prepareRes );
149  if ( prepareRes != SQLITE_OK )
150  {
151  return;
152  }
153 
154  QString cOpCode;
155  while ( statement.step() == SQLITE_ROW )
156  {
157  cOpCode = statement.columnAsText( 0 );
158  transforms.push_back( cOpCode.toInt() );
159  }
160 }
161 
162 QString QgsDatumTransform::datumTransformToProj( int datumTransform )
163 {
164  QString transformString;
165 
167  int openResult = database.open_v2( QgsApplication::srsDatabaseFilePath(), SQLITE_OPEN_READONLY, nullptr );
168  if ( openResult != SQLITE_OK )
169  {
170  return transformString;
171  }
172 
174  QString sql = QStringLiteral( "SELECT coord_op_method_code,p1,p2,p3,p4,p5,p6,p7 FROM tbl_datum_transform WHERE coord_op_code=%1" ).arg( datumTransform );
175  int prepareRes;
176  statement = database.prepare( sql, prepareRes );
177  if ( prepareRes != SQLITE_OK )
178  {
179  return transformString;
180  }
181 
182  if ( statement.step() == SQLITE_ROW )
183  {
184  //coord_op_methode_code
185  int methodCode = statement.columnAsInt64( 0 );
186  if ( methodCode == 9615 ) //ntv2
187  {
188  transformString = "+nadgrids=" + statement.columnAsText( 1 );
189  }
190  else if ( methodCode == 9603 || methodCode == 9606 || methodCode == 9607 )
191  {
192  transformString += QLatin1String( "+towgs84=" );
193  double p1 = statement.columnAsDouble( 1 );
194  double p2 = statement.columnAsDouble( 2 );
195  double p3 = statement.columnAsDouble( 3 );
196  double p4 = statement.columnAsDouble( 4 );
197  double p5 = statement.columnAsDouble( 5 );
198  double p6 = statement.columnAsDouble( 6 );
199  double p7 = statement.columnAsDouble( 7 );
200  if ( methodCode == 9603 ) //3 parameter transformation
201  {
202  transformString += QStringLiteral( "%1,%2,%3" ).arg( QString::number( p1 ), QString::number( p2 ), QString::number( p3 ) );
203  }
204  else //7 parameter transformation
205  {
206  transformString += QStringLiteral( "%1,%2,%3,%4,%5,%6,%7" ).arg( QString::number( p1 ), QString::number( p2 ), QString::number( p3 ), QString::number( p4 ), QString::number( p5 ), QString::number( p6 ), QString::number( p7 ) );
207  }
208  }
209  }
210 
211  return transformString;
212 }
213 
215 {
217  int openResult = database.open_v2( QgsApplication::srsDatabaseFilePath(), SQLITE_OPEN_READONLY, nullptr );
218  if ( openResult != SQLITE_OK )
219  {
220  return -1;
221  }
222 
224  QString sql = QStringLiteral( "SELECT coord_op_method_code,p1,p2,p3,p4,p5,p6,p7,coord_op_code FROM tbl_datum_transform" );
225  int prepareRes;
226  statement = database.prepare( sql, prepareRes );
227  if ( prepareRes != SQLITE_OK )
228  {
229  return -1;
230  }
231 
232  while ( statement.step() == SQLITE_ROW )
233  {
234  QString transformString;
235  //coord_op_methode_code
236  int methodCode = statement.columnAsInt64( 0 );
237  if ( methodCode == 9615 ) //ntv2
238  {
239  transformString = "+nadgrids=" + statement.columnAsText( 1 );
240  }
241  else if ( methodCode == 9603 || methodCode == 9606 || methodCode == 9607 )
242  {
243  transformString += QLatin1String( "+towgs84=" );
244  double p1 = statement.columnAsDouble( 1 );
245  double p2 = statement.columnAsDouble( 2 );
246  double p3 = statement.columnAsDouble( 3 );
247  double p4 = statement.columnAsDouble( 4 );
248  double p5 = statement.columnAsDouble( 5 );
249  double p6 = statement.columnAsDouble( 6 );
250  double p7 = statement.columnAsDouble( 7 );
251  if ( methodCode == 9603 ) //3 parameter transformation
252  {
253  transformString += QStringLiteral( "%1,%2,%3" ).arg( QString::number( p1 ), QString::number( p2 ), QString::number( p3 ) );
254  }
255  else //7 parameter transformation
256  {
257  transformString += QStringLiteral( "%1,%2,%3,%4,%5,%6,%7" ).arg( QString::number( p1 ), QString::number( p2 ), QString::number( p3 ), QString::number( p4 ), QString::number( p5 ), QString::number( p6 ), QString::number( p7 ) );
258  }
259  }
260 
261  if ( transformString.compare( string, Qt::CaseInsensitive ) == 0 )
262  {
263  return statement.columnAsInt64( 8 );
264  }
265  }
266 
267  return -1;
268 }
269 
271 {
273 
275  int openResult = database.open_v2( QgsApplication::srsDatabaseFilePath(), SQLITE_OPEN_READONLY, nullptr );
276  if ( openResult != SQLITE_OK )
277  {
278  return info;
279  }
280 
282  QString sql = QStringLiteral( "SELECT epsg_nr,source_crs_code,target_crs_code,remarks,scope,preferred,deprecated FROM tbl_datum_transform WHERE coord_op_code=%1" ).arg( datumTransform );
283  int prepareRes;
284  statement = database.prepare( sql, prepareRes );
285  if ( prepareRes != SQLITE_OK )
286  {
287  return info;
288  }
289 
290  int srcCrsId, destCrsId;
291  if ( statement.step() != SQLITE_ROW )
292  {
293  return info;
294  }
295 
296  info.datumTransformId = datumTransform;
297  info.epsgCode = statement.columnAsInt64( 0 );
298  srcCrsId = statement.columnAsInt64( 1 );
299  destCrsId = statement.columnAsInt64( 2 );
300  info.remarks = statement.columnAsText( 3 );
301  info.scope = statement.columnAsText( 4 );
302  info.preferred = statement.columnAsInt64( 5 ) != 0;
303  info.deprecated = statement.columnAsInt64( 6 ) != 0;
304 
305  QgsCoordinateReferenceSystem srcCrs = QgsCoordinateReferenceSystem::fromOgcWmsCrs( QStringLiteral( "EPSG:%1" ).arg( srcCrsId ) );
306  info.sourceCrsDescription = srcCrs.description();
307  info.sourceCrsAuthId = srcCrs.authid();
308  QgsCoordinateReferenceSystem destCrs = QgsCoordinateReferenceSystem::fromOgcWmsCrs( QStringLiteral( "EPSG:%1" ).arg( destCrsId ) );
309  info.destinationCrsDescription = destCrs.description();
310  info.destinationCrsAuthId = destCrs.authid();
311 
312  return info;
313 }
314 
315 #if PROJ_VERSION_MAJOR>=6
316 QgsDatumTransform::TransformDetails QgsDatumTransform::transformDetailsFromPj( PJ *op )
317 {
318  PJ_CONTEXT *pjContext = QgsProjContext::get();
319  TransformDetails details;
320  if ( !op )
321  return details;
322 
323  QgsProjUtils::proj_pj_unique_ptr normalized( proj_normalize_for_visualization( pjContext, op ) );
324  if ( normalized )
325  details.proj = QString( proj_as_proj_string( pjContext, normalized.get(), PJ_PROJ_5, nullptr ) );
326 
327  if ( details.proj.isEmpty() )
328  details.proj = QString( proj_as_proj_string( pjContext, op, PJ_PROJ_5, nullptr ) );
329 
330  if ( details.proj.isEmpty() )
331  return details;
332 
333  details.name = QString( proj_get_name( op ) );
334  details.accuracy = proj_coordoperation_get_accuracy( pjContext, op );
335  details.isAvailable = proj_coordoperation_is_instantiable( pjContext, op );
336 
337  details.authority = QString( proj_get_id_auth_name( op, 0 ) );
338  details.code = QString( proj_get_id_code( op, 0 ) );
339 
340  const char *areaOfUseName = nullptr;
341  double westLon = 0;
342  double southLat = 0;
343  double eastLon = 0;
344  double northLat = 0;
345  if ( proj_get_area_of_use( pjContext, op, &westLon, &southLat, &eastLon, &northLat, &areaOfUseName ) )
346  {
347  details.areaOfUse = QString( areaOfUseName );
348  // don't use the constructor which normalizes!
349  details.bounds.setXMinimum( westLon );
350  details.bounds.setYMinimum( southLat );
351  details.bounds.setXMaximum( eastLon );
352  details.bounds.setYMaximum( northLat );
353  }
354 
355  details.remarks = QString( proj_get_remarks( op ) );
356  details.scope = QString( proj_get_scope( op ) );
357 
358  for ( int j = 0; j < proj_coordoperation_get_grid_used_count( pjContext, op ); ++j )
359  {
360  const char *shortName = nullptr;
361  const char *fullName = nullptr;
362  const char *packageName = nullptr;
363  const char *url = nullptr;
364  int directDownload = 0;
365  int openLicense = 0;
366  int isAvailable = 0;
367  proj_coordoperation_get_grid_used( pjContext, op, j, &shortName, &fullName, &packageName, &url, &directDownload, &openLicense, &isAvailable );
368  GridDetails gridDetails;
369  gridDetails.shortName = QString( shortName );
370  gridDetails.fullName = QString( fullName );
371  gridDetails.packageName = QString( packageName );
372  gridDetails.url = QString( url );
373  gridDetails.directDownload = directDownload;
374  gridDetails.openLicense = openLicense;
375  gridDetails.isAvailable = isAvailable;
376 
377  details.grids.append( gridDetails );
378  }
379 
380  if ( proj_get_type( op ) == PJ_TYPE_CONCATENATED_OPERATION )
381  {
382  for ( int j = 0; j < proj_concatoperation_get_step_count( pjContext, op ); ++j )
383  {
384  QgsProjUtils::proj_pj_unique_ptr step( proj_concatoperation_get_step( pjContext, op, j ) );
385  if ( step )
386  {
387  SingleOperationDetails singleOpDetails;
388  singleOpDetails.remarks = QString( proj_get_remarks( step.get() ) );
389  singleOpDetails.scope = QString( proj_get_scope( step.get() ) );
390  singleOpDetails.authority = QString( proj_get_id_auth_name( step.get(), 0 ) );
391  singleOpDetails.code = QString( proj_get_id_code( step.get(), 0 ) );
392 
393  const char *areaOfUseName = nullptr;
394  if ( proj_get_area_of_use( pjContext, step.get(), nullptr, nullptr, nullptr, nullptr, &areaOfUseName ) )
395  {
396  singleOpDetails.areaOfUse = QString( areaOfUseName );
397  }
398  details.operationDetails.append( singleOpDetails );
399  }
400  }
401  }
402  else
403  {
404  SingleOperationDetails singleOpDetails;
405  singleOpDetails.remarks = QString( proj_get_remarks( op ) );
406  singleOpDetails.scope = QString( proj_get_scope( op ) );
407  singleOpDetails.authority = QString( proj_get_id_auth_name( op, 0 ) );
408  singleOpDetails.code = QString( proj_get_id_code( op, 0 ) );
409 
410  const char *areaOfUseName = nullptr;
411  if ( proj_get_area_of_use( pjContext, op, nullptr, nullptr, nullptr, nullptr, &areaOfUseName ) )
412  {
413  singleOpDetails.areaOfUse = QString( areaOfUseName );
414  }
415  details.operationDetails.append( singleOpDetails );
416  }
417 
418  return details;
419 }
420 #endif
QgsProjContext::get
static PJ_CONTEXT * get()
Returns a thread local instance of a proj context, safe for use in the current thread.
Definition: qgsprojutils.cpp:60
QgsDatumTransform::TransformInfo::destinationCrsDescription
QString destinationCrsDescription
Destination CRS description.
Definition: qgsdatumtransform.h:112
QgsCoordinateReferenceSystem::description
QString description() const
Returns the descriptive name of the CRS, e.g., "WGS 84" or "GDA 94 / Vicgrid94".
Definition: qgscoordinatereferencesystem.cpp:1326
sqlite3_database_unique_ptr::prepare
sqlite3_statement_unique_ptr prepare(const QString &sql, int &resultCode) const
Prepares a sql statement, returning the result.
Definition: qgssqliteutils.cpp:99
QgsCoordinateReferenceSystem::fromOgcWmsCrs
static QgsCoordinateReferenceSystem fromOgcWmsCrs(const QString &ogcCrs)
Creates a CRS from a given OGC WMS-format Coordinate Reference System string.
Definition: qgscoordinatereferencesystem.cpp:200
QgsDatumTransform::TransformPair
Contains datum transform information.
Definition: qgsdatumtransform.h:55
sqlite3_statement_unique_ptr::columnAsInt64
qlonglong columnAsInt64(int column) const
Gets column value from the current statement row as a long long integer (64 bits).
Definition: qgssqliteutils.cpp:73
QgsDatumTransform::TransformInfo::sourceCrsAuthId
QString sourceCrsAuthId
Source CRS auth ID.
Definition: qgsdatumtransform.h:103
QgsDatumTransform::TransformInfo::preferred
bool preferred
True if transform is the preferred transform to use for the source/destination CRS combination.
Definition: qgsdatumtransform.h:121
sqlite3_statement_unique_ptr::step
int step()
Steps to the next record in the statement, returning the sqlite3 result code.
Definition: qgssqliteutils.cpp:41
qgsapplication.h
QgsApplication::srsDatabaseFilePath
static QString srsDatabaseFilePath()
Returns the path to the srs.db file.
Definition: qgsapplication.cpp:1018
qgsdatumtransform.h
QgsDatumTransform::datumTransformations
static Q_DECL_DEPRECATED QList< QgsDatumTransform::TransformPair > datumTransformations(const QgsCoordinateReferenceSystem &source, const QgsCoordinateReferenceSystem &destination)
Returns a list of datum transformations which are available for the given source and destination CRS.
Definition: qgsdatumtransform.cpp:73
QgsDatumTransform::TransformInfo::remarks
QString remarks
Transform remarks.
Definition: qgsdatumtransform.h:115
QgsCoordinateReferenceSystem::authid
QString authid() const
Returns the authority identifier for the CRS.
Definition: qgscoordinatereferencesystem.cpp:1321
QgsDatumTransform::TransformInfo::deprecated
bool deprecated
True if transform is deprecated.
Definition: qgsdatumtransform.h:124
QgsDatumTransform::TransformInfo::sourceCrsDescription
QString sourceCrsDescription
Source CRS description.
Definition: qgsdatumtransform.h:109
QgsDatumTransform::TransformInfo::datumTransformId
int datumTransformId
Datum transform ID.
Definition: qgsdatumtransform.h:97
QgsCoordinateReferenceSystem::geographicCrsAuthId
QString geographicCrsAuthId() const
Returns auth id of related geographic CRS.
Definition: qgscoordinatereferencesystem.cpp:3551
sqlite3_database_unique_ptr::open_v2
int open_v2(const QString &path, int flags, const char *zVfs)
Opens the database at the specified file path.
Definition: qgssqliteutils.cpp:86
QgsDatumTransform::TransformDetails
Contains information about a coordinate transformation operation.
Definition: qgsdatumtransform.h:182
QgsCoordinateReferenceSystem
This class represents a coordinate reference system (CRS).
Definition: qgscoordinatereferencesystem.h:206
QgsDatumTransform::TransformInfo::destinationCrsAuthId
QString destinationCrsAuthId
Destination CRS auth ID.
Definition: qgsdatumtransform.h:106
sqlite3_statement_unique_ptr::columnAsText
QString columnAsText(int column) const
Returns the column value from the current statement row as a string.
Definition: qgssqliteutils.cpp:61
QgsDatumTransform::projStringToDatumTransformId
static Q_DECL_DEPRECATED int projStringToDatumTransformId(const QString &string)
Returns the datum transform ID corresponding to a specified proj string.
Definition: qgsdatumtransform.cpp:214
qgsprojutils.h
PJ_CONTEXT
void PJ_CONTEXT
Definition: qgsprojutils.h:151
qgssqliteutils.h
QgsDatumTransform::TransformInfo
Contains datum transform information.
Definition: qgsdatumtransform.h:95
sqlite3_statement_unique_ptr::columnAsDouble
double columnAsDouble(int column) const
Gets column value from the current statement row as a double.
Definition: qgssqliteutils.cpp:51
QgsDatumTransform::TransformDetails::proj
QString proj
Proj representation of transform operation.
Definition: qgsdatumtransform.h:184
sqlite3_database_unique_ptr
Unique pointer for sqlite3 databases, which automatically closes the database when the pointer goes o...
Definition: qgssqliteutils.h:119
QgsDatumTransform::TransformInfo::scope
QString scope
Scope of transform.
Definition: qgsdatumtransform.h:118
QgsDatumTransform::TransformInfo::epsgCode
int epsgCode
EPSG code for the transform, or 0 if not found in EPSG database.
Definition: qgsdatumtransform.h:100
qgscoordinatereferencesystem.h
QgsDatumTransform::datumTransformToProj
static Q_DECL_DEPRECATED QString datumTransformToProj(int datumTransformId)
Returns a proj string representing the specified datumTransformId datum transform ID.
Definition: qgsdatumtransform.cpp:162
QgsDatumTransform::datumTransformInfo
static Q_DECL_DEPRECATED QgsDatumTransform::TransformInfo datumTransformInfo(int datumTransformId)
Returns detailed information about the specified datumTransformId.
Definition: qgsdatumtransform.cpp:270
sqlite3_statement_unique_ptr
Unique pointer for sqlite3 prepared statements, which automatically finalizes the statement when the ...
Definition: qgssqliteutils.h:70
QgsDatumTransform::operations
static QList< QgsDatumTransform::TransformDetails > operations(const QgsCoordinateReferenceSystem &source, const QgsCoordinateReferenceSystem &destination, bool includeSuperseded=false)
Returns a list of coordinate operations available for transforming coordinates from the source to des...
Definition: qgsdatumtransform.cpp:28