QGIS API Documentation  3.9.0-Master (224899f119)
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 PROJ_VERSION_MAJOR>6 || (PROJ_VERSION_MAJOR==6 && PROJ_VERSION_MINOR>=2)
50  if ( includeSuperseded )
51  proj_operation_factory_context_set_discard_superseded( pjContext, operationContext, false );
52 #else
53  Q_UNUSED( includeSuperseded )
54 #endif
55  if ( PJ_OBJ_LIST *ops = proj_create_operations( pjContext, source.projObject(), destination.projObject(), operationContext ) )
56  {
57  int count = proj_list_get_count( ops );
58  for ( int i = 0; i < count; ++i )
59  {
60  QgsProjUtils::proj_pj_unique_ptr op( proj_list_get( pjContext, ops, i ) );
61  if ( !op )
62  continue;
63 
64  res.push_back( transformDetailsFromPj( op.get() ) );
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  details.name = QString( proj_get_name( op ) );
330  details.accuracy = proj_coordoperation_get_accuracy( pjContext, op );
331  details.isAvailable = proj_coordoperation_is_instantiable( pjContext, op );
332 
333  details.authority = QString( proj_get_id_auth_name( op, 0 ) );
334  details.code = QString( proj_get_id_code( op, 0 ) );
335 
336  const char *areaOfUseName = nullptr;
337  double westLon = 0;
338  double southLat = 0;
339  double eastLon = 0;
340  double northLat = 0;
341  if ( proj_get_area_of_use( pjContext, op, &westLon, &southLat, &eastLon, &northLat, &areaOfUseName ) )
342  {
343  details.areaOfUse = QString( areaOfUseName );
344  // don't use the constructor which normalizes!
345  details.bounds.setXMinimum( westLon );
346  details.bounds.setYMinimum( southLat );
347  details.bounds.setXMaximum( eastLon );
348  details.bounds.setYMaximum( northLat );
349  }
350 
351 #if PROJ_VERSION_MAJOR>6 || (PROJ_VERSION_MAJOR==6 && PROJ_VERSION_MINOR>=2)
352  details.remarks = QString( proj_get_remarks( op ) );
353  details.scope = QString( proj_get_scope( op ) );
354 #endif
355 
356  for ( int j = 0; j < proj_coordoperation_get_grid_used_count( pjContext, op ); ++j )
357  {
358  const char *shortName = nullptr;
359  const char *fullName = nullptr;
360  const char *packageName = nullptr;
361  const char *url = nullptr;
362  int directDownload = 0;
363  int openLicense = 0;
364  int isAvailable = 0;
365  proj_coordoperation_get_grid_used( pjContext, op, j, &shortName, &fullName, &packageName, &url, &directDownload, &openLicense, &isAvailable );
366  GridDetails gridDetails;
367  gridDetails.shortName = QString( shortName );
368  gridDetails.fullName = QString( fullName );
369  gridDetails.packageName = QString( packageName );
370  gridDetails.url = QString( url );
371  gridDetails.directDownload = directDownload;
372  gridDetails.openLicense = openLicense;
373  gridDetails.isAvailable = isAvailable;
374 
375  details.grids.append( gridDetails );
376  }
377 
378 #if PROJ_VERSION_MAJOR>6 || (PROJ_VERSION_MAJOR==6 && PROJ_VERSION_MINOR>=2)
379  for ( int j = 0; j < proj_concatoperation_get_step_count( pjContext, op ); ++j )
380  {
381  QgsProjUtils::proj_pj_unique_ptr step( proj_concatoperation_get_step( pjContext, op, j ) );
382  if ( step )
383  {
384  SingleOperationDetails singleOpDetails;
385  singleOpDetails.remarks = QString( proj_get_remarks( step.get() ) );
386  singleOpDetails.scope = QString( proj_get_scope( step.get() ) );
387  singleOpDetails.authority = QString( proj_get_id_auth_name( step.get(), 0 ) );
388  singleOpDetails.code = QString( proj_get_id_code( step.get(), 0 ) );
389 
390  const char *areaOfUseName = nullptr;
391  if ( proj_get_area_of_use( pjContext, step.get(), nullptr, nullptr, nullptr, nullptr, &areaOfUseName ) )
392  {
393  singleOpDetails.areaOfUse = QString( areaOfUseName );
394  }
395  details.operationDetails.append( singleOpDetails );
396  }
397  }
398 #endif
399 
400  return details;
401 }
402 #endif
QList< QgsDatumTransform::SingleOperationDetails > operationDetails
Contains information about the single operation steps used in the transform operation.
QString geographicCrsAuthId() const
Returns auth id of related geographic CRS.
QString remarks
Remarks for operation, from EPSG registry database.
QString name
Display name of transform operation.
QString destinationCrsDescription
Destination CRS description.
void setXMaximum(double x)
Set the maximum x value.
Definition: qgsrectangle.h:135
Unique pointer for sqlite3 prepared statements, which automatically finalizes the statement when the ...
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...
Contains information about a projection transformation grid file.
int datumTransformId
Datum transform ID.
QString packageName
Name of package the grid is included within.
bool deprecated
True if transform is deprecated.
QString scope
Scope of operation, from EPSG registry database.
QString fullName
Full name of transform grid.
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...
bool openLicense
true if grid is available under an open license
QString shortName
Short name of transform grid.
QString sourceCrsDescription
Source CRS description.
int step()
Steps to the next record in the statement, returning the sqlite3 result code.
bool isAvailable
true if grid is currently available for use
void setYMinimum(double y)
Set the minimum y value.
Definition: qgsrectangle.h:140
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.
QString url
Url to download grid from.
QString description() const
Returns the descriptive name of the CRS, e.g., "WGS 84" or "GDA 94 / Vicgrid94".
QString columnAsText(int column) const
Returns the column value from the current statement row as a string.
QString code
Authority code, e.g. "8447" (for EPSG:8447).
QString authority
Authority name, e.g. EPSG.
QString code
Identification code, e.g.
Contains datum transform information.
sqlite3_statement_unique_ptr prepare(const QString &sql, int &resultCode) const
Prepares a sql statement, returning the result.
QString destinationCrsAuthId
Destination CRS auth ID.
QgsRectangle bounds
Valid bounds for the coordinate operation.
bool preferred
True if transform is the preferred transform to use for the source/destination CRS combination...
int epsgCode
EPSG code for the transform, or 0 if not found in EPSG database.
QString sourceCrsAuthId
Source CRS auth ID.
QString authority
Authority name, e.g.
Unique pointer for sqlite3 databases, which automatically closes the database when the pointer goes o...
int open_v2(const QString &path, int flags, const char *zVfs)
Opens the database at the specified file path.
QString areaOfUse
Area of use string.
static QgsCoordinateReferenceSystem fromOgcWmsCrs(const QString &ogcCrs)
Creates a CRS from a given OGC WMS-format Coordinate Reference System string.
Contains datum transform information.
void PJ_CONTEXT
Definition: qgsprojutils.h:135
QList< QgsDatumTransform::GridDetails > grids
Contains a list of transform grids used by the operation.
void setYMaximum(double y)
Set the maximum y value.
Definition: qgsrectangle.h:145
This class represents a coordinate reference system (CRS).
bool isAvailable
true if operation is available.
static Q_DECL_DEPRECATED QString datumTransformToProj(int datumTransformId)
Returns a proj string representing the specified datumTransformId datum transform ID...
static QString srsDatabaseFilePath()
Returns the path to the srs.db file.
double accuracy
Transformation accuracy (in meters)
static Q_DECL_DEPRECATED int projStringToDatumTransformId(const QString &string)
Returns the datum transform ID corresponding to a specified proj string.
static PJ_CONTEXT * get()
Returns a thread local instance of a proj context, safe for use in the current thread.
bool directDownload
true if direct download of grid is possible
QString remarks
Transform remarks.
Contains information about a coordinate transformation operation.
QString proj
Proj representation of transform operation.
static Q_DECL_DEPRECATED QgsDatumTransform::TransformInfo datumTransformInfo(int datumTransformId)
Returns detailed information about the specified datumTransformId.
Contains information about a single coordinate operation.
double columnAsDouble(int column) const
Gets column value from the current statement row as a double.
qlonglong columnAsInt64(int column) const
Gets column value from the current statement row as a long long integer (64 bits).
QString scope
Scope of transform.
QString authid() const
Returns the authority identifier for the CRS.
void setXMinimum(double x)
Set the minimum x value.
Definition: qgsrectangle.h:130