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