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