QGIS API Documentation 3.38.0-Grenoble (exported)
Loading...
Searching...
No Matches
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
26QList<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
65QList< 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
129void 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
154QString 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.
PJ * projObject() const
Returns the underlying PROJ PJ object corresponding to the CRS, or nullptr if the CRS is invalid.
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.
void setYMinimum(double y)
Set the minimum y value.
void setXMinimum(double x)
Set the minimum x value.
void setYMaximum(double y)
Set the maximum y value.
void setXMaximum(double x)
Set the maximum x value.
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 projCtx_t PJ_CONTEXT
struct PJconsts PJ
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.