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