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