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