QGIS API Documentation  3.11.0-Master (358e466517)
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 
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 
53 {
54  d = new QgsCoordinateTransformPrivate();
55 }
56 
58 {
59  mContext = context;
60  d = new QgsCoordinateTransformPrivate( source, destination, mContext );
61 #ifdef QGISDEBUG
62  mHasContext = true;
63 #endif
64 
65  if ( !d->checkValidity() )
66  return;
67 
69 #if PROJ_VERSION_MAJOR>=6
70  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation ) )
71 #else
72  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mSourceDatumTransform, d->mDestinationDatumTransform ) )
73 #endif
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 PROJ_VERSION_MAJOR>=6
95  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation ) )
96 #else
97  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mSourceDatumTransform, d->mDestinationDatumTransform ) )
98 #endif
99  {
100  d->initialize();
101  addToCache();
102  }
104 }
105 
106 QgsCoordinateTransform::QgsCoordinateTransform( const QgsCoordinateReferenceSystem &source, const QgsCoordinateReferenceSystem &destination, int sourceDatumTransform, int destinationDatumTransform )
107 {
108  d = new QgsCoordinateTransformPrivate( source, destination, sourceDatumTransform, destinationDatumTransform );
109 #ifdef QGISDEBUG
110  mHasContext = true; // not strictly true, but we don't need to worry if datums have been explicitly set
111 #endif
112 
113  if ( !d->checkValidity() )
114  return;
115 
117 #if PROJ_VERSION_MAJOR>=6
118  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation ) )
119 #else
120  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mSourceDatumTransform, d->mDestinationDatumTransform ) )
121 #endif
122  {
123  d->initialize();
124  addToCache();
125  }
127 }
128 
130  : mContext( o.mContext )
131 #ifdef QGISDEBUG
132  , mHasContext( o.mHasContext )
133 #endif
134 {
135  d = o.d;
136 }
137 
139 {
140  d = o.d;
141 #ifdef QGISDEBUG
142  mHasContext = o.mHasContext;
143 #endif
144  mContext = o.mContext;
145  return *this;
146 }
147 
149 
151 {
152  d.detach();
153  d->mSourceCRS = crs;
154  if ( !d->checkValidity() )
155  return;
156 
157  d->calculateTransforms( mContext );
159 #if PROJ_VERSION_MAJOR>=6
160  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation ) )
161 #else
162  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mSourceDatumTransform, d->mDestinationDatumTransform ) )
163 #endif
164  {
165  d->initialize();
166  addToCache();
167  }
169 }
171 {
172  d.detach();
173  d->mDestCRS = crs;
174  if ( !d->checkValidity() )
175  return;
176 
177  d->calculateTransforms( mContext );
179 #if PROJ_VERSION_MAJOR>=6
180  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation ) )
181 #else
182  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mSourceDatumTransform, d->mDestinationDatumTransform ) )
183 #endif
184  {
185  d->initialize();
186  addToCache();
187  }
189 }
190 
192 {
193  d.detach();
194  mContext = context;
195 #ifdef QGISDEBUG
196  mHasContext = true;
197 #endif
198  if ( !d->checkValidity() )
199  return;
200 
201  d->calculateTransforms( mContext );
203 #if PROJ_VERSION_MAJOR>=6
204  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mProjCoordinateOperation ) )
205 #else
206  if ( !setFromCache( d->mSourceCRS, d->mDestCRS, d->mSourceDatumTransform, d->mDestinationDatumTransform ) )
207 #endif
208  {
209  d->initialize();
210  addToCache();
211  }
213 }
214 
216 {
217  return mContext;
218 }
219 
221 {
222  return d->mSourceCRS;
223 }
224 
226 {
227  return d->mDestCRS;
228 }
229 
231 {
232  if ( !d->mIsValid || d->mShortCircuit )
233  return point;
234 
235  // transform x
236  double x = point.x();
237  double y = point.y();
238  double z = 0.0;
239  try
240  {
241  transformCoords( 1, &x, &y, &z, direction );
242  }
243  catch ( const QgsCsException & )
244  {
245  // rethrow the exception
246  QgsDebugMsg( QStringLiteral( "rethrowing exception" ) );
247  throw;
248  }
249 
250  return QgsPointXY( x, y );
251 }
252 
253 
254 QgsPointXY QgsCoordinateTransform::transform( const double theX, const double theY = 0.0, TransformDirection direction ) const
255 {
256  try
257  {
258  return transform( QgsPointXY( theX, theY ), direction );
259  }
260  catch ( const QgsCsException & )
261  {
262  // rethrow the exception
263  QgsDebugMsg( QStringLiteral( "rethrowing exception" ) );
264  throw;
265  }
266 }
267 
269 {
270  if ( !d->mIsValid || d->mShortCircuit )
271  return rect;
272  // transform x
273  double x1 = rect.xMinimum();
274  double y1 = rect.yMinimum();
275  double x2 = rect.xMaximum();
276  double y2 = rect.yMaximum();
277 
278  // Number of points to reproject------+
279  // |
280  // V
281  try
282  {
283  double z = 0.0;
284  transformCoords( 1, &x1, &y1, &z, direction );
285  transformCoords( 1, &x2, &y2, &z, direction );
286  }
287  catch ( const QgsCsException & )
288  {
289  // rethrow the exception
290  QgsDebugMsg( QStringLiteral( "rethrowing exception" ) );
291  throw;
292  }
293 
294 #ifdef COORDINATE_TRANSFORM_VERBOSE
295  QgsDebugMsg( QStringLiteral( "Rect projection..." ) );
296  QgsDebugMsg( QStringLiteral( "Xmin : %1 --> %2" ).arg( rect.xMinimum() ).arg( x1 ) );
297  QgsDebugMsg( QStringLiteral( "Ymin : %1 --> %2" ).arg( rect.yMinimum() ).arg( y1 ) );
298  QgsDebugMsg( QStringLiteral( "Xmax : %1 --> %2" ).arg( rect.xMaximum() ).arg( x2 ) );
299  QgsDebugMsg( QStringLiteral( "Ymax : %1 --> %2" ).arg( rect.yMaximum() ).arg( y2 ) );
300 #endif
301  return QgsRectangle( x1, y1, x2, y2 );
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  QgsDebugMsg( QStringLiteral( "rethrowing exception" ) );
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  QgsDebugMsg( QStringLiteral( "rethrowing exception" ) );
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  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  QgsDebugMsg( QStringLiteral( "rethrowing exception" ) );
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  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  QgsDebugMsg( QStringLiteral( "rethrowing exception" ) );
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  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  double d = std::sqrt( ( rect.width() * rect.height() ) / std::pow( std::sqrt( static_cast< double >( nPoints ) ) - 1, 2.0 ) );
524  int nXPoints = std::min( static_cast< int >( std::ceil( rect.width() / d ) ) + 1, 1000 );
525  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  double dx = rect.width() / static_cast< double >( nXPoints - 1 );
541  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  QgsDebugMsg( QStringLiteral( "rethrowing exception" ) );
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.toProj4() ), 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.toProj4() ), QObject::tr( "CRS" ) );
637  return;
638  }
639 
640 #ifdef COORDINATE_TRANSFORM_VERBOSE
641  double xorg = *x;
642  double yorg = *y;
643  QgsDebugMsg( QStringLiteral( "[[[[[[ Number of points to transform: %1 ]]]]]]" ).arg( numPoints ) );
644 #endif
645 
646 #ifdef QGISDEBUG
647  if ( !mHasContext )
648  QgsDebugMsgLevel( QStringLiteral( "No QgsCoordinateTransformContext context set for transform" ), 4 );
649 #endif
650 
651  // use proj4 to do the transform
652 
653  // if the source/destination projection is lat/long, convert the points to radians
654  // prior to transforming
655  ProjData projData = d->threadLocalProjData();
656 
657  int projResult = 0;
658 #if PROJ_VERSION_MAJOR>=6
659  proj_errno_reset( projData );
660  proj_trans_generic( projData, direction == ForwardTransform ? PJ_FWD : PJ_INV,
661  x, sizeof( double ), numPoints,
662  y, sizeof( double ), numPoints,
663  z, sizeof( double ), numPoints,
664  nullptr, sizeof( double ), 0 );
665  // Try to - approximatively - emulate the behavior of pj_transform()...
666  // In the case of a single point transform, and a transformation error occurs,
667  // pj_transform() would return the errno. In cases of multiple point transform,
668  // it would continue (for non-transient errors, that is pipeline definition
669  // errors) and just set the resulting x,y to infinity. This is in fact a
670  // bit more subtle than that, and I'm not completely sure the logic in
671  // pj_transform() was really sane & fully bullet proof
672  // So here just check proj_errno() for single point transform
673  if ( numPoints == 1 )
674  {
675  projResult = proj_errno( projData );
676  }
677 #else
678  bool sourceIsLatLong = false;
679  bool destIsLatLong = false;
680 
681  projPJ sourceProj = projData.first;
682  projPJ destProj = projData.second;
683  sourceIsLatLong = pj_is_latlong( sourceProj );
684  destIsLatLong = pj_is_latlong( destProj );
685 
686  if ( ( destIsLatLong && ( direction == ReverseTransform ) )
687  || ( sourceIsLatLong && ( direction == ForwardTransform ) ) )
688  {
689  for ( int i = 0; i < numPoints; ++i )
690  {
691  x[i] *= DEG_TO_RAD;
692  y[i] *= DEG_TO_RAD;
693  }
694  }
695 #endif
696 
697 #if PROJ_VERSION_MAJOR<6
698  if ( direction == ReverseTransform )
699  {
700  projResult = pj_transform( destProj, sourceProj, numPoints, 0, x, y, z );
701  }
702  else
703  {
704  Q_ASSERT( sourceProj );
705  Q_ASSERT( destProj );
706  projResult = pj_transform( sourceProj, destProj, numPoints, 0, x, y, z );
707  }
708 #endif
709 
710  if ( projResult != 0 )
711  {
712  //something bad happened....
713  QString points;
714 
715  for ( int i = 0; i < numPoints; ++i )
716  {
717  if ( direction == ForwardTransform )
718  {
719  points += QStringLiteral( "(%1, %2)\n" ).arg( x[i], 0, 'f' ).arg( y[i], 0, 'f' );
720  }
721  else
722  {
723 #if PROJ_VERSION_MAJOR>=6
724  points += QStringLiteral( "(%1, %2)\n" ).arg( x[i], 0, 'f' ).arg( y[i], 0, 'f' );
725 #else
726  points += QStringLiteral( "(%1, %2)\n" ).arg( x[i] * RAD_TO_DEG, 0, 'f' ).arg( y[i] * RAD_TO_DEG, 0, 'f' );
727 #endif
728  }
729  }
730 
731  QString dir = ( direction == ForwardTransform ) ? QObject::tr( "forward transform" ) : QObject::tr( "inverse transform" );
732 
733 #if PROJ_VERSION_MAJOR>=6
734  QgsProjUtils::proj_pj_unique_ptr src( proj_get_source_crs( QgsProjContext::get(), projData ) );
735  QgsProjUtils::proj_pj_unique_ptr dest( proj_get_source_crs( QgsProjContext::get(), projData ) );
736  QString msg = QObject::tr( "%1 of\n"
737  "%2"
738  "PROJ: %3\n"
739  "Error: %4" )
740  .arg( dir,
741  points,
742  proj_as_proj_string( QgsProjContext::get(), projData, PJ_PROJ_5, nullptr ),
743  QString::fromUtf8( proj_errno_string( projResult ) ) );
744 #else
745  char *srcdef = pj_get_def( sourceProj, 0 );
746  char *dstdef = pj_get_def( destProj, 0 );
747 
748  QString msg = QObject::tr( "%1 of\n"
749  "%2"
750  "PROJ: %3 +to %4\n"
751  "Error: %5" )
752  .arg( dir,
753  points,
754  srcdef, dstdef,
755  QString::fromUtf8( pj_strerrno( projResult ) ) );
756 
757  pj_dalloc( srcdef );
758  pj_dalloc( dstdef );
759 #endif
760 
761  QgsDebugMsg( "Projection failed emitting invalid transform signal: " + msg );
762  QgsDebugMsg( QStringLiteral( "throwing exception" ) );
763 
764  throw QgsCsException( msg );
765  }
766 
767 #if PROJ_VERSION_MAJOR<6
768  // if the result is lat/long, convert the results from radians back
769  // to degrees
770  if ( ( destIsLatLong && ( direction == ForwardTransform ) )
771  || ( sourceIsLatLong && ( direction == ReverseTransform ) ) )
772  {
773  for ( int i = 0; i < numPoints; ++i )
774  {
775  x[i] *= RAD_TO_DEG;
776  y[i] *= RAD_TO_DEG;
777  }
778  }
779 #endif
780 #ifdef COORDINATE_TRANSFORM_VERBOSE
781  QgsDebugMsg( QStringLiteral( "[[[[[[ Projected %1, %2 to %3, %4 ]]]]]]" )
782  .arg( xorg, 0, 'g', 15 ).arg( yorg, 0, 'g', 15 )
783  .arg( *x, 0, 'g', 15 ).arg( *y, 0, 'g', 15 ) );
784 #endif
785 }
786 
788 {
789  return d->mIsValid;
790 }
791 
793 {
794  return !d->mIsValid || d->mShortCircuit;
795 }
796 
798 {
799  return d->mProjCoordinateOperation;
800 }
801 
802 void QgsCoordinateTransform::setCoordinateOperation( const QString &operation ) const
803 {
804  d.detach();
805  d->mProjCoordinateOperation = operation;
806 }
807 
808 const char *finder( const char *name )
809 {
810  QString proj;
811 #ifdef Q_OS_WIN
812  proj = QApplication::applicationDirPath()
813  + "/share/proj/" + QString( name );
814 #else
815  Q_UNUSED( name )
816 #endif
817  return proj.toUtf8();
818 }
819 
820 #if PROJ_VERSION_MAJOR>=6
821 bool QgsCoordinateTransform::setFromCache( const QgsCoordinateReferenceSystem &src, const QgsCoordinateReferenceSystem &dest, const QString &coordinateOperationProj )
822 {
823  if ( !src.isValid() || !dest.isValid() )
824  return false;
825 
826  const QString sourceKey = src.authid().isEmpty() ?
827  src.toWkt() : src.authid();
828  const QString destKey = dest.authid().isEmpty() ?
829  dest.toWkt() : dest.authid();
830 
831  if ( sourceKey.isEmpty() || destKey.isEmpty() )
832  return false;
833 
834  QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Read );
835  if ( sDisableCache )
836  return false;
837 
838  const QList< QgsCoordinateTransform > values = sTransforms.values( qMakePair( sourceKey, destKey ) );
839  for ( auto valIt = values.constBegin(); valIt != values.constEnd(); ++valIt )
840  {
841  if ( ( *valIt ).coordinateOperation() == coordinateOperationProj )
842  {
843  // need to save, and then restore the context... we don't want this to be cached or to use the values from the cache
845 #ifdef QGISDEBUG
846  bool hasContext = mHasContext;
847 #endif
848  *this = *valIt;
849  locker.unlock();
850 
851  mContext = context;
852 #ifdef QGISDEBUG
853  mHasContext = hasContext;
854 #endif
855 
856  return true;
857  }
858  }
859  return false;
860 }
861 #else
862 bool QgsCoordinateTransform::setFromCache( const QgsCoordinateReferenceSystem &src, const QgsCoordinateReferenceSystem &dest, int srcDatumTransform, int destDatumTransform )
863 {
864  if ( !src.isValid() || !dest.isValid() )
865  return false;
866 
867  const QString sourceKey = src.authid().isEmpty() ?
868  src.toWkt() : src.authid();
869  const QString destKey = dest.authid().isEmpty() ?
870  dest.toWkt() : dest.authid();
871 
872  if ( sourceKey.isEmpty() || destKey.isEmpty() )
873  return false;
874 
875  QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Read );
876  if ( sDisableCache )
877  return false;
878 
879  const QList< QgsCoordinateTransform > values = sTransforms.values( qMakePair( src.authid(), dest.authid() ) );
880  for ( auto valIt = values.constBegin(); valIt != values.constEnd(); ++valIt )
881  {
883  if ( ( *valIt ).sourceDatumTransformId() == srcDatumTransform &&
884  ( *valIt ).destinationDatumTransformId() == destDatumTransform )
885  {
886  // need to save, and then restore the context... we don't want this to be cached or to use the values from the cache
888 #ifdef QGISDEBUG
889  bool hasContext = mHasContext;
890 #endif
891  *this = *valIt;
892  locker.unlock();
893 
894  mContext = context;
895 #ifdef QGISDEBUG
896  mHasContext = hasContext;
897 #endif
898 
899  return true;
900  }
902  }
903  return false;
904 }
905 #endif
906 
907 void QgsCoordinateTransform::addToCache()
908 {
909  if ( !d->mSourceCRS.isValid() || !d->mDestCRS.isValid() )
910  return;
911 
912  const QString sourceKey = d->mSourceCRS.authid().isEmpty() ?
913  d->mSourceCRS.toWkt() : d->mSourceCRS.authid();
914  const QString destKey = d->mDestCRS.authid().isEmpty() ?
915  d->mDestCRS.toWkt() : d->mDestCRS.authid();
916 
917  if ( sourceKey.isEmpty() || destKey.isEmpty() )
918  return;
919 
920  QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Write );
921  if ( sDisableCache )
922  return;
923 
924  sTransforms.insertMulti( qMakePair( sourceKey, destKey ), *this );
925 }
926 
928 {
930  return d->mSourceDatumTransform;
932 }
933 
935 {
936  d.detach();
938  d->mSourceDatumTransform = dt;
940 }
941 
943 {
945  return d->mDestinationDatumTransform;
947 }
948 
950 {
951  d.detach();
953  d->mDestinationDatumTransform = dt;
955 }
956 
958 {
959  QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Write );
960  if ( sDisableCache )
961  return;
962 
963  if ( disableCache )
964  {
965  sDisableCache = true;
966  }
967 
968  sTransforms.clear();
969 }
970 
971 #if PROJ_VERSION_MAJOR>=6
972 void QgsCoordinateTransform::removeFromCacheObjectsBelongingToCurrentThread( void *pj_context )
973 {
974  // Not completely sure about object order destruction after main() has
975  // exited. So it is safer to check sDisableCache before using sCacheLock
976  // in case sCacheLock would have been destroyed before the current TLS
977  // QgsProjContext object that has called us...
978  if ( sDisableCache )
979  return;
980 
981  QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Write );
982  if ( sDisableCache )
983  return;
984 
985  for ( auto it = sTransforms.begin(); it != sTransforms.end(); )
986  {
987  auto &v = it.value();
988  if ( v.d->removeObjectsBelongingToCurrentThread( pj_context ) )
989  it = sTransforms.erase( it );
990  else
991  ++it;
992  }
993 }
994 #endif
995 
996 double QgsCoordinateTransform::scaleFactor( const QgsRectangle &ReferenceExtent ) const
997 {
998  QgsPointXY source1( ReferenceExtent.xMinimum(), ReferenceExtent.yMinimum() );
999  QgsPointXY source2( ReferenceExtent.xMaximum(), ReferenceExtent.yMaximum() );
1000  double distSourceUnits = std::sqrt( source1.sqrDist( source2 ) );
1001  QgsPointXY dest1 = transform( source1 );
1002  QgsPointXY dest2 = transform( source2 );
1003  double distDestUnits = std::sqrt( dest1.sqrDist( dest2 ) );
1004  return distDestUnits / distSourceUnits;
1005 }
1006 
1008 {
1009  QgsCoordinateTransformPrivate::setCustomMissingRequiredGridHandler( handler );
1010 }
1011 
1013 {
1014  QgsCoordinateTransformPrivate::setCustomMissingPreferredGridHandler( handler );
1015 }
1016 
1018 {
1019  QgsCoordinateTransformPrivate::setCustomCoordinateOperationCreationErrorHandler( handler );
1020 }
1021 
1023 {
1024  QgsCoordinateTransformPrivate::setCustomMissingGridUsedByContextHandler( handler );
1025 }
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:725
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...
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
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:426
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:90
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...
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:97
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...
#define Q_NOWARN_DEPRECATED_POP
Definition: qgis.h:726
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:
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.
static PJ_CONTEXT * get()
Returns a thread local instance of a proj context, safe for use in the current thread.
Custom exception class for Coordinate Reference System related exceptions.
Definition: qgsexception.h:65
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.