QGIS API Documentation  3.9.0-Master (224899f119)
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 = static_cast< int >( std::ceil( rect.width() / d ) ) + 1;
525  int nYPoints = static_cast< int >( std::ceil( rect.height() / d ) ) + 1;
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 
533  QVector<double> x( nXPoints * nYPoints );
534  QVector<double> y( nXPoints * nYPoints );
535  QVector<double> z( nXPoints * nYPoints );
536 
537  QgsDebugMsgLevel( QStringLiteral( "Entering transformBoundingBox..." ), 4 );
538 
539  // Populate the vectors
540 
541  double dx = rect.width() / static_cast< double >( nXPoints - 1 );
542  double dy = rect.height() / static_cast< double >( nYPoints - 1 );
543 
544  double pointY = rect.yMinimum();
545 
546  for ( int i = 0; i < nYPoints ; i++ )
547  {
548 
549  // Start at right edge
550  double pointX = rect.xMinimum();
551 
552  for ( int j = 0; j < nXPoints; j++ )
553  {
554  x[( i * nXPoints ) + j] = pointX;
555  y[( i * nXPoints ) + j] = pointY;
556  // and the height...
557  z[( i * nXPoints ) + j] = 0.0;
558  // QgsDebugMsg(QString("BBox coord: (%1, %2)").arg(x[(i*numP) + j]).arg(y[(i*numP) + j]));
559  pointX += dx;
560  }
561  pointY += dy;
562  }
563 
564  // Do transformation. Any exception generated must
565  // be handled in above layers.
566  try
567  {
568  transformCoords( nXPoints * nYPoints, x.data(), y.data(), z.data(), direction );
569  }
570  catch ( const QgsCsException & )
571  {
572  // rethrow the exception
573  QgsDebugMsg( QStringLiteral( "rethrowing exception" ) );
574  throw;
575  }
576 
577  // Calculate the bounding box and use that for the extent
578 
579  for ( int i = 0; i < nXPoints * nYPoints; i++ )
580  {
581  if ( !std::isfinite( x[i] ) || !std::isfinite( y[i] ) )
582  {
583  continue;
584  }
585 
586  if ( handle180Crossover )
587  {
588  //if crossing the date line, temporarily add 360 degrees to -ve longitudes
589  bb_rect.combineExtentWith( x[i] >= 0.0 ? x[i] : x[i] + 360.0, y[i] );
590  }
591  else
592  {
593  bb_rect.combineExtentWith( x[i], y[i] );
594  }
595  }
596 
597  if ( bb_rect.isNull() )
598  {
599  // something bad happened when reprojecting the filter rect... no finite points were left!
600  throw QgsCsException( QObject::tr( "Could not transform bounding box to target CRS" ) );
601  }
602 
603  if ( handle180Crossover )
604  {
605  //subtract temporary addition of 360 degrees from longitudes
606  if ( bb_rect.xMinimum() > 180.0 )
607  bb_rect.setXMinimum( bb_rect.xMinimum() - 360.0 );
608  if ( bb_rect.xMaximum() > 180.0 )
609  bb_rect.setXMaximum( bb_rect.xMaximum() - 360.0 );
610  }
611 
612  QgsDebugMsgLevel( "Projected extent: " + bb_rect.toString(), 4 );
613 
614  if ( bb_rect.isEmpty() )
615  {
616  QgsDebugMsgLevel( "Original extent: " + rect.toString(), 4 );
617  }
618 
619  return bb_rect;
620 }
621 
622 void QgsCoordinateTransform::transformCoords( int numPoints, double *x, double *y, double *z, TransformDirection direction ) const
623 {
624  if ( !d->mIsValid || d->mShortCircuit )
625  return;
626  // Refuse to transform the points if the srs's are invalid
627  if ( !d->mSourceCRS.isValid() )
628  {
629  QgsMessageLog::logMessage( QObject::tr( "The source spatial reference system (CRS) is not valid. "
630  "The coordinates can not be reprojected. The CRS is: %1" )
631  .arg( d->mSourceCRS.toProj4() ), QObject::tr( "CRS" ) );
632  return;
633  }
634  if ( !d->mDestCRS.isValid() )
635  {
636  QgsMessageLog::logMessage( QObject::tr( "The destination spatial reference system (CRS) is not valid. "
637  "The coordinates can not be reprojected. The CRS is: %1" ).arg( d->mDestCRS.toProj4() ), QObject::tr( "CRS" ) );
638  return;
639  }
640 
641 #ifdef COORDINATE_TRANSFORM_VERBOSE
642  double xorg = *x;
643  double yorg = *y;
644  QgsDebugMsg( QStringLiteral( "[[[[[[ Number of points to transform: %1 ]]]]]]" ).arg( numPoints ) );
645 #endif
646 
647 #ifdef QGISDEBUG
648  if ( !mHasContext )
649  QgsDebugMsgLevel( QStringLiteral( "No QgsCoordinateTransformContext context set for transform" ), 4 );
650 #endif
651 
652  // use proj4 to do the transform
653 
654  // if the source/destination projection is lat/long, convert the points to radians
655  // prior to transforming
656  ProjData projData = d->threadLocalProjData();
657 
658  int projResult = 0;
659 #if PROJ_VERSION_MAJOR>=6
660  proj_errno_reset( projData );
661  proj_trans_generic( projData, direction == ForwardTransform ? PJ_FWD : PJ_INV,
662  x, sizeof( double ), numPoints,
663  y, sizeof( double ), numPoints,
664  z, sizeof( double ), numPoints,
665  nullptr, sizeof( double ), 0 );
666  projResult = proj_errno( projData );
667 #else
668  bool sourceIsLatLong = false;
669  bool destIsLatLong = false;
670 
671  projPJ sourceProj = projData.first;
672  projPJ destProj = projData.second;
673  sourceIsLatLong = pj_is_latlong( sourceProj );
674  destIsLatLong = pj_is_latlong( destProj );
675 
676  if ( ( destIsLatLong && ( direction == ReverseTransform ) )
677  || ( sourceIsLatLong && ( direction == ForwardTransform ) ) )
678  {
679  for ( int i = 0; i < numPoints; ++i )
680  {
681  x[i] *= DEG_TO_RAD;
682  y[i] *= DEG_TO_RAD;
683  }
684  }
685 #endif
686 
687 #if PROJ_VERSION_MAJOR<6
688  if ( direction == ReverseTransform )
689  {
690  projResult = pj_transform( destProj, sourceProj, numPoints, 0, x, y, z );
691  }
692  else
693  {
694  Q_ASSERT( sourceProj );
695  Q_ASSERT( destProj );
696  projResult = pj_transform( sourceProj, destProj, numPoints, 0, x, y, z );
697  }
698 #endif
699 
700  if ( projResult != 0 )
701  {
702  //something bad happened....
703  QString points;
704 
705  for ( int i = 0; i < numPoints; ++i )
706  {
707  if ( direction == ForwardTransform )
708  {
709  points += QStringLiteral( "(%1, %2)\n" ).arg( x[i], 0, 'f' ).arg( y[i], 0, 'f' );
710  }
711  else
712  {
713 #if PROJ_VERSION_MAJOR>=6
714  points += QStringLiteral( "(%1, %2)\n" ).arg( x[i], 0, 'f' ).arg( y[i], 0, 'f' );
715 #else
716  points += QStringLiteral( "(%1, %2)\n" ).arg( x[i] * RAD_TO_DEG, 0, 'f' ).arg( y[i] * RAD_TO_DEG, 0, 'f' );
717 #endif
718  }
719  }
720 
721  QString dir = ( direction == ForwardTransform ) ? QObject::tr( "forward transform" ) : QObject::tr( "inverse transform" );
722 
723 #if PROJ_VERSION_MAJOR>=6
724  QgsProjUtils::proj_pj_unique_ptr src( proj_get_source_crs( QgsProjContext::get(), projData ) );
725  QgsProjUtils::proj_pj_unique_ptr dest( proj_get_source_crs( QgsProjContext::get(), projData ) );
726  QString msg = QObject::tr( "%1 of\n"
727  "%2"
728  "PROJ: %3\n"
729  "Error: %4" )
730  .arg( dir,
731  points,
732  proj_as_proj_string( QgsProjContext::get(), projData, PJ_PROJ_5, nullptr ),
733  QString::fromUtf8( proj_errno_string( projResult ) ) );
734 #else
735  char *srcdef = pj_get_def( sourceProj, 0 );
736  char *dstdef = pj_get_def( destProj, 0 );
737 
738  QString msg = QObject::tr( "%1 of\n"
739  "%2"
740  "PROJ: %3 +to %4\n"
741  "Error: %5" )
742  .arg( dir,
743  points,
744  srcdef, dstdef,
745  QString::fromUtf8( pj_strerrno( projResult ) ) );
746 
747  pj_dalloc( srcdef );
748  pj_dalloc( dstdef );
749 #endif
750 
751  QgsDebugMsg( "Projection failed emitting invalid transform signal: " + msg );
752  QgsDebugMsg( QStringLiteral( "throwing exception" ) );
753 
754  throw QgsCsException( msg );
755  }
756 
757 #if PROJ_VERSION_MAJOR<6
758  // if the result is lat/long, convert the results from radians back
759  // to degrees
760  if ( ( destIsLatLong && ( direction == ForwardTransform ) )
761  || ( sourceIsLatLong && ( direction == ReverseTransform ) ) )
762  {
763  for ( int i = 0; i < numPoints; ++i )
764  {
765  x[i] *= RAD_TO_DEG;
766  y[i] *= RAD_TO_DEG;
767  }
768  }
769 #endif
770 #ifdef COORDINATE_TRANSFORM_VERBOSE
771  QgsDebugMsg( QStringLiteral( "[[[[[[ Projected %1, %2 to %3, %4 ]]]]]]" )
772  .arg( xorg, 0, 'g', 15 ).arg( yorg, 0, 'g', 15 )
773  .arg( *x, 0, 'g', 15 ).arg( *y, 0, 'g', 15 ) );
774 #endif
775 }
776 
778 {
779  return d->mIsValid;
780 }
781 
783 {
784  return !d->mIsValid || d->mShortCircuit;
785 }
786 
788 {
789  return d->mProjCoordinateOperation;
790 }
791 
792 void QgsCoordinateTransform::setCoordinateOperation( const QString &operation ) const
793 {
794  d.detach();
795  d->mProjCoordinateOperation = operation;
796 }
797 
798 const char *finder( const char *name )
799 {
800  QString proj;
801 #ifdef Q_OS_WIN
802  proj = QApplication::applicationDirPath()
803  + "/share/proj/" + QString( name );
804 #else
805  Q_UNUSED( name )
806 #endif
807  return proj.toUtf8();
808 }
809 
810 #if PROJ_VERSION_MAJOR>=6
811 bool QgsCoordinateTransform::setFromCache( const QgsCoordinateReferenceSystem &src, const QgsCoordinateReferenceSystem &dest, const QString &coordinateOperationProj )
812 {
813  if ( !src.isValid() || !dest.isValid() )
814  return false;
815 
816  const QString sourceKey = src.authid().isEmpty() ?
817  src.toWkt() : src.authid();
818  const QString destKey = dest.authid().isEmpty() ?
819  dest.toWkt() : dest.authid();
820 
821  if ( sourceKey.isEmpty() || destKey.isEmpty() )
822  return false;
823 
824  QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Read );
825  if ( sDisableCache )
826  return false;
827 
828  const QList< QgsCoordinateTransform > values = sTransforms.values( qMakePair( sourceKey, destKey ) );
829  for ( auto valIt = values.constBegin(); valIt != values.constEnd(); ++valIt )
830  {
831  if ( ( *valIt ).coordinateOperation() == coordinateOperationProj )
832  {
833  // need to save, and then restore the context... we don't want this to be cached or to use the values from the cache
835 #ifdef QGISDEBUG
836  bool hasContext = mHasContext;
837 #endif
838  *this = *valIt;
839  locker.unlock();
840 
841  mContext = context;
842 #ifdef QGISDEBUG
843  mHasContext = hasContext;
844 #endif
845 
846  return true;
847  }
848  }
849  return false;
850 }
851 #else
852 bool QgsCoordinateTransform::setFromCache( const QgsCoordinateReferenceSystem &src, const QgsCoordinateReferenceSystem &dest, int srcDatumTransform, int destDatumTransform )
853 {
854  if ( !src.isValid() || !dest.isValid() )
855  return false;
856 
857  const QString sourceKey = src.authid().isEmpty() ?
858  src.toWkt() : src.authid();
859  const QString destKey = dest.authid().isEmpty() ?
860  dest.toWkt() : dest.authid();
861 
862  if ( sourceKey.isEmpty() || destKey.isEmpty() )
863  return false;
864 
865  QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Read );
866  if ( sDisableCache )
867  return false;
868 
869  const QList< QgsCoordinateTransform > values = sTransforms.values( qMakePair( src.authid(), dest.authid() ) );
870  for ( auto valIt = values.constBegin(); valIt != values.constEnd(); ++valIt )
871  {
873  if ( ( *valIt ).sourceDatumTransformId() == srcDatumTransform &&
874  ( *valIt ).destinationDatumTransformId() == destDatumTransform )
875  {
876  // need to save, and then restore the context... we don't want this to be cached or to use the values from the cache
878 #ifdef QGISDEBUG
879  bool hasContext = mHasContext;
880 #endif
881  *this = *valIt;
882  locker.unlock();
883 
884  mContext = context;
885 #ifdef QGISDEBUG
886  mHasContext = hasContext;
887 #endif
888 
889  return true;
890  }
892  }
893  return false;
894 }
895 #endif
896 
897 void QgsCoordinateTransform::addToCache()
898 {
899  if ( !d->mSourceCRS.isValid() || !d->mDestCRS.isValid() )
900  return;
901 
902  const QString sourceKey = d->mSourceCRS.authid().isEmpty() ?
903  d->mSourceCRS.toWkt() : d->mSourceCRS.authid();
904  const QString destKey = d->mDestCRS.authid().isEmpty() ?
905  d->mDestCRS.toWkt() : d->mDestCRS.authid();
906 
907  if ( sourceKey.isEmpty() || destKey.isEmpty() )
908  return;
909 
910  QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Write );
911  if ( sDisableCache )
912  return;
913 
914  sTransforms.insertMulti( qMakePair( sourceKey, destKey ), *this );
915 }
916 
918 {
920  return d->mSourceDatumTransform;
922 }
923 
925 {
926  d.detach();
928  d->mSourceDatumTransform = dt;
930 }
931 
933 {
935  return d->mDestinationDatumTransform;
937 }
938 
940 {
941  d.detach();
943  d->mDestinationDatumTransform = dt;
945 }
946 
948 {
949  QgsReadWriteLocker locker( sCacheLock, QgsReadWriteLocker::Write );
950  if ( sDisableCache )
951  return;
952 
953  if ( disableCache )
954  {
955  sDisableCache = true;
956  }
957 
958  sTransforms.clear();
959 }
960 
961 double QgsCoordinateTransform::scaleFactor( const QgsRectangle &ReferenceExtent ) const
962 {
963  QgsPointXY source1( ReferenceExtent.xMinimum(), ReferenceExtent.yMinimum() );
964  QgsPointXY source2( ReferenceExtent.xMaximum(), ReferenceExtent.yMaximum() );
965  double distSourceUnits = std::sqrt( source1.sqrDist( source2 ) );
966  QgsPointXY dest1 = transform( source1 );
967  QgsPointXY dest2 = transform( source2 );
968  double distDestUnits = std::sqrt( dest1.sqrDist( dest2 ) );
969  return distDestUnits / distSourceUnits;
970 }
971 
973 {
974  QgsCoordinateTransformPrivate::setCustomMissingRequiredGridHandler( handler );
975 }
976 
978 {
979  QgsCoordinateTransformPrivate::setCustomMissingPreferredGridHandler( handler );
980 }
981 
983 {
984  QgsCoordinateTransformPrivate::setCustomCoordinateOperationCreationErrorHandler( handler );
985 }
986 
988 {
989  QgsCoordinateTransformPrivate::setCustomMissingGridUsedByContextHandler( handler );
990 }
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 setXMaximum(double x)
Set the maximum x value.
Definition: qgsrectangle.h:135
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:634
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...
Reads and writes project states.
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
void combineExtentWith(const QgsRectangle &rect)
Expands the rectangle so that it covers both the original rectangle and the given rectangle...
Definition: qgsrectangle.h:359
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:635
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:
bool isNull() const
Test if the rectangle is null (all coordinates zero or after call to setMinimal()).
Definition: qgsrectangle.h:436
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.
void setXMinimum(double x)
Set the minimum x value.
Definition: qgsrectangle.h:130
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.