Branch data Line data Source code
1 : : /***************************************************************************
2 : : QgsCoordinateTransform.cpp - Coordinate Transforms
3 : : -------------------
4 : : begin : Dec 2004
5 : : copyright : (C) 2004 Tim Sutton
6 : : email : tim at linfiniti.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 "qgscoordinatetransform.h"
18 : : #include "qgscoordinatetransform_p.h"
19 : : #include "qgsapplication.h"
20 : : #include "qgsmessagelog.h"
21 : : #include "qgslogger.h"
22 : : #include "qgspointxy.h"
23 : : #include "qgsrectangle.h"
24 : : #include "qgsexception.h"
25 : : #include "qgsproject.h"
26 : : #include "qgsreadwritelocker.h"
27 : : #include "qgsvector3d.h"
28 : :
29 : : //qt includes
30 : : #include <QDomNode>
31 : : #include <QDomElement>
32 : : #include <QApplication>
33 : : #include <QPolygonF>
34 : : #include <QStringList>
35 : : #include <QVector>
36 : :
37 : : #include <proj.h>
38 : : #include "qgsprojutils.h"
39 : :
40 : : #include <sqlite3.h>
41 : : #include <qlogging.h>
42 : : #include <vector>
43 : : #include <algorithm>
44 : :
45 : : // if defined shows all information about transform to stdout
46 : : // #define COORDINATE_TRANSFORM_VERBOSE
47 : :
48 : 5 : QReadWriteLock QgsCoordinateTransform::sCacheLock;
49 : 5 : QMultiHash< QPair< QString, QString >, QgsCoordinateTransform > QgsCoordinateTransform::sTransforms; //same auth_id pairs might have different datum transformations
50 : : bool QgsCoordinateTransform::sDisableCache = false;
51 : :
52 : : std::function< void( const QgsCoordinateReferenceSystem &sourceCrs,
53 : : const QgsCoordinateReferenceSystem &destinationCrs,
54 : 5 : const QString &desiredOperation )> QgsCoordinateTransform::sFallbackOperationOccurredHandler = nullptr;
55 : :
56 : 318 : QgsCoordinateTransform::QgsCoordinateTransform()
57 : : {
58 : 318 : d = new QgsCoordinateTransformPrivate();
59 : 318 : }
60 : :
61 : 1523 : QgsCoordinateTransform::QgsCoordinateTransform( const QgsCoordinateReferenceSystem &source, const QgsCoordinateReferenceSystem &destination, const QgsCoordinateTransformContext &context )
62 : : {
63 : 1523 : mContext = context;
64 : 1523 : d = new QgsCoordinateTransformPrivate( source, destination, mContext );
65 : : #ifdef QGISDEBUG
66 : : mHasContext = true;
67 : : #endif
68 : :
69 : 1523 : if ( !d->checkValidity() )
70 : 1517 : return;
71 : :
72 : : Q_NOWARN_DEPRECATED_PUSH
73 : 6 : if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation, d->mAllowFallbackTransforms ) )
74 : : {
75 : 1 : d->initialize();
76 : 1 : addToCache();
77 : 1 : }
78 : : Q_NOWARN_DEPRECATED_POP
79 : 1523 : }
80 : :
81 : 10 : QgsCoordinateTransform::QgsCoordinateTransform( const QgsCoordinateReferenceSystem &source, const QgsCoordinateReferenceSystem &destination, const QgsProject *project )
82 : : {
83 : 10 : mContext = project ? project->transformContext() : QgsCoordinateTransformContext();
84 : 10 : d = new QgsCoordinateTransformPrivate( source, destination, mContext );
85 : : #ifdef QGISDEBUG
86 : : if ( project )
87 : : mHasContext = true;
88 : : #endif
89 : :
90 : 10 : if ( !d->checkValidity() )
91 : 4 : return;
92 : :
93 : : Q_NOWARN_DEPRECATED_PUSH
94 : 6 : if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation, d->mAllowFallbackTransforms ) )
95 : : {
96 : 1 : d->initialize();
97 : 1 : addToCache();
98 : 1 : }
99 : : Q_NOWARN_DEPRECATED_POP
100 : 10 : }
101 : :
102 : 0 : QgsCoordinateTransform::QgsCoordinateTransform( const QgsCoordinateReferenceSystem &source, const QgsCoordinateReferenceSystem &destination, int sourceDatumTransform, int destinationDatumTransform )
103 : : {
104 : 0 : d = new QgsCoordinateTransformPrivate( source, destination, sourceDatumTransform, destinationDatumTransform );
105 : : #ifdef QGISDEBUG
106 : : mHasContext = true; // not strictly true, but we don't need to worry if datums have been explicitly set
107 : : #endif
108 : :
109 : 0 : if ( !d->checkValidity() )
110 : 0 : return;
111 : :
112 : : Q_NOWARN_DEPRECATED_PUSH
113 : 0 : if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation, d->mAllowFallbackTransforms ) )
114 : : {
115 : 0 : d->initialize();
116 : 0 : addToCache();
117 : 0 : }
118 : : Q_NOWARN_DEPRECATED_POP
119 : 0 : }
120 : :
121 : 12 : QgsCoordinateTransform::QgsCoordinateTransform( const QgsCoordinateTransform &o )
122 : 12 : : mContext( o.mContext )
123 : : #ifdef QGISDEBUG
124 : : , mHasContext( o.mHasContext )
125 : : #endif
126 : 12 : , mLastError()
127 : : {
128 : 12 : d = o.d;
129 : 12 : }
130 : :
131 : 10 : QgsCoordinateTransform &QgsCoordinateTransform::operator=( const QgsCoordinateTransform &o ) //NOLINT
132 : : {
133 : 10 : d = o.d;
134 : : #ifdef QGISDEBUG
135 : : mHasContext = o.mHasContext;
136 : : #endif
137 : 10 : mContext = o.mContext;
138 : 10 : mLastError = QString();
139 : 10 : return *this;
140 : : }
141 : :
142 : 1863 : QgsCoordinateTransform::~QgsCoordinateTransform() {} //NOLINT
143 : :
144 : 0 : void QgsCoordinateTransform::setSourceCrs( const QgsCoordinateReferenceSystem &crs )
145 : : {
146 : 0 : d.detach();
147 : 0 : d->mSourceCRS = crs;
148 : 0 : if ( !d->checkValidity() )
149 : 0 : return;
150 : :
151 : 0 : d->calculateTransforms( mContext );
152 : : Q_NOWARN_DEPRECATED_PUSH
153 : 0 : if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation, d->mAllowFallbackTransforms ) )
154 : : {
155 : 0 : d->initialize();
156 : 0 : addToCache();
157 : 0 : }
158 : : Q_NOWARN_DEPRECATED_POP
159 : 0 : }
160 : 0 : void QgsCoordinateTransform::setDestinationCrs( const QgsCoordinateReferenceSystem &crs )
161 : : {
162 : 0 : d.detach();
163 : 0 : d->mDestCRS = crs;
164 : 0 : if ( !d->checkValidity() )
165 : 0 : return;
166 : :
167 : 0 : d->calculateTransforms( mContext );
168 : : Q_NOWARN_DEPRECATED_PUSH
169 : 0 : if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation, d->mAllowFallbackTransforms ) )
170 : : {
171 : 0 : d->initialize();
172 : 0 : addToCache();
173 : 0 : }
174 : : Q_NOWARN_DEPRECATED_POP
175 : 0 : }
176 : :
177 : 0 : void QgsCoordinateTransform::setContext( const QgsCoordinateTransformContext &context )
178 : : {
179 : 0 : d.detach();
180 : 0 : mContext = context;
181 : : #ifdef QGISDEBUG
182 : : mHasContext = true;
183 : : #endif
184 : 0 : if ( !d->checkValidity() )
185 : 0 : return;
186 : :
187 : 0 : d->calculateTransforms( mContext );
188 : : Q_NOWARN_DEPRECATED_PUSH
189 : 0 : if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation, d->mAllowFallbackTransforms ) )
190 : : {
191 : 0 : d->initialize();
192 : 0 : addToCache();
193 : 0 : }
194 : : Q_NOWARN_DEPRECATED_POP
195 : 0 : }
196 : :
197 : 0 : QgsCoordinateTransformContext QgsCoordinateTransform::context() const
198 : : {
199 : 0 : return mContext;
200 : : }
201 : :
202 : 0 : QgsCoordinateReferenceSystem QgsCoordinateTransform::sourceCrs() const
203 : : {
204 : 0 : return d->mSourceCRS;
205 : : }
206 : :
207 : 0 : QgsCoordinateReferenceSystem QgsCoordinateTransform::destinationCrs() const
208 : : {
209 : 0 : return d->mDestCRS;
210 : : }
211 : :
212 : 252 : QgsPointXY QgsCoordinateTransform::transform( const QgsPointXY &point, TransformDirection direction ) const
213 : : {
214 : 252 : if ( !d->mIsValid || d->mShortCircuit )
215 : 252 : return point;
216 : :
217 : : // transform x
218 : 0 : double x = point.x();
219 : 0 : double y = point.y();
220 : 0 : double z = 0.0;
221 : : try
222 : : {
223 : 0 : transformCoords( 1, &x, &y, &z, direction );
224 : 0 : }
225 : : catch ( const QgsCsException & )
226 : : {
227 : : // rethrow the exception
228 : 0 : QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
229 : 0 : throw;
230 : 0 : }
231 : :
232 : 0 : return QgsPointXY( x, y );
233 : 252 : }
234 : :
235 : :
236 : 0 : QgsPointXY QgsCoordinateTransform::transform( const double theX, const double theY = 0.0, TransformDirection direction ) const
237 : : {
238 : : try
239 : : {
240 : 0 : return transform( QgsPointXY( theX, theY ), direction );
241 : 0 : }
242 : : catch ( const QgsCsException & )
243 : : {
244 : : // rethrow the exception
245 : 0 : QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
246 : 0 : throw;
247 : 0 : }
248 : 0 : }
249 : :
250 : 267 : QgsRectangle QgsCoordinateTransform::transform( const QgsRectangle &rect, TransformDirection direction ) const
251 : : {
252 : 267 : if ( !d->mIsValid || d->mShortCircuit )
253 : 267 : return rect;
254 : : // transform x
255 : 0 : double x1 = rect.xMinimum();
256 : 0 : double y1 = rect.yMinimum();
257 : 0 : double x2 = rect.xMaximum();
258 : 0 : double y2 = rect.yMaximum();
259 : :
260 : : // Number of points to reproject------+
261 : : // |
262 : : // V
263 : : try
264 : : {
265 : 0 : double z = 0.0;
266 : 0 : transformCoords( 1, &x1, &y1, &z, direction );
267 : 0 : transformCoords( 1, &x2, &y2, &z, direction );
268 : 0 : }
269 : : catch ( const QgsCsException & )
270 : : {
271 : : // rethrow the exception
272 : 0 : QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
273 : 0 : throw;
274 : 0 : }
275 : :
276 : : #ifdef COORDINATE_TRANSFORM_VERBOSE
277 : : QgsDebugMsg( QStringLiteral( "Rect projection..." ) );
278 : : QgsDebugMsg( QStringLiteral( "Xmin : %1 --> %2" ).arg( rect.xMinimum() ).arg( x1 ) );
279 : : QgsDebugMsg( QStringLiteral( "Ymin : %1 --> %2" ).arg( rect.yMinimum() ).arg( y1 ) );
280 : : QgsDebugMsg( QStringLiteral( "Xmax : %1 --> %2" ).arg( rect.xMaximum() ).arg( x2 ) );
281 : : QgsDebugMsg( QStringLiteral( "Ymax : %1 --> %2" ).arg( rect.yMaximum() ).arg( y2 ) );
282 : : #endif
283 : 0 : return QgsRectangle( x1, y1, x2, y2 );
284 : 267 : }
285 : :
286 : 0 : QgsVector3D QgsCoordinateTransform::transform( const QgsVector3D &point, TransformDirection direction ) const
287 : : {
288 : 0 : double x = point.x();
289 : 0 : double y = point.y();
290 : 0 : double z = point.z();
291 : : try
292 : : {
293 : 0 : transformCoords( 1, &x, &y, &z, direction );
294 : 0 : }
295 : : catch ( const QgsCsException & )
296 : : {
297 : : // rethrow the exception
298 : 0 : QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
299 : 0 : throw;
300 : 0 : }
301 : 0 : return QgsVector3D( x, y, z );
302 : 0 : }
303 : :
304 : 2 : void QgsCoordinateTransform::transformInPlace( double &x, double &y, double &z,
305 : : TransformDirection direction ) const
306 : : {
307 : 2 : if ( !d->mIsValid || d->mShortCircuit )
308 : 0 : return;
309 : : #ifdef QGISDEBUG
310 : : // QgsDebugMsg(QString("Using transform in place %1 %2").arg(__FILE__).arg(__LINE__));
311 : : #endif
312 : : // transform x
313 : : try
314 : : {
315 : 2 : transformCoords( 1, &x, &y, &z, direction );
316 : 2 : }
317 : : catch ( const QgsCsException & )
318 : : {
319 : : // rethrow the exception
320 : 0 : QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
321 : 0 : throw;
322 : 0 : }
323 : 2 : }
324 : :
325 : 0 : void QgsCoordinateTransform::transformInPlace( float &x, float &y, double &z,
326 : : TransformDirection direction ) const
327 : : {
328 : 0 : double xd = static_cast< double >( x ), yd = static_cast< double >( y );
329 : 0 : transformInPlace( xd, yd, z, direction );
330 : 0 : x = xd;
331 : 0 : y = yd;
332 : 0 : }
333 : :
334 : 0 : void QgsCoordinateTransform::transformInPlace( float &x, float &y, float &z,
335 : : TransformDirection direction ) const
336 : : {
337 : 0 : if ( !d->mIsValid || d->mShortCircuit )
338 : 0 : return;
339 : : #ifdef QGISDEBUG
340 : : // QgsDebugMsg(QString("Using transform in place %1 %2").arg(__FILE__).arg(__LINE__));
341 : : #endif
342 : : // transform x
343 : : try
344 : : {
345 : 0 : double xd = x;
346 : 0 : double yd = y;
347 : 0 : double zd = z;
348 : 0 : transformCoords( 1, &xd, &yd, &zd, direction );
349 : 0 : x = xd;
350 : 0 : y = yd;
351 : 0 : z = zd;
352 : 0 : }
353 : : catch ( QgsCsException & )
354 : : {
355 : : // rethrow the exception
356 : 0 : QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
357 : 0 : throw;
358 : 0 : }
359 : 0 : }
360 : :
361 : 0 : void QgsCoordinateTransform::transformPolygon( QPolygonF &poly, TransformDirection direction ) const
362 : : {
363 : 0 : if ( !d->mIsValid || d->mShortCircuit )
364 : : {
365 : 0 : return;
366 : : }
367 : :
368 : : //create x, y arrays
369 : 0 : int nVertices = poly.size();
370 : :
371 : 0 : QVector<double> x( nVertices );
372 : 0 : QVector<double> y( nVertices );
373 : 0 : QVector<double> z( nVertices );
374 : 0 : double *destX = x.data();
375 : 0 : double *destY = y.data();
376 : 0 : double *destZ = z.data();
377 : :
378 : 0 : const QPointF *polyData = poly.constData();
379 : 0 : for ( int i = 0; i < nVertices; ++i )
380 : : {
381 : 0 : *destX++ = polyData->x();
382 : 0 : *destY++ = polyData->y();
383 : 0 : *destZ++ = 0;
384 : 0 : polyData++;
385 : 0 : }
386 : :
387 : 0 : QString err;
388 : : try
389 : : {
390 : 0 : transformCoords( nVertices, x.data(), y.data(), z.data(), direction );
391 : 0 : }
392 : : catch ( const QgsCsException &e )
393 : : {
394 : : // record the exception, but don't rethrow it until we've recorded the coordinates we *could* transform
395 : 0 : err = e.what();
396 : 0 : }
397 : :
398 : 0 : QPointF *destPoint = poly.data();
399 : 0 : const double *srcX = x.constData();
400 : 0 : const double *srcY = y.constData();
401 : 0 : for ( int i = 0; i < nVertices; ++i )
402 : : {
403 : 0 : destPoint->rx() = *srcX++;
404 : 0 : destPoint->ry() = *srcY++;
405 : 0 : destPoint++;
406 : 0 : }
407 : :
408 : : // rethrow the exception
409 : 0 : if ( !err.isEmpty() )
410 : 0 : throw QgsCsException( err );
411 : 0 : }
412 : :
413 : 0 : void QgsCoordinateTransform::transformInPlace(
414 : : QVector<double> &x, QVector<double> &y, QVector<double> &z,
415 : : TransformDirection direction ) const
416 : : {
417 : :
418 : 0 : if ( !d->mIsValid || d->mShortCircuit )
419 : 0 : return;
420 : :
421 : : Q_ASSERT( x.size() == y.size() );
422 : :
423 : : // Apparently, if one has a std::vector, it is valid to use the
424 : : // address of the first element in the vector as a pointer to an
425 : : // array of the vectors data, and hence easily interface with code
426 : : // that wants C-style arrays.
427 : :
428 : : try
429 : : {
430 : 0 : transformCoords( x.size(), &x[0], &y[0], &z[0], direction );
431 : 0 : }
432 : : catch ( const QgsCsException & )
433 : : {
434 : : // rethrow the exception
435 : 0 : QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
436 : 0 : throw;
437 : 0 : }
438 : 0 : }
439 : :
440 : :
441 : 0 : void QgsCoordinateTransform::transformInPlace(
442 : : QVector<float> &x, QVector<float> &y, QVector<float> &z,
443 : : TransformDirection direction ) const
444 : : {
445 : 0 : if ( !d->mIsValid || d->mShortCircuit )
446 : 0 : return;
447 : :
448 : : Q_ASSERT( x.size() == y.size() );
449 : :
450 : : // Apparently, if one has a std::vector, it is valid to use the
451 : : // address of the first element in the vector as a pointer to an
452 : : // array of the vectors data, and hence easily interface with code
453 : : // that wants C-style arrays.
454 : :
455 : : try
456 : : {
457 : : //copy everything to double vectors since proj needs double
458 : 0 : int vectorSize = x.size();
459 : 0 : QVector<double> xd( x.size() );
460 : 0 : QVector<double> yd( y.size() );
461 : 0 : QVector<double> zd( z.size() );
462 : :
463 : 0 : double *destX = xd.data();
464 : 0 : double *destY = yd.data();
465 : 0 : double *destZ = zd.data();
466 : :
467 : 0 : const float *srcX = x.constData();
468 : 0 : const float *srcY = y.constData();
469 : 0 : const float *srcZ = z.constData();
470 : :
471 : 0 : for ( int i = 0; i < vectorSize; ++i )
472 : : {
473 : 0 : *destX++ = static_cast< double >( *srcX++ );
474 : 0 : *destY++ = static_cast< double >( *srcY++ );
475 : 0 : *destZ++ = static_cast< double >( *srcZ++ );
476 : 0 : }
477 : :
478 : 0 : transformCoords( x.size(), &xd[0], &yd[0], &zd[0], direction );
479 : :
480 : : //copy back
481 : 0 : float *destFX = x.data();
482 : 0 : float *destFY = y.data();
483 : 0 : float *destFZ = z.data();
484 : 0 : const double *srcXD = xd.constData();
485 : 0 : const double *srcYD = yd.constData();
486 : 0 : const double *srcZD = zd.constData();
487 : 0 : for ( int i = 0; i < vectorSize; ++i )
488 : : {
489 : 0 : *destFX++ = static_cast< float >( *srcXD++ );
490 : 0 : *destFY++ = static_cast< float >( *srcYD++ );
491 : 0 : *destFZ++ = static_cast< float >( *srcZD++ );
492 : 0 : }
493 : 0 : }
494 : : catch ( QgsCsException & )
495 : : {
496 : : // rethrow the exception
497 : 0 : QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
498 : 0 : throw;
499 : 0 : }
500 : 0 : }
501 : :
502 : 4 : QgsRectangle QgsCoordinateTransform::transformBoundingBox( const QgsRectangle &rect, TransformDirection direction, const bool handle180Crossover ) const
503 : : {
504 : : // Calculate the bounding box of a QgsRectangle in the source CRS
505 : : // when projected to the destination CRS (or the inverse).
506 : : // This is done by looking at a number of points spread evenly
507 : : // across the rectangle
508 : :
509 : 4 : if ( !d->mIsValid || d->mShortCircuit )
510 : 4 : return rect;
511 : :
512 : 0 : if ( rect.isEmpty() )
513 : : {
514 : 0 : QgsPointXY p = transform( rect.xMinimum(), rect.yMinimum(), direction );
515 : 0 : return QgsRectangle( p, p );
516 : : }
517 : :
518 : : // 64 points (<=2.12) is not enough, see #13665, for EPSG:4326 -> EPSG:3574 (say that it is a hard one),
519 : : // are decent result from about 500 points and more. This method is called quite often, but
520 : : // even with 1000 points it takes < 1ms.
521 : : // TODO: how to effectively and precisely reproject bounding box?
522 : 0 : const int nPoints = 1000;
523 : 0 : double d = std::sqrt( ( rect.width() * rect.height() ) / std::pow( std::sqrt( static_cast< double >( nPoints ) ) - 1, 2.0 ) );
524 : 0 : int nXPoints = std::min( static_cast< int >( std::ceil( rect.width() / d ) ) + 1, 1000 );
525 : 0 : int nYPoints = std::min( static_cast< int >( std::ceil( rect.height() / d ) ) + 1, 1000 );
526 : :
527 : 0 : QgsRectangle bb_rect;
528 : 0 : bb_rect.setMinimal();
529 : :
530 : : // We're interfacing with C-style vectors in the
531 : : // end, so let's do C-style vectors here too.
532 : 0 : QVector<double> x( nXPoints * nYPoints );
533 : 0 : QVector<double> y( nXPoints * nYPoints );
534 : 0 : QVector<double> z( nXPoints * nYPoints );
535 : :
536 : 0 : QgsDebugMsgLevel( QStringLiteral( "Entering transformBoundingBox..." ), 4 );
537 : :
538 : : // Populate the vectors
539 : :
540 : 0 : double dx = rect.width() / static_cast< double >( nXPoints - 1 );
541 : 0 : double dy = rect.height() / static_cast< double >( nYPoints - 1 );
542 : :
543 : 0 : double pointY = rect.yMinimum();
544 : :
545 : 0 : for ( int i = 0; i < nYPoints ; i++ )
546 : : {
547 : :
548 : : // Start at right edge
549 : 0 : double pointX = rect.xMinimum();
550 : :
551 : 0 : for ( int j = 0; j < nXPoints; j++ )
552 : : {
553 : 0 : x[( i * nXPoints ) + j] = pointX;
554 : 0 : y[( i * nXPoints ) + j] = pointY;
555 : : // and the height...
556 : 0 : z[( i * nXPoints ) + j] = 0.0;
557 : : // QgsDebugMsg(QString("BBox coord: (%1, %2)").arg(x[(i*numP) + j]).arg(y[(i*numP) + j]));
558 : 0 : pointX += dx;
559 : 0 : }
560 : 0 : pointY += dy;
561 : 0 : }
562 : :
563 : : // Do transformation. Any exception generated must
564 : : // be handled in above layers.
565 : : try
566 : : {
567 : 0 : transformCoords( nXPoints * nYPoints, x.data(), y.data(), z.data(), direction );
568 : 0 : }
569 : : catch ( const QgsCsException & )
570 : : {
571 : : // rethrow the exception
572 : 0 : QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
573 : 0 : throw;
574 : 0 : }
575 : :
576 : : // Calculate the bounding box and use that for the extent
577 : :
578 : 0 : for ( int i = 0; i < nXPoints * nYPoints; i++ )
579 : : {
580 : 0 : if ( !std::isfinite( x[i] ) || !std::isfinite( y[i] ) )
581 : : {
582 : 0 : continue;
583 : : }
584 : :
585 : 0 : if ( handle180Crossover )
586 : : {
587 : : //if crossing the date line, temporarily add 360 degrees to -ve longitudes
588 : 0 : bb_rect.combineExtentWith( x[i] >= 0.0 ? x[i] : x[i] + 360.0, y[i] );
589 : 0 : }
590 : : else
591 : : {
592 : 0 : bb_rect.combineExtentWith( x[i], y[i] );
593 : : }
594 : 0 : }
595 : :
596 : 0 : if ( bb_rect.isNull() )
597 : : {
598 : : // something bad happened when reprojecting the filter rect... no finite points were left!
599 : 0 : throw QgsCsException( QObject::tr( "Could not transform bounding box to target CRS" ) );
600 : : }
601 : :
602 : 0 : if ( handle180Crossover )
603 : : {
604 : : //subtract temporary addition of 360 degrees from longitudes
605 : 0 : if ( bb_rect.xMinimum() > 180.0 )
606 : 0 : bb_rect.setXMinimum( bb_rect.xMinimum() - 360.0 );
607 : 0 : if ( bb_rect.xMaximum() > 180.0 )
608 : 0 : bb_rect.setXMaximum( bb_rect.xMaximum() - 360.0 );
609 : 0 : }
610 : :
611 : 0 : QgsDebugMsgLevel( "Projected extent: " + bb_rect.toString(), 4 );
612 : :
613 : 0 : if ( bb_rect.isEmpty() )
614 : : {
615 : 0 : QgsDebugMsgLevel( "Original extent: " + rect.toString(), 4 );
616 : 0 : }
617 : :
618 : 0 : return bb_rect;
619 : 4 : }
620 : :
621 : 122 : void QgsCoordinateTransform::transformCoords( int numPoints, double *x, double *y, double *z, TransformDirection direction ) const
622 : : {
623 : 122 : if ( !d->mIsValid || d->mShortCircuit )
624 : 96 : return;
625 : : // Refuse to transform the points if the srs's are invalid
626 : 26 : if ( !d->mSourceCRS.isValid() )
627 : : {
628 : 0 : QgsMessageLog::logMessage( QObject::tr( "The source spatial reference system (CRS) is not valid. "
629 : : "The coordinates can not be reprojected. The CRS is: %1" )
630 : 0 : .arg( d->mSourceCRS.toProj() ), QObject::tr( "CRS" ) );
631 : 0 : return;
632 : : }
633 : 26 : if ( !d->mDestCRS.isValid() )
634 : : {
635 : 0 : QgsMessageLog::logMessage( QObject::tr( "The destination spatial reference system (CRS) is not valid. "
636 : 0 : "The coordinates can not be reprojected. The CRS is: %1" ).arg( d->mDestCRS.toProj() ), QObject::tr( "CRS" ) );
637 : 0 : return;
638 : : }
639 : :
640 : 26 : std::vector< int > zNanPositions;
641 : 100 : for ( int i = 0; i < numPoints; i++ )
642 : : {
643 : 74 : if ( std::isnan( z[i] ) )
644 : : {
645 : 0 : zNanPositions.push_back( i );
646 : 0 : z[i] = 0.0;
647 : 0 : }
648 : 74 : }
649 : :
650 : 26 : std::vector< double > xprev( numPoints );
651 : 26 : memcpy( xprev.data(), x, sizeof( double ) * numPoints );
652 : 26 : std::vector< double > yprev( numPoints );
653 : 26 : memcpy( yprev.data(), y, sizeof( double ) * numPoints );
654 : 26 : std::vector< double > zprev( numPoints );
655 : 26 : memcpy( zprev.data(), z, sizeof( double ) * numPoints );
656 : :
657 : : #ifdef COORDINATE_TRANSFORM_VERBOSE
658 : : double xorg = *x;
659 : : double yorg = *y;
660 : : QgsDebugMsg( QStringLiteral( "[[[[[[ Number of points to transform: %1 ]]]]]]" ).arg( numPoints ) );
661 : : #endif
662 : :
663 : : #ifdef QGISDEBUG
664 : : if ( !mHasContext )
665 : : QgsDebugMsgLevel( QStringLiteral( "No QgsCoordinateTransformContext context set for transform" ), 4 );
666 : : #endif
667 : :
668 : : // use proj4 to do the transform
669 : :
670 : : // if the source/destination projection is lat/long, convert the points to radians
671 : : // prior to transforming
672 : 26 : ProjData projData = d->threadLocalProjData();
673 : :
674 : 26 : int projResult = 0;
675 : :
676 : 26 : proj_errno_reset( projData );
677 : 26 : proj_trans_generic( projData, ( direction == ForwardTransform && !d->mIsReversed ) || ( direction == ReverseTransform && d->mIsReversed ) ? PJ_FWD : PJ_INV,
678 : 26 : x, sizeof( double ), numPoints,
679 : 26 : y, sizeof( double ), numPoints,
680 : 26 : z, sizeof( double ), numPoints,
681 : : nullptr, sizeof( double ), 0 );
682 : : // Try to - approximatively - emulate the behavior of pj_transform()...
683 : : // In the case of a single point transform, and a transformation error occurs,
684 : : // pj_transform() would return the errno. In cases of multiple point transform,
685 : : // it would continue (for non-transient errors, that is pipeline definition
686 : : // errors) and just set the resulting x,y to infinity. This is in fact a
687 : : // bit more subtle than that, and I'm not completely sure the logic in
688 : : // pj_transform() was really sane & fully bullet proof
689 : : // So here just check proj_errno() for single point transform
690 : 26 : int actualRes = 0;
691 : 1889 : if ( numPoints == 1 )
692 : 1863 : {
693 : 1865 : projResult = proj_errno( projData );
694 : 2 : actualRes = projResult;
695 : 2 : }
696 : : else
697 : : {
698 : 24 : actualRes = proj_errno( projData );
699 : : }
700 : 26 : if ( actualRes == 0 )
701 : : {
702 : : // proj_errno is sometimes not an accurate method to test for transform failures - so we need to
703 : : // manually scan for nan values
704 : 126 : if ( std::any_of( x, x + numPoints, []( double v ) { return std::isinf( v ); } )
705 : 100 : || std::any_of( y, y + numPoints, []( double v ) { return std::isinf( v ); } )
706 : 100 : || std::any_of( z, z + numPoints, []( double v ) { return std::isinf( v ); } ) )
707 : : {
708 : 0 : actualRes = 1;
709 : 0 : }
710 : 26 : }
711 : :
712 : 26 : mFallbackOperationOccurred = false;
713 : 26 : if ( actualRes != 0
714 : 26 : && ( d->mAvailableOpCount > 1 || d->mAvailableOpCount == -1 ) // only use fallbacks if more than one operation is possible -- otherwise we've already tried it and it failed
715 : 0 : && ( d->mAllowFallbackTransforms || mBallparkTransformsAreAppropriate ) )
716 : : {
717 : : // fail #1 -- try with getting proj to auto-pick an appropriate coordinate operation for the points
718 : 0 : if ( PJ *transform = d->threadLocalFallbackProjData() )
719 : : {
720 : 0 : projResult = 0;
721 : 0 : proj_errno_reset( transform );
722 : 0 : proj_trans_generic( transform, direction == ForwardTransform ? PJ_FWD : PJ_INV,
723 : 0 : xprev.data(), sizeof( double ), numPoints,
724 : 0 : yprev.data(), sizeof( double ), numPoints,
725 : 0 : zprev.data(), sizeof( double ), numPoints,
726 : : nullptr, sizeof( double ), 0 );
727 : : // Try to - approximatively - emulate the behavior of pj_transform()...
728 : : // In the case of a single point transform, and a transformation error occurs,
729 : : // pj_transform() would return the errno. In cases of multiple point transform,
730 : : // it would continue (for non-transient errors, that is pipeline definition
731 : : // errors) and just set the resulting x,y to infinity. This is in fact a
732 : : // bit more subtle than that, and I'm not completely sure the logic in
733 : : // pj_transform() was really sane & fully bullet proof
734 : : // So here just check proj_errno() for single point transform
735 : 0 : if ( numPoints == 1 )
736 : : {
737 : : // hmm - something very odd here. We can't trust proj_errno( transform ), as that's giving us incorrect error numbers
738 : : // (such as "failed to load datum shift file", which is definitely incorrect for a default proj created operation!)
739 : : // so we resort to testing values ourselves...
740 : 0 : projResult = std::isinf( xprev[0] ) || std::isinf( yprev[0] ) || std::isinf( zprev[0] ) ? 1 : 0;
741 : 0 : }
742 : :
743 : 0 : if ( projResult == 0 )
744 : : {
745 : 0 : memcpy( x, xprev.data(), sizeof( double ) * numPoints );
746 : 0 : memcpy( y, yprev.data(), sizeof( double ) * numPoints );
747 : 0 : memcpy( z, zprev.data(), sizeof( double ) * numPoints );
748 : 0 : mFallbackOperationOccurred = true;
749 : 0 : }
750 : :
751 : 0 : if ( !mBallparkTransformsAreAppropriate && !mDisableFallbackHandler && sFallbackOperationOccurredHandler )
752 : : {
753 : 0 : sFallbackOperationOccurredHandler( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation );
754 : : #if 0
755 : : const QString warning = QStringLiteral( "A fallback coordinate operation was used between %1 and %2" ).arg( d->mSourceCRS.authid(),
756 : : d->mDestCRS.authid() );
757 : : qWarning( "%s", warning.toLatin1().constData() );
758 : : #endif
759 : 0 : }
760 : 0 : }
761 : 0 : }
762 : :
763 : 26 : for ( const int &pos : zNanPositions )
764 : : {
765 : 0 : z[pos] = std::numeric_limits<double>::quiet_NaN();
766 : : }
767 : :
768 : 26 : if ( projResult != 0 )
769 : : {
770 : : //something bad happened....
771 : 0 : QString points;
772 : :
773 : 0 : for ( int i = 0; i < numPoints; ++i )
774 : : {
775 : 0 : if ( direction == ForwardTransform )
776 : : {
777 : 0 : points += QStringLiteral( "(%1, %2)\n" ).arg( x[i], 0, 'f' ).arg( y[i], 0, 'f' );
778 : 0 : }
779 : : else
780 : : {
781 : 0 : points += QStringLiteral( "(%1, %2)\n" ).arg( x[i], 0, 'f' ).arg( y[i], 0, 'f' );
782 : : }
783 : 0 : }
784 : :
785 : 0 : QString dir = ( direction == ForwardTransform ) ? QObject::tr( "forward transform" ) : QObject::tr( "inverse transform" );
786 : :
787 : 0 : QString msg = QObject::tr( "%1 of\n"
788 : : "%2"
789 : : "Error: %3" )
790 : 0 : .arg( dir,
791 : : points,
792 : 0 : projResult < 0 ? QString::fromUtf8( proj_errno_string( projResult ) ) : QObject::tr( "Fallback transform failed" ) );
793 : :
794 : :
795 : : // don't flood console with thousands of duplicate transform error messages
796 : 0 : if ( msg != mLastError )
797 : : {
798 : 0 : QgsDebugMsg( "Projection failed emitting invalid transform signal: " + msg );
799 : 0 : mLastError = msg;
800 : 0 : }
801 : 0 : QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
802 : :
803 : 0 : throw QgsCsException( msg );
804 : 0 : }
805 : :
806 : : #ifdef COORDINATE_TRANSFORM_VERBOSE
807 : : QgsDebugMsg( QStringLiteral( "[[[[[[ Projected %1, %2 to %3, %4 ]]]]]]" )
808 : : .arg( xorg, 0, 'g', 15 ).arg( yorg, 0, 'g', 15 )
809 : : .arg( *x, 0, 'g', 15 ).arg( *y, 0, 'g', 15 ) );
810 : : #endif
811 : 122 : }
812 : :
813 : 2222 : bool QgsCoordinateTransform::isValid() const
814 : : {
815 : 2222 : return d->mIsValid;
816 : : }
817 : :
818 : 717 : bool QgsCoordinateTransform::isShortCircuited() const
819 : : {
820 : 717 : return !d->mIsValid || d->mShortCircuit;
821 : : }
822 : :
823 : 10 : QString QgsCoordinateTransform::coordinateOperation() const
824 : : {
825 : 10 : return d->mProjCoordinateOperation;
826 : : }
827 : :
828 : 0 : QgsDatumTransform::TransformDetails QgsCoordinateTransform::instantiatedCoordinateOperationDetails() const
829 : : {
830 : 0 : ProjData projData = d->threadLocalProjData();
831 : 0 : return QgsDatumTransform::transformDetailsFromPj( projData );
832 : : }
833 : :
834 : 0 : void QgsCoordinateTransform::setCoordinateOperation( const QString &operation ) const
835 : : {
836 : 0 : d.detach();
837 : 0 : d->mProjCoordinateOperation = operation;
838 : 0 : d->mShouldReverseCoordinateOperation = false;
839 : 0 : }
840 : :
841 : 0 : void QgsCoordinateTransform::setAllowFallbackTransforms( bool allowed )
842 : : {
843 : 0 : d.detach();
844 : 0 : d->mAllowFallbackTransforms = allowed;
845 : 0 : }
846 : :
847 : 10 : bool QgsCoordinateTransform::allowFallbackTransforms() const
848 : : {
849 : 10 : return d->mAllowFallbackTransforms;
850 : : }
851 : :
852 : 0 : void QgsCoordinateTransform::setBallparkTransformsAreAppropriate( bool appropriate )
853 : : {
854 : 0 : mBallparkTransformsAreAppropriate = appropriate;
855 : 0 : }
856 : :
857 : 0 : void QgsCoordinateTransform::disableFallbackOperationHandler( bool disabled )
858 : : {
859 : 0 : mDisableFallbackHandler = disabled;
860 : 0 : }
861 : :
862 : 0 : bool QgsCoordinateTransform::fallbackOperationOccurred() const
863 : : {
864 : 0 : return mFallbackOperationOccurred;
865 : : }
866 : :
867 : 0 : const char *finder( const char *name )
868 : : {
869 : 0 : QString proj;
870 : : #ifdef Q_OS_WIN
871 : : proj = QApplication::applicationDirPath()
872 : : + "/share/proj/" + QString( name );
873 : : #else
874 : : Q_UNUSED( name )
875 : : #endif
876 : 0 : return proj.toUtf8();
877 : 0 : }
878 : :
879 : 12 : bool QgsCoordinateTransform::setFromCache( const QgsCoordinateReferenceSystem &src, const QgsCoordinateReferenceSystem &dest, const QString &coordinateOperationProj, bool allowFallback )
880 : : {
881 : 12 : if ( !src.isValid() || !dest.isValid() )
882 : 0 : return false;
883 : :
884 : 24 : const QString sourceKey = src.authid().isEmpty() ?
885 : 12 : src.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED ) : src.authid();
886 : 24 : const QString destKey = dest.authid().isEmpty() ?
887 : 12 : dest.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED ) : dest.authid();
888 : :
889 : 12 : if ( sourceKey.isEmpty() || destKey.isEmpty() )
890 : 0 : return false;
891 : :
892 : 12 : QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Read );
893 : 12 : if ( sDisableCache )
894 : 0 : return false;
895 : :
896 : 12 : const QList< QgsCoordinateTransform > values = sTransforms.values( qMakePair( sourceKey, destKey ) );
897 : 12 : for ( auto valIt = values.constBegin(); valIt != values.constEnd(); ++valIt )
898 : : {
899 : 10 : if ( ( *valIt ).coordinateOperation() == coordinateOperationProj && ( *valIt ).allowFallbackTransforms() == allowFallback )
900 : : {
901 : : // need to save, and then restore the context... we don't want this to be cached or to use the values from the cache
902 : 10 : QgsCoordinateTransformContext context = mContext;
903 : : #ifdef QGISDEBUG
904 : : bool hasContext = mHasContext;
905 : : #endif
906 : 10 : *this = *valIt;
907 : 10 : locker.unlock();
908 : :
909 : 10 : mContext = context;
910 : : #ifdef QGISDEBUG
911 : : mHasContext = hasContext;
912 : : #endif
913 : :
914 : 10 : return true;
915 : 10 : }
916 : 0 : }
917 : 2 : return false;
918 : 12 : }
919 : :
920 : 2 : void QgsCoordinateTransform::addToCache()
921 : : {
922 : 2 : if ( !d->mSourceCRS.isValid() || !d->mDestCRS.isValid() )
923 : 0 : return;
924 : :
925 : 4 : const QString sourceKey = d->mSourceCRS.authid().isEmpty() ?
926 : 2 : d->mSourceCRS.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED ) : d->mSourceCRS.authid();
927 : 4 : const QString destKey = d->mDestCRS.authid().isEmpty() ?
928 : 2 : d->mDestCRS.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED ) : d->mDestCRS.authid();
929 : :
930 : 2 : if ( sourceKey.isEmpty() || destKey.isEmpty() )
931 : 0 : return;
932 : :
933 : 2 : QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Write );
934 : 2 : if ( sDisableCache )
935 : 0 : return;
936 : :
937 : 2 : sTransforms.insert( qMakePair( sourceKey, destKey ), *this );
938 : 2 : }
939 : :
940 : 0 : int QgsCoordinateTransform::sourceDatumTransformId() const
941 : : {
942 : : Q_NOWARN_DEPRECATED_PUSH
943 : 0 : return d->mSourceDatumTransform;
944 : : Q_NOWARN_DEPRECATED_POP
945 : : }
946 : :
947 : 0 : void QgsCoordinateTransform::setSourceDatumTransformId( int dt )
948 : : {
949 : 0 : d.detach();
950 : : Q_NOWARN_DEPRECATED_PUSH
951 : 0 : d->mSourceDatumTransform = dt;
952 : : Q_NOWARN_DEPRECATED_POP
953 : 0 : }
954 : :
955 : 0 : int QgsCoordinateTransform::destinationDatumTransformId() const
956 : : {
957 : : Q_NOWARN_DEPRECATED_PUSH
958 : 0 : return d->mDestinationDatumTransform;
959 : : Q_NOWARN_DEPRECATED_POP
960 : : }
961 : :
962 : 0 : void QgsCoordinateTransform::setDestinationDatumTransformId( int dt )
963 : : {
964 : 0 : d.detach();
965 : : Q_NOWARN_DEPRECATED_PUSH
966 : 0 : d->mDestinationDatumTransform = dt;
967 : : Q_NOWARN_DEPRECATED_POP
968 : 0 : }
969 : :
970 : 8 : void QgsCoordinateTransform::invalidateCache( bool disableCache )
971 : : {
972 : 8 : QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Write );
973 : 8 : if ( sDisableCache )
974 : 3 : return;
975 : :
976 : 5 : if ( disableCache )
977 : : {
978 : 5 : sDisableCache = true;
979 : 5 : }
980 : :
981 : 5 : sTransforms.clear();
982 : 8 : }
983 : :
984 : 3 : void QgsCoordinateTransform::removeFromCacheObjectsBelongingToCurrentThread( void *pj_context )
985 : : {
986 : : // Not completely sure about object order destruction after main() has
987 : : // exited. So it is safer to check sDisableCache before using sCacheLock
988 : : // in case sCacheLock would have been destroyed before the current TLS
989 : : // QgsProjContext object that has called us...
990 : 3 : if ( sDisableCache )
991 : 3 : return;
992 : :
993 : 0 : QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Write );
994 : : // cppcheck-suppress identicalConditionAfterEarlyExit
995 : 0 : if ( sDisableCache )
996 : 0 : return;
997 : :
998 : 0 : for ( auto it = sTransforms.begin(); it != sTransforms.end(); )
999 : : {
1000 : 0 : auto &v = it.value();
1001 : 0 : if ( v.d->removeObjectsBelongingToCurrentThread( pj_context ) )
1002 : 0 : it = sTransforms.erase( it );
1003 : : else
1004 : 0 : ++it;
1005 : : }
1006 : 3 : }
1007 : :
1008 : 88 : double QgsCoordinateTransform::scaleFactor( const QgsRectangle &ReferenceExtent ) const
1009 : : {
1010 : 88 : QgsPointXY source1( ReferenceExtent.xMinimum(), ReferenceExtent.yMinimum() );
1011 : 88 : QgsPointXY source2( ReferenceExtent.xMaximum(), ReferenceExtent.yMaximum() );
1012 : 88 : double distSourceUnits = std::sqrt( source1.sqrDist( source2 ) );
1013 : 88 : QgsPointXY dest1 = transform( source1 );
1014 : 88 : QgsPointXY dest2 = transform( source2 );
1015 : 88 : double distDestUnits = std::sqrt( dest1.sqrDist( dest2 ) );
1016 : 88 : return distDestUnits / distSourceUnits;
1017 : : }
1018 : :
1019 : 0 : void QgsCoordinateTransform::setCustomMissingRequiredGridHandler( const std::function<void ( const QgsCoordinateReferenceSystem &, const QgsCoordinateReferenceSystem &, const QgsDatumTransform::GridDetails & )> &handler )
1020 : : {
1021 : 0 : QgsCoordinateTransformPrivate::setCustomMissingRequiredGridHandler( handler );
1022 : 0 : }
1023 : :
1024 : 0 : void QgsCoordinateTransform::setCustomMissingPreferredGridHandler( const std::function<void ( const QgsCoordinateReferenceSystem &, const QgsCoordinateReferenceSystem &, const QgsDatumTransform::TransformDetails &, const QgsDatumTransform::TransformDetails & )> &handler )
1025 : : {
1026 : 0 : QgsCoordinateTransformPrivate::setCustomMissingPreferredGridHandler( handler );
1027 : 0 : }
1028 : :
1029 : 0 : void QgsCoordinateTransform::setCustomCoordinateOperationCreationErrorHandler( const std::function<void ( const QgsCoordinateReferenceSystem &, const QgsCoordinateReferenceSystem &, const QString & )> &handler )
1030 : : {
1031 : 0 : QgsCoordinateTransformPrivate::setCustomCoordinateOperationCreationErrorHandler( handler );
1032 : 0 : }
1033 : :
1034 : 0 : void QgsCoordinateTransform::setCustomMissingGridUsedByContextHandler( const std::function<void ( const QgsCoordinateReferenceSystem &, const QgsCoordinateReferenceSystem &, const QgsDatumTransform::TransformDetails & )> &handler )
1035 : : {
1036 : 0 : QgsCoordinateTransformPrivate::setCustomMissingGridUsedByContextHandler( handler );
1037 : 0 : }
1038 : :
1039 : 0 : void QgsCoordinateTransform::setFallbackOperationOccurredHandler( const std::function<void ( const QgsCoordinateReferenceSystem &, const QgsCoordinateReferenceSystem &, const QString & )> &handler )
1040 : : {
1041 : 0 : sFallbackOperationOccurredHandler = handler;
1042 : 0 : }
|
