Branch data Line data Source code
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 : : #include "qgsprojutils.h"
24 : : #include <proj.h>
25 : :
26 : 0 : QList<QgsDatumTransform::TransformDetails> QgsDatumTransform::operations( const QgsCoordinateReferenceSystem &source, const QgsCoordinateReferenceSystem &destination, bool includeSuperseded )
27 : : {
28 : 0 : QList< QgsDatumTransform::TransformDetails > res;
29 : 0 : if ( !source.projObject() || !destination.projObject() )
30 : 0 : return res;
31 : :
32 : 0 : PJ_CONTEXT *pjContext = QgsProjContext::get();
33 : :
34 : 0 : 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 : 0 : 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 : 0 : proj_operation_factory_context_set_spatial_criterion( pjContext, operationContext, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION );
41 : :
42 : 0 : if ( includeSuperseded )
43 : 0 : proj_operation_factory_context_set_discard_superseded( pjContext, operationContext, false );
44 : :
45 : 0 : if ( PJ_OBJ_LIST *ops = proj_create_operations( pjContext, source.projObject(), destination.projObject(), operationContext ) )
46 : : {
47 : 0 : int count = proj_list_get_count( ops );
48 : 0 : for ( int i = 0; i < count; ++i )
49 : : {
50 : 0 : QgsProjUtils::proj_pj_unique_ptr op( proj_list_get( pjContext, ops, i ) );
51 : 0 : if ( !op )
52 : 0 : continue;
53 : :
54 : 0 : QgsDatumTransform::TransformDetails details = transformDetailsFromPj( op.get() );
55 : 0 : if ( !details.proj.isEmpty() )
56 : 0 : res.push_back( details );
57 : :
58 : 0 : }
59 : 0 : proj_list_destroy( ops );
60 : 0 : }
61 : 0 : proj_operation_factory_context_destroy( operationContext );
62 : 0 : return res;
63 : 0 : }
64 : :
65 : 0 : QList< QgsDatumTransform::TransformPair > QgsDatumTransform::datumTransformations( const QgsCoordinateReferenceSystem &srcCRS, const QgsCoordinateReferenceSystem &destCRS )
66 : : {
67 : 0 : QList< QgsDatumTransform::TransformPair > transformations;
68 : :
69 : 0 : QString srcGeoId = srcCRS.geographicCrsAuthId();
70 : 0 : QString destGeoId = destCRS.geographicCrsAuthId();
71 : :
72 : 0 : if ( srcGeoId.isEmpty() || destGeoId.isEmpty() )
73 : : {
74 : 0 : return transformations;
75 : : }
76 : :
77 : 0 : QStringList srcSplit = srcGeoId.split( ':' );
78 : 0 : QStringList destSplit = destGeoId.split( ':' );
79 : :
80 : 0 : if ( srcSplit.size() < 2 || destSplit.size() < 2 )
81 : : {
82 : 0 : return transformations;
83 : : }
84 : :
85 : 0 : int srcAuthCode = srcSplit.at( 1 ).toInt();
86 : 0 : int destAuthCode = destSplit.at( 1 ).toInt();
87 : :
88 : 0 : if ( srcAuthCode == destAuthCode )
89 : : {
90 : 0 : return transformations; //crs have the same datum
91 : : }
92 : :
93 : 0 : QList<int> directTransforms;
94 : 0 : 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 : 0 : QList<int> reverseDirectTransforms;
97 : 0 : 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 : 0 : QList<int> srcToWgs84;
100 : 0 : 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 : 0 : QList<int> destToWgs84;
103 : 0 : 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 : 0 : for ( int transform : std::as_const( directTransforms ) )
108 : : {
109 : 0 : transformations.push_back( QgsDatumTransform::TransformPair( transform, -1 ) );
110 : : }
111 : :
112 : : //add direct datum transformations
113 : 0 : for ( int transform : std::as_const( reverseDirectTransforms ) )
114 : : {
115 : 0 : transformations.push_back( QgsDatumTransform::TransformPair( -1, transform ) );
116 : : }
117 : :
118 : 0 : for ( int srcTransform : std::as_const( srcToWgs84 ) )
119 : : {
120 : 0 : for ( int destTransform : std::as_const( destToWgs84 ) )
121 : : {
122 : 0 : transformations.push_back( QgsDatumTransform::TransformPair( srcTransform, destTransform ) );
123 : : }
124 : : }
125 : :
126 : 0 : return transformations;
127 : 0 : }
128 : :
129 : 0 : void QgsDatumTransform::searchDatumTransform( const QString &sql, QList< int > &transforms )
130 : : {
131 : 0 : sqlite3_database_unique_ptr database;
132 : 0 : int openResult = database.open_v2( QgsApplication::srsDatabaseFilePath(), SQLITE_OPEN_READONLY, nullptr );
133 : 0 : if ( openResult != SQLITE_OK )
134 : : {
135 : 0 : return;
136 : : }
137 : :
138 : 0 : sqlite3_statement_unique_ptr statement;
139 : : int prepareRes;
140 : 0 : statement = database.prepare( sql, prepareRes );
141 : 0 : if ( prepareRes != SQLITE_OK )
142 : : {
143 : 0 : return;
144 : : }
145 : :
146 : 0 : QString cOpCode;
147 : 0 : while ( statement.step() == SQLITE_ROW )
148 : : {
149 : 0 : cOpCode = statement.columnAsText( 0 );
150 : 0 : transforms.push_back( cOpCode.toInt() );
151 : : }
152 : 0 : }
153 : :
154 : 0 : QString QgsDatumTransform::datumTransformToProj( int datumTransform )
155 : : {
156 : 0 : QString transformString;
157 : :
158 : 0 : sqlite3_database_unique_ptr database;
159 : 0 : int openResult = database.open_v2( QgsApplication::srsDatabaseFilePath(), SQLITE_OPEN_READONLY, nullptr );
160 : 0 : if ( openResult != SQLITE_OK )
161 : : {
162 : 0 : return transformString;
163 : : }
164 : :
165 : 0 : sqlite3_statement_unique_ptr statement;
166 : 0 : 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 : 0 : statement = database.prepare( sql, prepareRes );
169 : 0 : if ( prepareRes != SQLITE_OK )
170 : : {
171 : 0 : return transformString;
172 : : }
173 : :
174 : 0 : if ( statement.step() == SQLITE_ROW )
175 : : {
176 : : //coord_op_methode_code
177 : 0 : int methodCode = statement.columnAsInt64( 0 );
178 : 0 : if ( methodCode == 9615 ) //ntv2
179 : : {
180 : 0 : transformString = "+nadgrids=" + statement.columnAsText( 1 );
181 : 0 : }
182 : 0 : else if ( methodCode == 9603 || methodCode == 9606 || methodCode == 9607 )
183 : : {
184 : 0 : transformString += QLatin1String( "+towgs84=" );
185 : 0 : double p1 = statement.columnAsDouble( 1 );
186 : 0 : double p2 = statement.columnAsDouble( 2 );
187 : 0 : double p3 = statement.columnAsDouble( 3 );
188 : 0 : double p4 = statement.columnAsDouble( 4 );
189 : 0 : double p5 = statement.columnAsDouble( 5 );
190 : 0 : double p6 = statement.columnAsDouble( 6 );
191 : 0 : double p7 = statement.columnAsDouble( 7 );
192 : 0 : if ( methodCode == 9603 ) //3 parameter transformation
193 : : {
194 : 0 : transformString += QStringLiteral( "%1,%2,%3" ).arg( QString::number( p1 ), QString::number( p2 ), QString::number( p3 ) );
195 : 0 : }
196 : : else //7 parameter transformation
197 : : {
198 : 0 : 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 : 0 : }
201 : 0 : }
202 : :
203 : 0 : return transformString;
204 : 0 : }
205 : :
206 : 0 : int QgsDatumTransform::projStringToDatumTransformId( const QString &string )
207 : : {
208 : 0 : sqlite3_database_unique_ptr database;
209 : 0 : int openResult = database.open_v2( QgsApplication::srsDatabaseFilePath(), SQLITE_OPEN_READONLY, nullptr );
210 : 0 : if ( openResult != SQLITE_OK )
211 : : {
212 : 0 : return -1;
213 : : }
214 : :
215 : 0 : sqlite3_statement_unique_ptr statement;
216 : 0 : 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 : 0 : statement = database.prepare( sql, prepareRes );
219 : 0 : if ( prepareRes != SQLITE_OK )
220 : : {
221 : 0 : return -1;
222 : : }
223 : :
224 : 0 : while ( statement.step() == SQLITE_ROW )
225 : : {
226 : 0 : QString transformString;
227 : : //coord_op_methode_code
228 : 0 : int methodCode = statement.columnAsInt64( 0 );
229 : 0 : if ( methodCode == 9615 ) //ntv2
230 : : {
231 : 0 : transformString = "+nadgrids=" + statement.columnAsText( 1 );
232 : 0 : }
233 : 0 : else if ( methodCode == 9603 || methodCode == 9606 || methodCode == 9607 )
234 : : {
235 : 0 : transformString += QLatin1String( "+towgs84=" );
236 : 0 : double p1 = statement.columnAsDouble( 1 );
237 : 0 : double p2 = statement.columnAsDouble( 2 );
238 : 0 : double p3 = statement.columnAsDouble( 3 );
239 : 0 : double p4 = statement.columnAsDouble( 4 );
240 : 0 : double p5 = statement.columnAsDouble( 5 );
241 : 0 : double p6 = statement.columnAsDouble( 6 );
242 : 0 : double p7 = statement.columnAsDouble( 7 );
243 : 0 : if ( methodCode == 9603 ) //3 parameter transformation
244 : : {
245 : 0 : transformString += QStringLiteral( "%1,%2,%3" ).arg( QString::number( p1 ), QString::number( p2 ), QString::number( p3 ) );
246 : 0 : }
247 : : else //7 parameter transformation
248 : : {
249 : 0 : 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 : 0 : }
252 : :
253 : 0 : if ( transformString.compare( string, Qt::CaseInsensitive ) == 0 )
254 : : {
255 : 0 : return statement.columnAsInt64( 8 );
256 : : }
257 : 0 : }
258 : :
259 : 0 : return -1;
260 : 0 : }
261 : :
262 : 0 : QgsDatumTransform::TransformInfo QgsDatumTransform::datumTransformInfo( int datumTransform )
263 : : {
264 : 0 : QgsDatumTransform::TransformInfo info;
265 : :
266 : 0 : sqlite3_database_unique_ptr database;
267 : 0 : int openResult = database.open_v2( QgsApplication::srsDatabaseFilePath(), SQLITE_OPEN_READONLY, nullptr );
268 : 0 : if ( openResult != SQLITE_OK )
269 : : {
270 : 0 : return info;
271 : : }
272 : :
273 : 0 : sqlite3_statement_unique_ptr statement;
274 : 0 : 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 : 0 : statement = database.prepare( sql, prepareRes );
277 : 0 : if ( prepareRes != SQLITE_OK )
278 : : {
279 : 0 : return info;
280 : : }
281 : :
282 : : int srcCrsId, destCrsId;
283 : 0 : if ( statement.step() != SQLITE_ROW )
284 : : {
285 : 0 : return info;
286 : : }
287 : :
288 : 0 : info.datumTransformId = datumTransform;
289 : 0 : info.epsgCode = statement.columnAsInt64( 0 );
290 : 0 : srcCrsId = statement.columnAsInt64( 1 );
291 : 0 : destCrsId = statement.columnAsInt64( 2 );
292 : 0 : info.remarks = statement.columnAsText( 3 );
293 : 0 : info.scope = statement.columnAsText( 4 );
294 : 0 : info.preferred = statement.columnAsInt64( 5 ) != 0;
295 : 0 : info.deprecated = statement.columnAsInt64( 6 ) != 0;
296 : :
297 : 0 : QgsCoordinateReferenceSystem srcCrs = QgsCoordinateReferenceSystem::fromOgcWmsCrs( QStringLiteral( "EPSG:%1" ).arg( srcCrsId ) );
298 : 0 : info.sourceCrsDescription = srcCrs.description();
299 : 0 : info.sourceCrsAuthId = srcCrs.authid();
300 : 0 : QgsCoordinateReferenceSystem destCrs = QgsCoordinateReferenceSystem::fromOgcWmsCrs( QStringLiteral( "EPSG:%1" ).arg( destCrsId ) );
301 : 0 : info.destinationCrsDescription = destCrs.description();
302 : 0 : info.destinationCrsAuthId = destCrs.authid();
303 : :
304 : 0 : return info;
305 : 0 : }
306 : :
307 : 0 : QgsDatumTransform::TransformDetails QgsDatumTransform::transformDetailsFromPj( PJ *op )
308 : : {
309 : 0 : PJ_CONTEXT *pjContext = QgsProjContext::get();
310 : 0 : TransformDetails details;
311 : 0 : if ( !op )
312 : 0 : return details;
313 : :
314 : 0 : QgsProjUtils::proj_pj_unique_ptr normalized( proj_normalize_for_visualization( pjContext, op ) );
315 : 0 : if ( normalized )
316 : 0 : details.proj = QString( proj_as_proj_string( pjContext, normalized.get(), PJ_PROJ_5, nullptr ) );
317 : :
318 : 0 : if ( details.proj.isEmpty() )
319 : 0 : details.proj = QString( proj_as_proj_string( pjContext, op, PJ_PROJ_5, nullptr ) );
320 : :
321 : 0 : if ( details.proj.isEmpty() )
322 : 0 : return details;
323 : :
324 : 0 : details.name = QString( proj_get_name( op ) );
325 : 0 : details.accuracy = proj_coordoperation_get_accuracy( pjContext, op );
326 : 0 : details.isAvailable = proj_coordoperation_is_instantiable( pjContext, op );
327 : :
328 : 0 : details.authority = QString( proj_get_id_auth_name( op, 0 ) );
329 : 0 : details.code = QString( proj_get_id_code( op, 0 ) );
330 : :
331 : 0 : const char *areaOfUseName = nullptr;
332 : 0 : double westLon = 0;
333 : 0 : double southLat = 0;
334 : 0 : double eastLon = 0;
335 : 0 : double northLat = 0;
336 : 0 : if ( proj_get_area_of_use( pjContext, op, &westLon, &southLat, &eastLon, &northLat, &areaOfUseName ) )
337 : : {
338 : 0 : details.areaOfUse = QString( areaOfUseName );
339 : : // don't use the constructor which normalizes!
340 : 0 : details.bounds.setXMinimum( westLon );
341 : 0 : details.bounds.setYMinimum( southLat );
342 : 0 : details.bounds.setXMaximum( eastLon );
343 : 0 : details.bounds.setYMaximum( northLat );
344 : 0 : }
345 : :
346 : 0 : details.remarks = QString( proj_get_remarks( op ) );
347 : 0 : details.scope = QString( proj_get_scope( op ) );
348 : :
349 : 0 : for ( int j = 0; j < proj_coordoperation_get_grid_used_count( pjContext, op ); ++j )
350 : : {
351 : 0 : const char *shortName = nullptr;
352 : 0 : const char *fullName = nullptr;
353 : 0 : const char *packageName = nullptr;
354 : 0 : const char *url = nullptr;
355 : 0 : int directDownload = 0;
356 : 0 : int openLicense = 0;
357 : 0 : int isAvailable = 0;
358 : 0 : proj_coordoperation_get_grid_used( pjContext, op, j, &shortName, &fullName, &packageName, &url, &directDownload, &openLicense, &isAvailable );
359 : 0 : GridDetails gridDetails;
360 : 0 : gridDetails.shortName = QString( shortName );
361 : 0 : gridDetails.fullName = QString( fullName );
362 : 0 : gridDetails.packageName = QString( packageName );
363 : 0 : gridDetails.url = QString( url );
364 : 0 : gridDetails.directDownload = directDownload;
365 : 0 : gridDetails.openLicense = openLicense;
366 : 0 : gridDetails.isAvailable = isAvailable;
367 : :
368 : 0 : details.grids.append( gridDetails );
369 : 0 : }
370 : :
371 : 0 : if ( proj_get_type( op ) == PJ_TYPE_CONCATENATED_OPERATION )
372 : : {
373 : 0 : for ( int j = 0; j < proj_concatoperation_get_step_count( pjContext, op ); ++j )
374 : : {
375 : 0 : QgsProjUtils::proj_pj_unique_ptr step( proj_concatoperation_get_step( pjContext, op, j ) );
376 : 0 : if ( step )
377 : : {
378 : 0 : SingleOperationDetails singleOpDetails;
379 : 0 : singleOpDetails.remarks = QString( proj_get_remarks( step.get() ) );
380 : 0 : singleOpDetails.scope = QString( proj_get_scope( step.get() ) );
381 : 0 : singleOpDetails.authority = QString( proj_get_id_auth_name( step.get(), 0 ) );
382 : 0 : singleOpDetails.code = QString( proj_get_id_code( step.get(), 0 ) );
383 : :
384 : 0 : const char *areaOfUseName = nullptr;
385 : 0 : if ( proj_get_area_of_use( pjContext, step.get(), nullptr, nullptr, nullptr, nullptr, &areaOfUseName ) )
386 : : {
387 : 0 : singleOpDetails.areaOfUse = QString( areaOfUseName );
388 : 0 : }
389 : 0 : details.operationDetails.append( singleOpDetails );
390 : 0 : }
391 : 0 : }
392 : 0 : }
393 : : else
394 : : {
395 : 0 : SingleOperationDetails singleOpDetails;
396 : 0 : singleOpDetails.remarks = QString( proj_get_remarks( op ) );
397 : 0 : singleOpDetails.scope = QString( proj_get_scope( op ) );
398 : 0 : singleOpDetails.authority = QString( proj_get_id_auth_name( op, 0 ) );
399 : 0 : singleOpDetails.code = QString( proj_get_id_code( op, 0 ) );
400 : :
401 : 0 : const char *areaOfUseName = nullptr;
402 : 0 : if ( proj_get_area_of_use( pjContext, op, nullptr, nullptr, nullptr, nullptr, &areaOfUseName ) )
403 : : {
404 : 0 : singleOpDetails.areaOfUse = QString( areaOfUseName );
405 : 0 : }
406 : 0 : details.operationDetails.append( singleOpDetails );
407 : 0 : }
408 : :
409 : 0 : return details;
410 : 0 : }
|
