QGIS API Documentation 3.41.0-Master (fda2aa46e9a)
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 std::vector<double> x( nXPoints * static_cast< std::size_t >( nYPoints ) );
636 std::vector<double> y( nXPoints * static_cast< std::size_t >( nYPoints ) );
637 std::vector<double> z( nXPoints * static_cast< std::size_t >( nYPoints ) );
638
639 QgsDebugMsgLevel( QStringLiteral( "Entering transformBoundingBox..." ), 4 );
640
641 // Populate the vectors
642
643 const double dx = rect.width() / static_cast< double >( nXPoints - 1 );
644 const double dy = ( yMax - yMin ) / static_cast< double >( nYPoints - 1 );
645
646 double pointY = yMin;
647
648 for ( int i = 0; i < nYPoints ; i++ )
649 {
650
651 // Start at right edge
652 double pointX = rect.xMinimum();
653
654 for ( int j = 0; j < nXPoints; j++ )
655 {
656 x[( i * nXPoints ) + j] = pointX;
657 y[( i * nXPoints ) + j] = pointY;
658 // and the height...
659 z[( i * nXPoints ) + j] = 0.0;
660 // QgsDebugMsgLevel(QString("BBox coord: (%1, %2)").arg(x[(i*numP) + j]).arg(y[(i*numP) + j]), 2);
661 pointX += dx;
662 }
663 pointY += dy;
664 }
665
666 // Do transformation. Any exception generated must
667 // be handled in above layers.
668 try
669 {
670 transformCoords( nXPoints * nYPoints, x.data(), y.data(), z.data(), direction );
671 }
672 catch ( const QgsCsException & )
673 {
674 // rethrow the exception
675 QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
676 throw;
677 }
678
679 // 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°
680 bool doHandle180Crossover = false;
681 if ( nXPoints > 0 )
682 {
683 const double xMin = std::fmod( x[0], 180.0 );
684 const double xMax = std::fmod( x[nXPoints - 1], 180.0 );
685 if ( handle180Crossover
686 && ( ( direction == Qgis::TransformDirection::Forward && d->mDestCRS.isGeographic() ) ||
687 ( direction == Qgis::TransformDirection::Reverse && d->mSourceCRS.isGeographic() ) )
688 && xMin > 0.0 && xMin <= 180.0 && xMax < 0.0 && xMax >= -180.0 )
689 {
690 doHandle180Crossover = true;
691 }
692 }
693
694 // Calculate the bounding box and use that for the extent
695 for ( int i = 0; i < nXPoints * nYPoints; i++ )
696 {
697 if ( !std::isfinite( x[i] ) || !std::isfinite( y[i] ) )
698 {
699 continue;
700 }
701
702 if ( doHandle180Crossover )
703 {
704 //if crossing the date line, temporarily add 360 degrees to -ve longitudes
705 bb_rect.combineExtentWith( x[i] >= 0.0 ? x[i] : x[i] + 360.0, y[i] );
706 }
707 else
708 {
709 bb_rect.combineExtentWith( x[i], y[i] );
710 }
711 }
712
713 if ( bb_rect.isNull() )
714 {
715 // something bad happened when reprojecting the filter rect... no finite points were left!
716 throw QgsCsException( QObject::tr( "Could not transform bounding box to target CRS" ) );
717 }
718
719 if ( doHandle180Crossover )
720 {
721 //subtract temporary addition of 360 degrees from longitudes
722 if ( bb_rect.xMinimum() > 180.0 )
723 bb_rect.setXMinimum( bb_rect.xMinimum() - 360.0 );
724 if ( bb_rect.xMaximum() > 180.0 )
725 bb_rect.setXMaximum( bb_rect.xMaximum() - 360.0 );
726 }
727
728 QgsDebugMsgLevel( "Projected extent: " + bb_rect.toString(), 4 );
729
730 if ( bb_rect.isEmpty() )
731 {
732 QgsDebugMsgLevel( "Original extent: " + rect.toString(), 4 );
733 }
734
735 return bb_rect;
736}
737
738void QgsCoordinateTransform::transformCoords( int numPoints, double *x, double *y, double *z, Qgis::TransformDirection direction ) const
739{
740 if ( !d->mIsValid || d->mShortCircuit )
741 return;
742 // Refuse to transform the points if the srs's are invalid
743 if ( !d->mSourceCRS.isValid() )
744 {
745 QgsMessageLog::logMessage( QObject::tr( "The source spatial reference system (CRS) is not valid. "
746 "The coordinates can not be reprojected. The CRS is: %1" )
747 .arg( d->mSourceCRS.toProj() ), QObject::tr( "CRS" ) );
748 return;
749 }
750 if ( !d->mDestCRS.isValid() )
751 {
752 QgsMessageLog::logMessage( QObject::tr( "The destination spatial reference system (CRS) is not valid. "
753 "The coordinates can not be reprojected. The CRS is: %1" ).arg( d->mDestCRS.toProj() ), QObject::tr( "CRS" ) );
754 return;
755 }
756
757 std::vector< int > zNanPositions;
758 for ( int i = 0; i < numPoints; i++ )
759 {
760 if ( std::isnan( z[i] ) )
761 {
762 zNanPositions.push_back( i );
763 z[i] = 0.0;
764 }
765 }
766
767 std::vector< double > xprev( numPoints );
768 memcpy( xprev.data(), x, sizeof( double ) * numPoints );
769 std::vector< double > yprev( numPoints );
770 memcpy( yprev.data(), y, sizeof( double ) * numPoints );
771 std::vector< double > zprev( numPoints );
772 memcpy( zprev.data(), z, sizeof( double ) * numPoints );
773
774 const bool useTime = !std::isnan( d->mDefaultTime );
775 std::vector< double > t( useTime ? numPoints : 0, d->mDefaultTime );
776
777#ifdef COORDINATE_TRANSFORM_VERBOSE
778 double xorg = *x;
779 double yorg = *y;
780 QgsDebugMsgLevel( QStringLiteral( "[[[[[[ Number of points to transform: %1 ]]]]]]" ).arg( numPoints ), 2 );
781#endif
782
783#ifdef QGISDEBUG
784 if ( !mHasContext )
785 QgsDebugMsgLevel( QStringLiteral( "No QgsCoordinateTransformContext context set for transform" ), 4 );
786#endif
787
788 // use proj4 to do the transform
789
790 // if the source/destination projection is lat/long, convert the points to radians
791 // prior to transforming
792 ProjData projData = d->threadLocalProjData();
793
794 int projResult = 0;
795
796 proj_errno_reset( projData );
797 proj_trans_generic( projData, ( direction == Qgis::TransformDirection::Forward && !d->mIsReversed ) || ( direction == Qgis::TransformDirection::Reverse && d->mIsReversed ) ? PJ_FWD : PJ_INV,
798 x, sizeof( double ), numPoints,
799 y, sizeof( double ), numPoints,
800 z, sizeof( double ), numPoints,
801 useTime ? t.data() : nullptr, sizeof( double ), useTime ? numPoints : 0 );
802 // Try to - approximately - emulate the behavior of pj_transform()...
803 // In the case of a single point transform, and a transformation error occurs,
804 // pj_transform() would return the errno. In cases of multiple point transform,
805 // it would continue (for non-transient errors, that is pipeline definition
806 // errors) and just set the resulting x,y to infinity. This is in fact a
807 // bit more subtle than that, and I'm not completely sure the logic in
808 // pj_transform() was really sane & fully bullet proof
809 // So here just check proj_errno() for single point transform
810 int actualRes = 0;
811 if ( numPoints == 1 )
812 {
813 projResult = proj_errno( projData );
814 actualRes = projResult;
815 }
816 else
817 {
818 actualRes = proj_errno( projData );
819 }
820 if ( actualRes == 0 )
821 {
822 // proj_errno is sometimes not an accurate method to test for transform failures - so we need to
823 // manually scan for nan values
824 if ( std::any_of( x, x + numPoints, []( double v ) { return std::isinf( v ); } )
825 || std::any_of( y, y + numPoints, []( double v ) { return std::isinf( v ); } )
826 || std::any_of( z, z + numPoints, []( double v ) { return std::isinf( v ); } ) )
827 {
828 actualRes = 1;
829 }
830 }
831
832 mFallbackOperationOccurred = false;
833 if ( actualRes != 0
834 && ( 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
835 && ( d->mAllowFallbackTransforms || mBallparkTransformsAreAppropriate ) )
836 {
837 // fail #1 -- try with getting proj to auto-pick an appropriate coordinate operation for the points
838 if ( PJ *transform = d->threadLocalFallbackProjData() )
839 {
840 projResult = 0;
841 proj_errno_reset( transform );
842 proj_trans_generic( transform, direction == Qgis::TransformDirection::Forward ? PJ_FWD : PJ_INV,
843 xprev.data(), sizeof( double ), numPoints,
844 yprev.data(), sizeof( double ), numPoints,
845 zprev.data(), sizeof( double ), numPoints,
846 useTime ? t.data() : nullptr, sizeof( double ), useTime ? numPoints : 0 );
847 // Try to - approximately - emulate the behavior of pj_transform()...
848 // In the case of a single point transform, and a transformation error occurs,
849 // pj_transform() would return the errno. In cases of multiple point transform,
850 // it would continue (for non-transient errors, that is pipeline definition
851 // errors) and just set the resulting x,y to infinity. This is in fact a
852 // bit more subtle than that, and I'm not completely sure the logic in
853 // pj_transform() was really sane & fully bullet proof
854 // So here just check proj_errno() for single point transform
855 if ( numPoints == 1 )
856 {
857 // hmm - something very odd here. We can't trust proj_errno( transform ), as that's giving us incorrect error numbers
858 // (such as "failed to load datum shift file", which is definitely incorrect for a default proj created operation!)
859 // so we resort to testing values ourselves...
860 projResult = std::isinf( xprev[0] ) || std::isinf( yprev[0] ) || std::isinf( zprev[0] ) ? 1 : 0;
861 }
862
863 if ( projResult == 0 )
864 {
865 memcpy( x, xprev.data(), sizeof( double ) * numPoints );
866 memcpy( y, yprev.data(), sizeof( double ) * numPoints );
867 memcpy( z, zprev.data(), sizeof( double ) * numPoints );
868 mFallbackOperationOccurred = true;
869 }
870
871 if ( !mBallparkTransformsAreAppropriate && !mDisableFallbackHandler && sFallbackOperationOccurredHandler )
872 {
873 sFallbackOperationOccurredHandler( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation );
874#if 0
875 const QString warning = QStringLiteral( "A fallback coordinate operation was used between %1 and %2" ).arg( d->mSourceCRS.authid(),
876 d->mDestCRS.authid() );
877 qWarning( "%s", warning.toLatin1().constData() );
878#endif
879 }
880 }
881 }
882
883 for ( const int &pos : zNanPositions )
884 {
885 z[pos] = std::numeric_limits<double>::quiet_NaN();
886 }
887
888 if ( projResult != 0 )
889 {
890 //something bad happened....
891 QString points;
892
893 for ( int i = 0; i < numPoints; ++i )
894 {
895 points += QStringLiteral( "(%1, %2)\n" ).arg( x[i], 0, 'f' ).arg( y[i], 0, 'f' );
896 }
897
898 const QString dir = ( direction == Qgis::TransformDirection::Forward ) ? QObject::tr( "forward transform" ) : QObject::tr( "inverse transform" );
899
900 const QString msg = QObject::tr( "%1 of\n"
901 "%2"
902 "Error: %3" )
903 .arg( dir,
904 points,
905 projResult < 0 ? QString::fromUtf8( proj_errno_string( projResult ) ) : QObject::tr( "Fallback transform failed" ) );
906
907
908 // don't flood console with thousands of duplicate transform error messages
909 if ( msg != mLastError )
910 {
911 QgsDebugError( "Projection failed emitting invalid transform signal: " + msg );
912 mLastError = msg;
913 }
914 QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
915
916 throw QgsCsException( msg );
917 }
918
919#ifdef COORDINATE_TRANSFORM_VERBOSE
920 QgsDebugMsgLevel( QStringLiteral( "[[[[[[ Projected %1, %2 to %3, %4 ]]]]]]" )
921 .arg( xorg, 0, 'g', 15 ).arg( yorg, 0, 'g', 15 )
922 .arg( *x, 0, 'g', 15 ).arg( *y, 0, 'g', 15 ), 2 );
923#endif
924}
925
927{
928 return d->mIsValid;
929}
930
932{
933 return !d->mIsValid || d->mShortCircuit;
934}
935
937{
938 return d->mIsValid && d->mHasVerticalComponent;
939}
940
942{
943 return d->mProjCoordinateOperation;
944}
945
947{
948 ProjData projData = d->threadLocalProjData();
950}
951
952void QgsCoordinateTransform::setCoordinateOperation( const QString &operation ) const
953{
954 d.detach();
955 d->mProjCoordinateOperation = operation;
956 d->mShouldReverseCoordinateOperation = false;
957}
958
960{
961 d.detach();
962 d->mAllowFallbackTransforms = allowed;
963}
964
966{
967 return d->mAllowFallbackTransforms;
968}
969
971{
972 mBallparkTransformsAreAppropriate = appropriate;
973}
974
976{
977 mDisableFallbackHandler = disabled;
978}
979
981{
982 return mFallbackOperationOccurred;
983}
984
985const char *finder( const char *name )
986{
987 QString proj;
988#ifdef Q_OS_WIN
989 proj = QApplication::applicationDirPath()
990 + "/share/proj/" + QString( name );
991#else
992 Q_UNUSED( name )
993#endif
994 return proj.toUtf8();
995}
996
997bool QgsCoordinateTransform::setFromCache( const QgsCoordinateReferenceSystem &src, const QgsCoordinateReferenceSystem &dest, const QString &coordinateOperationProj, bool allowFallback )
998{
999 if ( !src.isValid() || !dest.isValid() )
1000 return false;
1001
1002 const QString sourceKey = src.authid().isEmpty() ?
1004 const QString destKey = dest.authid().isEmpty() ?
1006
1007 if ( sourceKey.isEmpty() || destKey.isEmpty() )
1008 return false;
1009
1010 QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Read );
1011 if ( sDisableCache )
1012 return false;
1013
1014 const QList< QgsCoordinateTransform > values = sTransforms.values( qMakePair( sourceKey, destKey ) );
1015 for ( auto valIt = values.constBegin(); valIt != values.constEnd(); ++valIt )
1016 {
1017 if ( ( *valIt ).coordinateOperation() == coordinateOperationProj
1018 && ( *valIt ).allowFallbackTransforms() == allowFallback
1019 && qgsNanCompatibleEquals( src.coordinateEpoch(), ( *valIt ).sourceCrs().coordinateEpoch() )
1020 && qgsNanCompatibleEquals( dest.coordinateEpoch(), ( *valIt ).destinationCrs().coordinateEpoch() )
1021 )
1022 {
1023 // need to save, and then restore the context... we don't want this to be cached or to use the values from the cache
1024 const QgsCoordinateTransformContext context = mContext;
1025#ifdef QGISDEBUG
1026 const bool hasContext = mHasContext;
1027#endif
1028 *this = *valIt;
1029 locker.unlock();
1030
1031 mContext = context;
1032#ifdef QGISDEBUG
1033 mHasContext = hasContext;
1034#endif
1035
1036 return true;
1037 }
1038 }
1039 return false;
1040}
1041
1042void QgsCoordinateTransform::addToCache()
1043{
1044 if ( !d->mSourceCRS.isValid() || !d->mDestCRS.isValid() )
1045 return;
1046
1047 const QString sourceKey = d->mSourceCRS.authid().isEmpty() ?
1048 d->mSourceCRS.toWkt( Qgis::CrsWktVariant::Preferred ) : d->mSourceCRS.authid();
1049 const QString destKey = d->mDestCRS.authid().isEmpty() ?
1050 d->mDestCRS.toWkt( Qgis::CrsWktVariant::Preferred ) : d->mDestCRS.authid();
1051
1052 if ( sourceKey.isEmpty() || destKey.isEmpty() )
1053 return;
1054
1055 const QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Write );
1056 if ( sDisableCache )
1057 return;
1058
1059 sTransforms.insert( qMakePair( sourceKey, destKey ), *this );
1060}
1061
1063{
1065 return d->mSourceDatumTransform;
1067}
1068
1070{
1071 d.detach();
1073 d->mSourceDatumTransform = dt;
1075}
1076
1078{
1080 return d->mDestinationDatumTransform;
1082}
1083
1085{
1086 d.detach();
1088 d->mDestinationDatumTransform = dt;
1090}
1091
1093{
1094 const QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Write );
1095 if ( sDisableCache )
1096 return;
1097
1098 if ( disableCache )
1099 {
1100 sDisableCache = true;
1101 }
1102
1103 sTransforms.clear();
1104}
1105
1106void QgsCoordinateTransform::removeFromCacheObjectsBelongingToCurrentThread( void *pj_context )
1107{
1108 // Not completely sure about object order destruction after main() has
1109 // exited. So it is safer to check sDisableCache before using sCacheLock
1110 // in case sCacheLock would have been destroyed before the current TLS
1111 // QgsProjContext object that has called us...
1112 if ( sDisableCache )
1113 return;
1114
1115 const QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Write );
1116 // cppcheck-suppress identicalConditionAfterEarlyExit
1117 if ( sDisableCache )
1118 return;
1119
1120 for ( auto it = sTransforms.begin(); it != sTransforms.end(); )
1121 {
1122 auto &v = it.value();
1123 if ( v.d->removeObjectsBelongingToCurrentThread( pj_context ) )
1124 it = sTransforms.erase( it );
1125 else
1126 ++it;
1127 }
1128}
1129
1130double QgsCoordinateTransform::scaleFactor( const QgsRectangle &ReferenceExtent ) const
1131{
1132 const QgsPointXY source1( ReferenceExtent.xMinimum(), ReferenceExtent.yMinimum() );
1133 const QgsPointXY source2( ReferenceExtent.xMaximum(), ReferenceExtent.yMaximum() );
1134 const double distSourceUnits = std::sqrt( source1.sqrDist( source2 ) );
1135 const QgsPointXY dest1 = transform( source1 );
1136 const QgsPointXY dest2 = transform( source2 );
1137 const double distDestUnits = std::sqrt( dest1.sqrDist( dest2 ) );
1138 return distDestUnits / distSourceUnits;
1139}
1140
1142{
1143 QgsCoordinateTransformPrivate::setCustomMissingRequiredGridHandler( handler );
1144}
1145
1147{
1148 QgsCoordinateTransformPrivate::setCustomMissingPreferredGridHandler( handler );
1149}
1150
1152{
1153 QgsCoordinateTransformPrivate::setCustomCoordinateOperationCreationErrorHandler( handler );
1154}
1155
1157{
1158 QgsCoordinateTransformPrivate::setCustomMissingGridUsedByContextHandler( handler );
1159}
1160
1161void QgsCoordinateTransform::setFallbackOperationOccurredHandler( const std::function<void ( const QgsCoordinateReferenceSystem &, const QgsCoordinateReferenceSystem &, const QString & )> &handler )
1162{
1163 sFallbackOperationOccurredHandler = handler;
1164}
1165
1167{
1168 QgsCoordinateTransformPrivate::setDynamicCrsToDynamicCrsWarningHandler( handler );
1169}
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:6494
#define Q_NOWARN_DEPRECATED_PUSH
Definition qgis.h:6493
bool qgsNanCompatibleEquals(double a, double b)
Compare two doubles, treating nan values as equal.
Definition qgis.h:5879
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.