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