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