QGIS API Documentation  3.21.0-Master (909859188c)
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 
29 //qt includes
30 #include <QDomNode>
31 #include <QDomElement>
32 #include <QApplication>
33 #include <QPolygonF>
34 #include <QStringList>
35 #include <QVector>
36 
37 #include <proj.h>
38 #include "qgsprojutils.h"
39 
40 #include <sqlite3.h>
41 #include <qlogging.h>
42 #include <vector>
43 #include <algorithm>
44 
45 // if defined shows all information about transform to stdout
46 // #define COORDINATE_TRANSFORM_VERBOSE
47 
48 QReadWriteLock QgsCoordinateTransform::sCacheLock;
49 QMultiHash< QPair< QString, QString >, QgsCoordinateTransform > QgsCoordinateTransform::sTransforms; //same auth_id pairs might have different datum transformations
50 bool QgsCoordinateTransform::sDisableCache = false;
51 
52 std::function< void( const QgsCoordinateReferenceSystem &sourceCrs,
53  const QgsCoordinateReferenceSystem &destinationCrs,
54  const QString &desiredOperation )> QgsCoordinateTransform::sFallbackOperationOccurredHandler = nullptr;
55 
57 {
58  d = new QgsCoordinateTransformPrivate();
59 }
60 
62 {
63  mContext = context;
64  d = new QgsCoordinateTransformPrivate( source, destination, mContext );
65 #ifdef QGISDEBUG
66  mHasContext = true;
67 #endif
68 
69  if ( !d->checkValidity() )
70  return;
71 
73  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation, d->mAllowFallbackTransforms ) )
74  {
75  d->initialize();
76  addToCache();
77  }
79 }
80 
82 {
83  mContext = project ? project->transformContext() : QgsCoordinateTransformContext();
84  d = new QgsCoordinateTransformPrivate( source, destination, mContext );
85 #ifdef QGISDEBUG
86  if ( project )
87  mHasContext = true;
88 #endif
89 
90  if ( !d->checkValidity() )
91  return;
92 
94  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation, d->mAllowFallbackTransforms ) )
95  {
96  d->initialize();
97  addToCache();
98  }
100 }
101 
102 QgsCoordinateTransform::QgsCoordinateTransform( const QgsCoordinateReferenceSystem &source, const QgsCoordinateReferenceSystem &destination, int sourceDatumTransform, int destinationDatumTransform )
103 {
104  d = new QgsCoordinateTransformPrivate( source, destination, sourceDatumTransform, destinationDatumTransform );
105 #ifdef QGISDEBUG
106  mHasContext = true; // not strictly true, but we don't need to worry if datums have been explicitly set
107 #endif
108 
109  if ( !d->checkValidity() )
110  return;
111 
113  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation, d->mAllowFallbackTransforms ) )
114  {
115  d->initialize();
116  addToCache();
117  }
119 }
120 
122  : mContext( o.mContext )
123 #ifdef QGISDEBUG
124  , mHasContext( o.mHasContext )
125 #endif
126  , mLastError()
127 {
128  d = o.d;
129 }
130 
132 {
133  d = o.d;
134 #ifdef QGISDEBUG
135  mHasContext = o.mHasContext;
136 #endif
137  mContext = o.mContext;
138  mLastError = QString();
139  return *this;
140 }
141 
143 
145 {
146  d.detach();
147  d->mSourceCRS = crs;
148  if ( !d->checkValidity() )
149  return;
150 
151  d->calculateTransforms( mContext );
153  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation, d->mAllowFallbackTransforms ) )
154  {
155  d->initialize();
156  addToCache();
157  }
159 }
161 {
162  d.detach();
163  d->mDestCRS = crs;
164  if ( !d->checkValidity() )
165  return;
166 
167  d->calculateTransforms( mContext );
169  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation, d->mAllowFallbackTransforms ) )
170  {
171  d->initialize();
172  addToCache();
173  }
175 }
176 
178 {
179  d.detach();
180  mContext = context;
181 #ifdef QGISDEBUG
182  mHasContext = true;
183 #endif
184  if ( !d->checkValidity() )
185  return;
186 
187  d->calculateTransforms( mContext );
189  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation, d->mAllowFallbackTransforms ) )
190  {
191  d->initialize();
192  addToCache();
193  }
195 }
196 
198 {
199  return mContext;
200 }
201 
203 {
204  return d->mSourceCRS;
205 }
206 
208 {
209  return d->mDestCRS;
210 }
211 
213 {
214  if ( !d->mIsValid || d->mShortCircuit )
215  return point;
216 
217  // transform x
218  double x = point.x();
219  double y = point.y();
220  double z = 0.0;
221  try
222  {
223  transformCoords( 1, &x, &y, &z, direction );
224  }
225  catch ( const QgsCsException & )
226  {
227  // rethrow the exception
228  QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
229  throw;
230  }
231 
232  return QgsPointXY( x, y );
233 }
234 
235 
236 QgsPointXY QgsCoordinateTransform::transform( const double theX, const double theY = 0.0, TransformDirection direction ) const
237 {
238  try
239  {
240  return transform( QgsPointXY( theX, theY ), direction );
241  }
242  catch ( const QgsCsException & )
243  {
244  // rethrow the exception
245  QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
246  throw;
247  }
248 }
249 
251 {
252  if ( !d->mIsValid || d->mShortCircuit )
253  return rect;
254  // transform x
255  double x1 = rect.xMinimum();
256  double y1 = rect.yMinimum();
257  double x2 = rect.xMaximum();
258  double y2 = rect.yMaximum();
259 
260  // Number of points to reproject------+
261  // |
262  // V
263  try
264  {
265  double z = 0.0;
266  transformCoords( 1, &x1, &y1, &z, direction );
267  transformCoords( 1, &x2, &y2, &z, direction );
268  }
269  catch ( const QgsCsException & )
270  {
271  // rethrow the exception
272  QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
273  throw;
274  }
275 
276 #ifdef COORDINATE_TRANSFORM_VERBOSE
277  QgsDebugMsg( QStringLiteral( "Rect projection..." ) );
278  QgsDebugMsg( QStringLiteral( "Xmin : %1 --> %2" ).arg( rect.xMinimum() ).arg( x1 ) );
279  QgsDebugMsg( QStringLiteral( "Ymin : %1 --> %2" ).arg( rect.yMinimum() ).arg( y1 ) );
280  QgsDebugMsg( QStringLiteral( "Xmax : %1 --> %2" ).arg( rect.xMaximum() ).arg( x2 ) );
281  QgsDebugMsg( QStringLiteral( "Ymax : %1 --> %2" ).arg( rect.yMaximum() ).arg( y2 ) );
282 #endif
283  return QgsRectangle( x1, y1, x2, y2 );
284 }
285 
287 {
288  double x = point.x();
289  double y = point.y();
290  double z = point.z();
291  try
292  {
293  transformCoords( 1, &x, &y, &z, direction );
294  }
295  catch ( const QgsCsException & )
296  {
297  // rethrow the exception
298  QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
299  throw;
300  }
301  return QgsVector3D( x, y, z );
302 }
303 
304 void QgsCoordinateTransform::transformInPlace( double &x, double &y, double &z,
305  TransformDirection direction ) const
306 {
307  if ( !d->mIsValid || d->mShortCircuit )
308  return;
309 #ifdef QGISDEBUG
310 // QgsDebugMsg(QString("Using transform in place %1 %2").arg(__FILE__).arg(__LINE__));
311 #endif
312  // transform x
313  try
314  {
315  transformCoords( 1, &x, &y, &z, direction );
316  }
317  catch ( const QgsCsException & )
318  {
319  // rethrow the exception
320  QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
321  throw;
322  }
323 }
324 
325 void QgsCoordinateTransform::transformInPlace( float &x, float &y, double &z,
326  TransformDirection direction ) const
327 {
328  double xd = static_cast< double >( x ), yd = static_cast< double >( y );
329  transformInPlace( xd, yd, z, direction );
330  x = xd;
331  y = yd;
332 }
333 
334 void QgsCoordinateTransform::transformInPlace( float &x, float &y, float &z,
335  TransformDirection direction ) const
336 {
337  if ( !d->mIsValid || d->mShortCircuit )
338  return;
339 #ifdef QGISDEBUG
340  // QgsDebugMsg(QString("Using transform in place %1 %2").arg(__FILE__).arg(__LINE__));
341 #endif
342  // transform x
343  try
344  {
345  double xd = x;
346  double yd = y;
347  double zd = z;
348  transformCoords( 1, &xd, &yd, &zd, direction );
349  x = xd;
350  y = yd;
351  z = zd;
352  }
353  catch ( QgsCsException & )
354  {
355  // rethrow the exception
356  QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
357  throw;
358  }
359 }
360 
361 void QgsCoordinateTransform::transformPolygon( QPolygonF &poly, TransformDirection direction ) const
362 {
363  if ( !d->mIsValid || d->mShortCircuit )
364  {
365  return;
366  }
367 
368  //create x, y arrays
369  const int nVertices = poly.size();
370 
371  QVector<double> x( nVertices );
372  QVector<double> y( nVertices );
373  QVector<double> z( nVertices );
374  double *destX = x.data();
375  double *destY = y.data();
376  double *destZ = z.data();
377 
378  const QPointF *polyData = poly.constData();
379  for ( int i = 0; i < nVertices; ++i )
380  {
381  *destX++ = polyData->x();
382  *destY++ = polyData->y();
383  *destZ++ = 0;
384  polyData++;
385  }
386 
387  QString err;
388  try
389  {
390  transformCoords( nVertices, x.data(), y.data(), z.data(), direction );
391  }
392  catch ( const QgsCsException &e )
393  {
394  // record the exception, but don't rethrow it until we've recorded the coordinates we *could* transform
395  err = e.what();
396  }
397 
398  QPointF *destPoint = poly.data();
399  const double *srcX = x.constData();
400  const double *srcY = y.constData();
401  for ( int i = 0; i < nVertices; ++i )
402  {
403  destPoint->rx() = *srcX++;
404  destPoint->ry() = *srcY++;
405  destPoint++;
406  }
407 
408  // rethrow the exception
409  if ( !err.isEmpty() )
410  throw QgsCsException( err );
411 }
412 
414  QVector<double> &x, QVector<double> &y, QVector<double> &z,
415  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 
442  QVector<float> &x, QVector<float> &y, QVector<float> &z,
443  TransformDirection direction ) const
444 {
445  if ( !d->mIsValid || d->mShortCircuit )
446  return;
447 
448  Q_ASSERT( x.size() == y.size() );
449 
450  // Apparently, if one has a std::vector, it is valid to use the
451  // address of the first element in the vector as a pointer to an
452  // array of the vectors data, and hence easily interface with code
453  // that wants C-style arrays.
454 
455  try
456  {
457  //copy everything to double vectors since proj needs double
458  const int vectorSize = x.size();
459  QVector<double> xd( x.size() );
460  QVector<double> yd( y.size() );
461  QVector<double> zd( z.size() );
462 
463  double *destX = xd.data();
464  double *destY = yd.data();
465  double *destZ = zd.data();
466 
467  const float *srcX = x.constData();
468  const float *srcY = y.constData();
469  const float *srcZ = z.constData();
470 
471  for ( int i = 0; i < vectorSize; ++i )
472  {
473  *destX++ = static_cast< double >( *srcX++ );
474  *destY++ = static_cast< double >( *srcY++ );
475  *destZ++ = static_cast< double >( *srcZ++ );
476  }
477 
478  transformCoords( x.size(), &xd[0], &yd[0], &zd[0], direction );
479 
480  //copy back
481  float *destFX = x.data();
482  float *destFY = y.data();
483  float *destFZ = z.data();
484  const double *srcXD = xd.constData();
485  const double *srcYD = yd.constData();
486  const double *srcZD = zd.constData();
487  for ( int i = 0; i < vectorSize; ++i )
488  {
489  *destFX++ = static_cast< float >( *srcXD++ );
490  *destFY++ = static_cast< float >( *srcYD++ );
491  *destFZ++ = static_cast< float >( *srcZD++ );
492  }
493  }
494  catch ( QgsCsException & )
495  {
496  // rethrow the exception
497  QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
498  throw;
499  }
500 }
501 
502 QgsRectangle QgsCoordinateTransform::transformBoundingBox( const QgsRectangle &rect, TransformDirection direction, const bool handle180Crossover ) const
503 {
504  // Calculate the bounding box of a QgsRectangle in the source CRS
505  // when projected to the destination CRS (or the inverse).
506  // This is done by looking at a number of points spread evenly
507  // across the rectangle
508 
509  if ( !d->mIsValid || d->mShortCircuit )
510  return rect;
511 
512  if ( rect.isEmpty() )
513  {
514  const QgsPointXY p = transform( rect.xMinimum(), rect.yMinimum(), direction );
515  return QgsRectangle( p, p );
516  }
517 
518  // 64 points (<=2.12) is not enough, see #13665, for EPSG:4326 -> EPSG:3574 (say that it is a hard one),
519  // are decent result from about 500 points and more. This method is called quite often, but
520  // even with 1000 points it takes < 1ms.
521  // TODO: how to effectively and precisely reproject bounding box?
522  const int nPoints = 1000;
523  const double d = std::sqrt( ( rect.width() * rect.height() ) / std::pow( std::sqrt( static_cast< double >( nPoints ) ) - 1, 2.0 ) );
524  const int nXPoints = std::min( static_cast< int >( std::ceil( rect.width() / d ) ) + 1, 1000 );
525  const int nYPoints = std::min( static_cast< int >( std::ceil( rect.height() / d ) ) + 1, 1000 );
526 
527  QgsRectangle bb_rect;
528  bb_rect.setMinimal();
529 
530  // We're interfacing with C-style vectors in the
531  // end, so let's do C-style vectors here too.
532  QVector<double> x( nXPoints * nYPoints );
533  QVector<double> y( nXPoints * nYPoints );
534  QVector<double> z( nXPoints * nYPoints );
535 
536  QgsDebugMsgLevel( QStringLiteral( "Entering transformBoundingBox..." ), 4 );
537 
538  // Populate the vectors
539 
540  const double dx = rect.width() / static_cast< double >( nXPoints - 1 );
541  const double dy = rect.height() / static_cast< double >( nYPoints - 1 );
542 
543  double pointY = rect.yMinimum();
544 
545  for ( int i = 0; i < nYPoints ; i++ )
546  {
547 
548  // Start at right edge
549  double pointX = rect.xMinimum();
550 
551  for ( int j = 0; j < nXPoints; j++ )
552  {
553  x[( i * nXPoints ) + j] = pointX;
554  y[( i * nXPoints ) + j] = pointY;
555  // and the height...
556  z[( i * nXPoints ) + j] = 0.0;
557  // QgsDebugMsg(QString("BBox coord: (%1, %2)").arg(x[(i*numP) + j]).arg(y[(i*numP) + j]));
558  pointX += dx;
559  }
560  pointY += dy;
561  }
562 
563  // Do transformation. Any exception generated must
564  // be handled in above layers.
565  try
566  {
567  transformCoords( nXPoints * nYPoints, x.data(), y.data(), z.data(), direction );
568  }
569  catch ( const QgsCsException & )
570  {
571  // rethrow the exception
572  QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
573  throw;
574  }
575 
576  // Calculate the bounding box and use that for the extent
577 
578  for ( int i = 0; i < nXPoints * nYPoints; i++ )
579  {
580  if ( !std::isfinite( x[i] ) || !std::isfinite( y[i] ) )
581  {
582  continue;
583  }
584 
585  if ( handle180Crossover )
586  {
587  //if crossing the date line, temporarily add 360 degrees to -ve longitudes
588  bb_rect.combineExtentWith( x[i] >= 0.0 ? x[i] : x[i] + 360.0, y[i] );
589  }
590  else
591  {
592  bb_rect.combineExtentWith( x[i], y[i] );
593  }
594  }
595 
596  if ( bb_rect.isNull() )
597  {
598  // something bad happened when reprojecting the filter rect... no finite points were left!
599  throw QgsCsException( QObject::tr( "Could not transform bounding box to target CRS" ) );
600  }
601 
602  if ( handle180Crossover )
603  {
604  //subtract temporary addition of 360 degrees from longitudes
605  if ( bb_rect.xMinimum() > 180.0 )
606  bb_rect.setXMinimum( bb_rect.xMinimum() - 360.0 );
607  if ( bb_rect.xMaximum() > 180.0 )
608  bb_rect.setXMaximum( bb_rect.xMaximum() - 360.0 );
609  }
610 
611  QgsDebugMsgLevel( "Projected extent: " + bb_rect.toString(), 4 );
612 
613  if ( bb_rect.isEmpty() )
614  {
615  QgsDebugMsgLevel( "Original extent: " + rect.toString(), 4 );
616  }
617 
618  return bb_rect;
619 }
620 
621 void QgsCoordinateTransform::transformCoords( int numPoints, double *x, double *y, double *z, TransformDirection direction ) const
622 {
623  if ( !d->mIsValid || d->mShortCircuit )
624  return;
625  // Refuse to transform the points if the srs's are invalid
626  if ( !d->mSourceCRS.isValid() )
627  {
628  QgsMessageLog::logMessage( QObject::tr( "The source spatial reference system (CRS) is not valid. "
629  "The coordinates can not be reprojected. The CRS is: %1" )
630  .arg( d->mSourceCRS.toProj() ), QObject::tr( "CRS" ) );
631  return;
632  }
633  if ( !d->mDestCRS.isValid() )
634  {
635  QgsMessageLog::logMessage( QObject::tr( "The destination spatial reference system (CRS) is not valid. "
636  "The coordinates can not be reprojected. The CRS is: %1" ).arg( d->mDestCRS.toProj() ), QObject::tr( "CRS" ) );
637  return;
638  }
639 
640  std::vector< int > zNanPositions;
641  for ( int i = 0; i < numPoints; i++ )
642  {
643  if ( std::isnan( z[i] ) )
644  {
645  zNanPositions.push_back( i );
646  z[i] = 0.0;
647  }
648  }
649 
650  std::vector< double > xprev( numPoints );
651  memcpy( xprev.data(), x, sizeof( double ) * numPoints );
652  std::vector< double > yprev( numPoints );
653  memcpy( yprev.data(), y, sizeof( double ) * numPoints );
654  std::vector< double > zprev( numPoints );
655  memcpy( zprev.data(), z, sizeof( double ) * numPoints );
656 
657  const bool useTime = !std::isnan( d->mDefaultTime );
658  std::vector< double > t( useTime ? numPoints : 0, d->mDefaultTime );
659 
660 #ifdef COORDINATE_TRANSFORM_VERBOSE
661  double xorg = *x;
662  double yorg = *y;
663  QgsDebugMsg( QStringLiteral( "[[[[[[ Number of points to transform: %1 ]]]]]]" ).arg( numPoints ) );
664 #endif
665 
666 #ifdef QGISDEBUG
667  if ( !mHasContext )
668  QgsDebugMsgLevel( QStringLiteral( "No QgsCoordinateTransformContext context set for transform" ), 4 );
669 #endif
670 
671  // use proj4 to do the transform
672 
673  // if the source/destination projection is lat/long, convert the points to radians
674  // prior to transforming
675  ProjData projData = d->threadLocalProjData();
676 
677  int projResult = 0;
678 
679  proj_errno_reset( projData );
680  proj_trans_generic( projData, ( direction == ForwardTransform && !d->mIsReversed ) || ( direction == ReverseTransform && d->mIsReversed ) ? PJ_FWD : PJ_INV,
681  x, sizeof( double ), numPoints,
682  y, sizeof( double ), numPoints,
683  z, sizeof( double ), numPoints,
684  useTime ? t.data() : nullptr, sizeof( double ), useTime ? numPoints : 0 );
685  // Try to - approximatively - emulate the behavior of pj_transform()...
686  // In the case of a single point transform, and a transformation error occurs,
687  // pj_transform() would return the errno. In cases of multiple point transform,
688  // it would continue (for non-transient errors, that is pipeline definition
689  // errors) and just set the resulting x,y to infinity. This is in fact a
690  // bit more subtle than that, and I'm not completely sure the logic in
691  // pj_transform() was really sane & fully bullet proof
692  // So here just check proj_errno() for single point transform
693  int actualRes = 0;
694  if ( numPoints == 1 )
695  {
696  projResult = proj_errno( projData );
697  actualRes = projResult;
698  }
699  else
700  {
701  actualRes = proj_errno( projData );
702  }
703  if ( actualRes == 0 )
704  {
705  // proj_errno is sometimes not an accurate method to test for transform failures - so we need to
706  // manually scan for nan values
707  if ( std::any_of( x, x + numPoints, []( double v ) { return std::isinf( v ); } )
708  || std::any_of( y, y + numPoints, []( double v ) { return std::isinf( v ); } )
709  || std::any_of( z, z + numPoints, []( double v ) { return std::isinf( v ); } ) )
710  {
711  actualRes = 1;
712  }
713  }
714 
715  mFallbackOperationOccurred = false;
716  if ( actualRes != 0
717  && ( 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
718  && ( d->mAllowFallbackTransforms || mBallparkTransformsAreAppropriate ) )
719  {
720  // fail #1 -- try with getting proj to auto-pick an appropriate coordinate operation for the points
721  if ( PJ *transform = d->threadLocalFallbackProjData() )
722  {
723  projResult = 0;
724  proj_errno_reset( transform );
725  proj_trans_generic( transform, direction == ForwardTransform ? PJ_FWD : PJ_INV,
726  xprev.data(), sizeof( double ), numPoints,
727  yprev.data(), sizeof( double ), numPoints,
728  zprev.data(), sizeof( double ), numPoints,
729  useTime ? t.data() : nullptr, sizeof( double ), useTime ? numPoints : 0 );
730  // Try to - approximatively - emulate the behavior of pj_transform()...
731  // In the case of a single point transform, and a transformation error occurs,
732  // pj_transform() would return the errno. In cases of multiple point transform,
733  // it would continue (for non-transient errors, that is pipeline definition
734  // errors) and just set the resulting x,y to infinity. This is in fact a
735  // bit more subtle than that, and I'm not completely sure the logic in
736  // pj_transform() was really sane & fully bullet proof
737  // So here just check proj_errno() for single point transform
738  if ( numPoints == 1 )
739  {
740  // hmm - something very odd here. We can't trust proj_errno( transform ), as that's giving us incorrect error numbers
741  // (such as "failed to load datum shift file", which is definitely incorrect for a default proj created operation!)
742  // so we resort to testing values ourselves...
743  projResult = std::isinf( xprev[0] ) || std::isinf( yprev[0] ) || std::isinf( zprev[0] ) ? 1 : 0;
744  }
745 
746  if ( projResult == 0 )
747  {
748  memcpy( x, xprev.data(), sizeof( double ) * numPoints );
749  memcpy( y, yprev.data(), sizeof( double ) * numPoints );
750  memcpy( z, zprev.data(), sizeof( double ) * numPoints );
751  mFallbackOperationOccurred = true;
752  }
753 
754  if ( !mBallparkTransformsAreAppropriate && !mDisableFallbackHandler && sFallbackOperationOccurredHandler )
755  {
756  sFallbackOperationOccurredHandler( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation );
757 #if 0
758  const QString warning = QStringLiteral( "A fallback coordinate operation was used between %1 and %2" ).arg( d->mSourceCRS.authid(),
759  d->mDestCRS.authid() );
760  qWarning( "%s", warning.toLatin1().constData() );
761 #endif
762  }
763  }
764  }
765 
766  for ( const int &pos : zNanPositions )
767  {
768  z[pos] = std::numeric_limits<double>::quiet_NaN();
769  }
770 
771  if ( projResult != 0 )
772  {
773  //something bad happened....
774  QString points;
775 
776  for ( int i = 0; i < numPoints; ++i )
777  {
778  if ( direction == ForwardTransform )
779  {
780  points += QStringLiteral( "(%1, %2)\n" ).arg( x[i], 0, 'f' ).arg( y[i], 0, 'f' );
781  }
782  else
783  {
784  points += QStringLiteral( "(%1, %2)\n" ).arg( x[i], 0, 'f' ).arg( y[i], 0, 'f' );
785  }
786  }
787 
788  const QString dir = ( direction == ForwardTransform ) ? QObject::tr( "forward transform" ) : QObject::tr( "inverse transform" );
789 
790  const QString msg = QObject::tr( "%1 of\n"
791  "%2"
792  "Error: %3" )
793  .arg( dir,
794  points,
795  projResult < 0 ? QString::fromUtf8( proj_errno_string( projResult ) ) : QObject::tr( "Fallback transform failed" ) );
796 
797 
798  // don't flood console with thousands of duplicate transform error messages
799  if ( msg != mLastError )
800  {
801  QgsDebugMsg( "Projection failed emitting invalid transform signal: " + msg );
802  mLastError = msg;
803  }
804  QgsDebugMsgLevel( QStringLiteral( "rethrowing exception" ), 2 );
805 
806  throw QgsCsException( msg );
807  }
808 
809 #ifdef COORDINATE_TRANSFORM_VERBOSE
810  QgsDebugMsg( QStringLiteral( "[[[[[[ Projected %1, %2 to %3, %4 ]]]]]]" )
811  .arg( xorg, 0, 'g', 15 ).arg( yorg, 0, 'g', 15 )
812  .arg( *x, 0, 'g', 15 ).arg( *y, 0, 'g', 15 ) );
813 #endif
814 }
815 
817 {
818  return d->mIsValid;
819 }
820 
822 {
823  return !d->mIsValid || d->mShortCircuit;
824 }
825 
827 {
828  return d->mProjCoordinateOperation;
829 }
830 
832 {
833  ProjData projData = d->threadLocalProjData();
835 }
836 
837 void QgsCoordinateTransform::setCoordinateOperation( const QString &operation ) const
838 {
839  d.detach();
840  d->mProjCoordinateOperation = operation;
841  d->mShouldReverseCoordinateOperation = false;
842 }
843 
845 {
846  d.detach();
847  d->mAllowFallbackTransforms = allowed;
848 }
849 
851 {
852  return d->mAllowFallbackTransforms;
853 }
854 
856 {
857  mBallparkTransformsAreAppropriate = appropriate;
858 }
859 
861 {
862  mDisableFallbackHandler = disabled;
863 }
864 
866 {
867  return mFallbackOperationOccurred;
868 }
869 
870 const char *finder( const char *name )
871 {
872  QString proj;
873 #ifdef Q_OS_WIN
874  proj = QApplication::applicationDirPath()
875  + "/share/proj/" + QString( name );
876 #else
877  Q_UNUSED( name )
878 #endif
879  return proj.toUtf8();
880 }
881 
882 bool QgsCoordinateTransform::setFromCache( const QgsCoordinateReferenceSystem &src, const QgsCoordinateReferenceSystem &dest, const QString &coordinateOperationProj, bool allowFallback )
883 {
884  if ( !src.isValid() || !dest.isValid() )
885  return false;
886 
887  const QString sourceKey = src.authid().isEmpty() ?
889  const QString destKey = dest.authid().isEmpty() ?
891 
892  if ( sourceKey.isEmpty() || destKey.isEmpty() )
893  return false;
894 
895  QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Read );
896  if ( sDisableCache )
897  return false;
898 
899  const QList< QgsCoordinateTransform > values = sTransforms.values( qMakePair( sourceKey, destKey ) );
900  for ( auto valIt = values.constBegin(); valIt != values.constEnd(); ++valIt )
901  {
902  if ( ( *valIt ).coordinateOperation() == coordinateOperationProj
903  && ( *valIt ).allowFallbackTransforms() == allowFallback
904  && qgsNanCompatibleEquals( src.coordinateEpoch(), ( *valIt ).sourceCrs().coordinateEpoch() )
905  && qgsNanCompatibleEquals( dest.coordinateEpoch(), ( *valIt ).destinationCrs().coordinateEpoch() )
906  )
907  {
908  // need to save, and then restore the context... we don't want this to be cached or to use the values from the cache
909  const QgsCoordinateTransformContext context = mContext;
910 #ifdef QGISDEBUG
911  const bool hasContext = mHasContext;
912 #endif
913  *this = *valIt;
914  locker.unlock();
915 
916  mContext = context;
917 #ifdef QGISDEBUG
918  mHasContext = hasContext;
919 #endif
920 
921  return true;
922  }
923  }
924  return false;
925 }
926 
927 void QgsCoordinateTransform::addToCache()
928 {
929  if ( !d->mSourceCRS.isValid() || !d->mDestCRS.isValid() )
930  return;
931 
932  const QString sourceKey = d->mSourceCRS.authid().isEmpty() ?
933  d->mSourceCRS.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED ) : d->mSourceCRS.authid();
934  const QString destKey = d->mDestCRS.authid().isEmpty() ?
935  d->mDestCRS.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED ) : d->mDestCRS.authid();
936 
937  if ( sourceKey.isEmpty() || destKey.isEmpty() )
938  return;
939 
940  const QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Write );
941  if ( sDisableCache )
942  return;
943 
944  sTransforms.insert( qMakePair( sourceKey, destKey ), *this );
945 }
946 
948 {
950  return d->mSourceDatumTransform;
952 }
953 
955 {
956  d.detach();
958  d->mSourceDatumTransform = dt;
960 }
961 
963 {
965  return d->mDestinationDatumTransform;
967 }
968 
970 {
971  d.detach();
973  d->mDestinationDatumTransform = dt;
975 }
976 
978 {
979  const QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Write );
980  if ( sDisableCache )
981  return;
982 
983  if ( disableCache )
984  {
985  sDisableCache = true;
986  }
987 
988  sTransforms.clear();
989 }
990 
991 void QgsCoordinateTransform::removeFromCacheObjectsBelongingToCurrentThread( void *pj_context )
992 {
993  // Not completely sure about object order destruction after main() has
994  // exited. So it is safer to check sDisableCache before using sCacheLock
995  // in case sCacheLock would have been destroyed before the current TLS
996  // QgsProjContext object that has called us...
997  if ( sDisableCache )
998  return;
999 
1000  const QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Write );
1001  // cppcheck-suppress identicalConditionAfterEarlyExit
1002  if ( sDisableCache )
1003  return;
1004 
1005  for ( auto it = sTransforms.begin(); it != sTransforms.end(); )
1006  {
1007  auto &v = it.value();
1008  if ( v.d->removeObjectsBelongingToCurrentThread( pj_context ) )
1009  it = sTransforms.erase( it );
1010  else
1011  ++it;
1012  }
1013 }
1014 
1015 double QgsCoordinateTransform::scaleFactor( const QgsRectangle &ReferenceExtent ) const
1016 {
1017  const QgsPointXY source1( ReferenceExtent.xMinimum(), ReferenceExtent.yMinimum() );
1018  const QgsPointXY source2( ReferenceExtent.xMaximum(), ReferenceExtent.yMaximum() );
1019  const double distSourceUnits = std::sqrt( source1.sqrDist( source2 ) );
1020  const QgsPointXY dest1 = transform( source1 );
1021  const QgsPointXY dest2 = transform( source2 );
1022  const double distDestUnits = std::sqrt( dest1.sqrDist( dest2 ) );
1023  return distDestUnits / distSourceUnits;
1024 }
1025 
1027 {
1028  QgsCoordinateTransformPrivate::setCustomMissingRequiredGridHandler( handler );
1029 }
1030 
1032 {
1033  QgsCoordinateTransformPrivate::setCustomMissingPreferredGridHandler( handler );
1034 }
1035 
1037 {
1038  QgsCoordinateTransformPrivate::setCustomCoordinateOperationCreationErrorHandler( handler );
1039 }
1040 
1042 {
1043  QgsCoordinateTransformPrivate::setCustomMissingGridUsedByContextHandler( handler );
1044 }
1045 
1046 void QgsCoordinateTransform::setFallbackOperationOccurredHandler( const std::function<void ( const QgsCoordinateReferenceSystem &, const QgsCoordinateReferenceSystem &, const QString & )> &handler )
1047 {
1048  sFallbackOperationOccurredHandler = handler;
1049 }
1050 
1052 {
1053  QgsCoordinateTransformPrivate::setDynamicCrsToDynamicCrsWarningHandler( handler );
1054 }
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.
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...
QgsPointXY transform(const QgsPointXY &point, TransformDirection direction=ForwardTransform) const SIP_THROW(QgsCsException)
Transform the point from the source CRS to the destination CRS.
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.
void transformPolygon(QPolygonF &polygon, TransformDirection direction=ForwardTransform) const SIP_THROW(QgsCsException)
Transforms a polygon to the destination coordinate system.
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.
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.
TransformDirection
Enum used to indicate the direction (forward or inverse) of the transform.
@ ForwardTransform
Transform from source to destination CRS.
@ ReverseTransform
Transform from destination to source CRS.
QgsRectangle transformBoundingBox(const QgsRectangle &rectangle, TransformDirection direction=ForwardTransform, bool handle180Crossover=false) const SIP_THROW(QgsCsException)
Transforms a rectangle from the source CRS to the destination CRS.
void transformCoords(int numPoint, double *x, double *y, double *z, TransformDirection direction=ForwardTransform) const SIP_THROW(QgsCsException)
Transform an array of coordinates to the destination CRS.
QgsCoordinateTransform & operator=(const QgsCoordinateTransform &o)
Assignment operator.
void transformInPlace(double &x, double &y, double &z, TransformDirection direction=ForwardTransform) const SIP_THROW(QgsCsException)
Transforms an array of x, y and z double coordinates in place, from the source CRS to the destination...
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.
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:1604
#define Q_NOWARN_DEPRECATED_PUSH
Definition: qgis.h:1603
bool qgsNanCompatibleEquals(double a, double b)
Compare two doubles, treating nan values as equal.
Definition: qgis.h:1092
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.