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"
18 #include "qgscoordinatereferencesystem.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
operations(const QgsCoordinateReferenceSystem & source,const QgsCoordinateReferenceSystem & destination,bool includeSuperseded)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
datumTransformations(const QgsCoordinateReferenceSystem & srcCRS,const QgsCoordinateReferenceSystem & destCRS)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
searchDatumTransform(const QString & sql,QList<int> & transforms)137 void QgsDatumTransform::searchDatumTransform( const QString &sql, QList< int > &transforms )
138 {
139 sqlite3_database_unique_ptr database;
140 int openResult = database.open_v2( QgsApplication::srsDatabaseFilePath(), SQLITE_OPEN_READONLY, nullptr );
141 if ( openResult != SQLITE_OK )
142 {
143 return;
144 }
145
146 sqlite3_statement_unique_ptr statement;
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
datumTransformToProj(int datumTransform)162 QString QgsDatumTransform::datumTransformToProj( int datumTransform )
163 {
164 QString transformString;
165
166 sqlite3_database_unique_ptr database;
167 int openResult = database.open_v2( QgsApplication::srsDatabaseFilePath(), SQLITE_OPEN_READONLY, nullptr );
168 if ( openResult != SQLITE_OK )
169 {
170 return transformString;
171 }
172
173 sqlite3_statement_unique_ptr statement;
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
projStringToDatumTransformId(const QString & string)214 int QgsDatumTransform::projStringToDatumTransformId( const QString &string )
215 {
216 sqlite3_database_unique_ptr database;
217 int openResult = database.open_v2( QgsApplication::srsDatabaseFilePath(), SQLITE_OPEN_READONLY, nullptr );
218 if ( openResult != SQLITE_OK )
219 {
220 return -1;
221 }
222
223 sqlite3_statement_unique_ptr statement;
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
datumTransformInfo(int datumTransform)270 QgsDatumTransform::TransformInfo QgsDatumTransform::datumTransformInfo( int datumTransform )
271 {
272 QgsDatumTransform::TransformInfo info;
273
274 sqlite3_database_unique_ptr database;
275 int openResult = database.open_v2( QgsApplication::srsDatabaseFilePath(), SQLITE_OPEN_READONLY, nullptr );
276 if ( openResult != SQLITE_OK )
277 {
278 return info;
279 }
280
281 sqlite3_statement_unique_ptr statement;
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
transformDetailsFromPj(PJ * op)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
421