QGIS API Documentation  2.18.21-Las Palmas (9fba24a)
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"
18 #include "qgsapplication.h"
19 #include "qgscrscache.h"
20 #include "qgsmessagelog.h"
21 #include "qgslogger.h"
22 
23 //qt includes
24 #include <QDomNode>
25 #include <QDomElement>
26 #include <QApplication>
27 #include <QPolygonF>
28 #include <QStringList>
29 #include <QVector>
30 
31 extern "C"
32 {
33 #include <proj_api.h>
34 }
35 #include <sqlite3.h>
36 
37 // if defined shows all information about transform to stdout
38 // #define COORDINATE_TRANSFORM_VERBOSE
39 
40 QThreadStorage< QgsCoordinateTransform::QgsProjContextStore* > QgsCoordinateTransform::mProjContext;
41 
43  : QObject()
44  , mShortCircuit( false )
45  , mInitialisedFlag( false )
46  , mSourceDatumTransform( -1 )
47  , mDestinationDatumTransform( -1 )
48 {
49  setFinder();
50 }
51 
53  : QObject()
54  , mShortCircuit( false )
55  , mInitialisedFlag( false )
56  , mSourceDatumTransform( -1 )
57  , mDestinationDatumTransform( -1 )
58 {
59  setFinder();
60  mSourceCRS = source;
61  mDestCRS = dest;
62  initialise();
63 }
64 
65 QgsCoordinateTransform::QgsCoordinateTransform( long theSourceSrsId, long theDestSrsId )
66  : QObject()
67  , mInitialisedFlag( false )
68  , mSourceCRS( QgsCRSCache::instance()->crsBySrsId( theSourceSrsId ) )
69  , mDestCRS( QgsCRSCache::instance()->crsBySrsId( theDestSrsId ) )
70  , mSourceDatumTransform( -1 )
71  , mDestinationDatumTransform( -1 )
72 {
73  initialise();
74 }
75 
76 QgsCoordinateTransform::QgsCoordinateTransform( const QString& theSourceCRS, const QString& theDestCRS )
77  : QObject()
78  , mInitialisedFlag( false )
79  , mSourceDatumTransform( -1 )
80  , mDestinationDatumTransform( -1 )
81 {
82  setFinder();
83  mSourceCRS = QgsCRSCache::instance()->crsByWkt( theSourceCRS );
84  mDestCRS = QgsCRSCache::instance()->crsByWkt( theDestCRS );
85  // initialize the coordinate system data structures
86  //XXX Who spells initialize initialise?
87  //XXX A: Its the queen's english....
88  //XXX : Long live the queen! Lets get on with the initialization...
89  initialise();
90 }
91 
93  const QString& theDestWkt,
94  QgsCoordinateReferenceSystem::CrsType theSourceCRSType )
95  : QObject()
96  , mInitialisedFlag( false )
97  , mSourceDatumTransform( -1 )
98  , mDestinationDatumTransform( -1 )
99 {
100  setFinder();
101 
102  mSourceCRS.createFromId( theSourceSrid, theSourceCRSType );
103  mDestCRS = QgsCRSCache::instance()->crsByWkt( theDestWkt );
104  // initialize the coordinate system data structures
105  //XXX Who spells initialize initialise?
106  //XXX A: Its the queen's english....
107  //XXX : Long live the queen! Lets get on with the initialization...
108  initialise();
109 }
110 
112 {
113  freeProj();
114 }
115 
117 {
121  tr->initialise();
122  return tr;
123 }
124 
126 {
127  mSourceCRS = theCRS;
128  initialise();
129 }
131 {
132  mDestCRS = theCRS;
133  initialise();
134 }
135 
137 {
139  mDestCRS = QgsCRSCache::instance()->crsBySrsId( theCRSID );
140  initialise();
141 }
142 
143 // XXX This whole function is full of multiple return statements!!!
144 // And probably shouldn't be a void
146 {
147  // XXX Warning - multiple return paths in this block!!
148  if ( !mSourceCRS.isValid() )
149  {
150  //mSourceCRS = defaultWkt;
151  // Pass through with no projection since we have no idea what the layer
152  // coordinates are and projecting them may not be appropriate
153  mShortCircuit = true;
154  QgsDebugMsg( "SourceCRS seemed invalid!" );
155  return;
156  }
157 
158  if ( !mDestCRS.isValid() )
159  {
160  //No destination projection is set so we set the default output projection to
161  //be the same as input proj.
162  mDestCRS = QgsCRSCache::instance()->crsByOgcWmsCrs( mSourceCRS.authid() );
163  }
164 
165  bool useDefaultDatumTransform = ( mSourceDatumTransform == - 1 && mDestinationDatumTransform == -1 );
166 
167  freeProj();
168 
169  mSourceProjString = mSourceCRS.toProj4();
170  if ( !useDefaultDatumTransform )
171  {
172  mSourceProjString = stripDatumTransform( mSourceProjString );
173  }
174  if ( mSourceDatumTransform != -1 )
175  {
176  mSourceProjString += ( ' ' + datumTransformString( mSourceDatumTransform ) );
177  }
178 
179  mDestProjString = mDestCRS.toProj4();
180  if ( !useDefaultDatumTransform )
181  {
182  mDestProjString = stripDatumTransform( mDestProjString );
183  }
184  if ( mDestinationDatumTransform != -1 )
185  {
186  mDestProjString += ( ' ' + datumTransformString( mDestinationDatumTransform ) );
187  }
188 
189  if ( !useDefaultDatumTransform )
190  {
191  addNullGridShifts( mSourceProjString, mDestProjString );
192  }
193 
194  // create proj projections for current thread
195  QPair< projPJ, projPJ > res = threadLocalProjData();
196 
197 #ifdef COORDINATE_TRANSFORM_VERBOSE
198  QgsDebugMsg( "From proj : " + mSourceCRS.toProj4() );
199  QgsDebugMsg( "To proj : " + mDestCRS.toProj4() );
200 #endif
201 
202  mInitialisedFlag = true;
203  if ( !res.second )
204  {
205  mInitialisedFlag = false;
206  }
207  if ( !res.first )
208  {
209  mInitialisedFlag = false;
210  }
211 #ifdef COORDINATE_TRANSFORM_VERBOSE
212  if ( mInitialisedFlag )
213  {
214  QgsDebugMsg( "------------------------------------------------------------" );
215  QgsDebugMsg( "The OGR Coordinate transformation for this layer was set to" );
216  QgsLogger::debug<QgsCoordinateReferenceSystem>( "Input", mSourceCRS, __FILE__, __FUNCTION__, __LINE__ );
217  QgsLogger::debug<QgsCoordinateReferenceSystem>( "Output", mDestCRS, __FILE__, __FUNCTION__, __LINE__ );
218  QgsDebugMsg( "------------------------------------------------------------" );
219  }
220  else
221  {
222  QgsDebugMsg( "------------------------------------------------------------" );
223  QgsDebugMsg( "The OGR Coordinate transformation FAILED TO INITIALISE!" );
224  QgsDebugMsg( "------------------------------------------------------------" );
225  }
226 #else
227  if ( !mInitialisedFlag )
228  {
229  QgsDebugMsg( "Coordinate transformation failed to initialize!" );
230  }
231 #endif
232 
233  //XXX todo overload == operator for QgsCoordinateReferenceSystem
234  //at the moment srs.parameters contains the whole proj def...soon it wont...
235  //if (mSourceCRS->toProj4() == mDestCRS->toProj4())
236  if ( mSourceCRS == mDestCRS )
237  {
238  // If the source and destination projection are the same, set the short
239  // circuit flag (no transform takes place)
240  mShortCircuit = true;
241  QgsDebugMsgLevel( "Source/Dest CRS equal, shortcircuit is set.", 3 );
242  }
243  else
244  {
245  // Transform must take place
246  mShortCircuit = false;
247  QgsDebugMsgLevel( "Source/Dest CRS UNequal, shortcircuit is NOt set.", 3 );
248  }
249 
250 }
251 
252 //
253 //
254 // TRANSFORMERS BELOW THIS POINT .........
255 //
256 //
257 //
258 
259 
261 {
262  if ( mShortCircuit || !mInitialisedFlag )
263  return thePoint;
264  // transform x
265  double x = thePoint.x();
266  double y = thePoint.y();
267  double z = 0.0;
268  try
269  {
270  transformCoords( 1, &x, &y, &z, direction );
271  }
272  catch ( const QgsCsException & )
273  {
274  // rethrow the exception
275  QgsDebugMsg( "rethrowing exception" );
276  throw;
277  }
278 
279  return QgsPoint( x, y );
280 }
281 
282 
283 QgsPoint QgsCoordinateTransform::transform( const double theX, const double theY = 0.0, TransformDirection direction ) const
284 {
285  try
286  {
287  return transform( QgsPoint( theX, theY ), direction );
288  }
289  catch ( const QgsCsException & )
290  {
291  // rethrow the exception
292  QgsDebugMsg( "rethrowing exception" );
293  throw;
294  }
295 }
296 
298 {
299  if ( mShortCircuit || !mInitialisedFlag )
300  return theRect;
301  // transform x
302  double x1 = theRect.xMinimum();
303  double y1 = theRect.yMinimum();
304  double x2 = theRect.xMaximum();
305  double y2 = theRect.yMaximum();
306 
307  // Number of points to reproject------+
308  // |
309  // V
310  try
311  {
312  double z = 0.0;
313  transformCoords( 1, &x1, &y1, &z, direction );
314  transformCoords( 1, &x2, &y2, &z, direction );
315  }
316  catch ( const QgsCsException & )
317  {
318  // rethrow the exception
319  QgsDebugMsg( "rethrowing exception" );
320  throw;
321  }
322 
323 #ifdef COORDINATE_TRANSFORM_VERBOSE
324  QgsDebugMsg( "Rect projection..." );
325  QgsDebugMsg( QString( "Xmin : %1 --> %2" ).arg( theRect.xMinimum() ).arg( x1 ) );
326  QgsDebugMsg( QString( "Ymin : %1 --> %2" ).arg( theRect.yMinimum() ).arg( y1 ) );
327  QgsDebugMsg( QString( "Xmax : %1 --> %2" ).arg( theRect.xMaximum() ).arg( x2 ) );
328  QgsDebugMsg( QString( "Ymax : %1 --> %2" ).arg( theRect.yMaximum() ).arg( y2 ) );
329 #endif
330  return QgsRectangle( x1, y1, x2, y2 );
331 }
332 
333 void QgsCoordinateTransform::transformInPlace( double& x, double& y, double& z,
334  TransformDirection direction ) const
335 {
336  if ( mShortCircuit || !mInitialisedFlag )
337  return;
338 #ifdef QGISDEBUG
339 // QgsDebugMsg(QString("Using transform in place %1 %2").arg(__FILE__).arg(__LINE__));
340 #endif
341  // transform x
342  try
343  {
344  transformCoords( 1, &x, &y, &z, direction );
345  }
346  catch ( const QgsCsException & )
347  {
348  // rethrow the exception
349  QgsDebugMsg( "rethrowing exception" );
350  throw;
351  }
352 }
353 
354 void QgsCoordinateTransform::transformInPlace( float& x, float& y, double& z,
355  TransformDirection direction ) const
356 {
357  double xd = static_cast< double >( x ), yd = static_cast< double >( y );
358  transformInPlace( xd, yd, z, direction );
359  x = xd;
360  y = yd;
361 }
362 
363 void QgsCoordinateTransform::transformInPlace( float& x, float& y, float& z,
364  TransformDirection direction ) const
365 {
366  if ( mShortCircuit || !mInitialisedFlag )
367  return;
368 #ifdef QGISDEBUG
369  // QgsDebugMsg(QString("Using transform in place %1 %2").arg(__FILE__).arg(__LINE__));
370 #endif
371  // transform x
372  try
373  {
374  double xd = x;
375  double yd = y;
376  double zd = z;
377  transformCoords( 1, &xd, &yd, &zd, direction );
378  x = xd;
379  y = yd;
380  z = zd;
381  }
382  catch ( QgsCsException & )
383  {
384  // rethrow the exception
385  QgsDebugMsg( "rethrowing exception" );
386  throw;
387  }
388 }
389 
391 {
392  if ( mShortCircuit || !mInitialisedFlag )
393  {
394  return;
395  }
396 
397  //create x, y arrays
398  int nVertices = poly.size();
399 
400  QVector<double> x( nVertices );
401  QVector<double> y( nVertices );
402  QVector<double> z( nVertices );
403 
404  for ( int i = 0; i < nVertices; ++i )
405  {
406  const QPointF& pt = poly.at( i );
407  x[i] = pt.x();
408  y[i] = pt.y();
409  z[i] = 0;
410  }
411 
412  try
413  {
414  transformCoords( nVertices, x.data(), y.data(), z.data(), direction );
415  }
416  catch ( const QgsCsException & )
417  {
418  // rethrow the exception
419  QgsDebugMsg( "rethrowing exception" );
420  throw;
421  }
422 
423  for ( int i = 0; i < nVertices; ++i )
424  {
425  QPointF& pt = poly[i];
426  pt.rx() = x[i];
427  pt.ry() = y[i];
428  }
429 }
430 
433  TransformDirection direction ) const
434 {
435  if ( mShortCircuit || !mInitialisedFlag )
436  return;
437 
438  Q_ASSERT( x.size() == y.size() );
439 
440  // Apparently, if one has a std::vector, it is valid to use the
441  // address of the first element in the vector as a pointer to an
442  // array of the vectors data, and hence easily interface with code
443  // that wants C-style arrays.
444 
445  try
446  {
447  transformCoords( x.size(), &x[0], &y[0], &z[0], direction );
448  }
449  catch ( const QgsCsException & )
450  {
451  // rethrow the exception
452  QgsDebugMsg( "rethrowing exception" );
453  throw;
454  }
455 }
456 
457 
460  TransformDirection direction ) const
461 {
462  if ( mShortCircuit || !mInitialisedFlag )
463  return;
464 
465  Q_ASSERT( x.size() == y.size() );
466 
467  // Apparently, if one has a std::vector, it is valid to use the
468  // address of the first element in the vector as a pointer to an
469  // array of the vectors data, and hence easily interface with code
470  // that wants C-style arrays.
471 
472  try
473  {
474  //copy everything to double vectors since proj needs double
475  int vectorSize = x.size();
476  QVector<double> xd( x.size() );
477  QVector<double> yd( y.size() );
478  QVector<double> zd( z.size() );
479  for ( int i = 0; i < vectorSize; ++i )
480  {
481  xd[i] = x[i];
482  yd[i] = y[i];
483  zd[i] = z[i];
484  }
485  transformCoords( x.size(), &xd[0], &yd[0], &zd[0], direction );
486 
487  //copy back
488  for ( int i = 0; i < vectorSize; ++i )
489  {
490  x[i] = xd[i];
491  y[i] = yd[i];
492  z[i] = zd[i];
493  }
494  }
495  catch ( QgsCsException & )
496  {
497  // rethrow the exception
498  QgsDebugMsg( "rethrowing exception" );
499  throw;
500  }
501 }
502 
503 QgsRectangle QgsCoordinateTransform::transformBoundingBox( const QgsRectangle &rect, TransformDirection direction, const bool handle180Crossover ) const
504 {
505  // Calculate the bounding box of a QgsRectangle in the source CRS
506  // when projected to the destination CRS (or the inverse).
507  // This is done by looking at a number of points spread evenly
508  // across the rectangle
509 
510  if ( mShortCircuit || !mInitialisedFlag )
511  return rect;
512 
513  if ( rect.isEmpty() )
514  {
515  QgsPoint p = transform( rect.xMinimum(), rect.yMinimum(), direction );
516  return QgsRectangle( p, p );
517  }
518 
519  // 64 points (<=2.12) is not enough, see #13665, for EPSG:4326 -> EPSG:3574 (say that it is a hard one),
520  // are decent result from about 500 points and more. This method is called quite often, but
521  // even with 1000 points it takes < 1ms
522  // TODO: how to effectively and precisely reproject bounding box?
523  const int nPoints = 1000;
524  double d = sqrt(( rect.width() * rect.height() ) / pow( sqrt( static_cast< double >( nPoints ) ) - 1, 2.0 ) );
525  int nXPoints = static_cast< int >( ceil( rect.width() / d ) ) + 1;
526  int nYPoints = static_cast< int >( ceil( rect.height() / d ) ) + 1;
527 
528  QgsRectangle bb_rect;
529  bb_rect.setMinimal();
530 
531  // We're interfacing with C-style vectors in the
532  // end, so let's do C-style vectors here too.
533 
534  QVector<double> x( nXPoints * nYPoints );
535  QVector<double> y( nXPoints * nYPoints );
536  QVector<double> z( nXPoints * nYPoints );
537 
538  QgsDebugMsg( "Entering transformBoundingBox..." );
539 
540  // Populate the vectors
541 
542  double dx = rect.width() / static_cast< double >( nXPoints - 1 );
543  double dy = rect.height() / static_cast< double >( nYPoints - 1 );
544 
545  double pointY = rect.yMinimum();
546 
547  for ( int i = 0; i < nYPoints ; i++ )
548  {
549 
550  // Start at right edge
551  double pointX = rect.xMinimum();
552 
553  for ( int j = 0; j < nXPoints; j++ )
554  {
555  x[( i*nXPoints ) + j] = pointX;
556  y[( i*nXPoints ) + j] = pointY;
557  // and the height...
558  z[( i*nXPoints ) + j] = 0.0;
559  // QgsDebugMsg(QString("BBox coord: (%1, %2)").arg(x[(i*numP) + j]).arg(y[(i*numP) + j]));
560  pointX += dx;
561  }
562  pointY += dy;
563  }
564 
565  // Do transformation. Any exception generated must
566  // be handled in above layers.
567  try
568  {
569  transformCoords( nXPoints * nYPoints, x.data(), y.data(), z.data(), direction );
570  }
571  catch ( const QgsCsException & )
572  {
573  // rethrow the exception
574  QgsDebugMsg( "rethrowing exception" );
575  throw;
576  }
577 
578  // Calculate the bounding box and use that for the extent
579 
580  for ( int i = 0; i < nXPoints * nYPoints; i++ )
581  {
582  if ( !qIsFinite( x[i] ) || !qIsFinite( y[i] ) )
583  {
584  continue;
585  }
586 
587  if ( handle180Crossover )
588  {
589  //if crossing the date line, temporarily add 360 degrees to -ve longitudes
590  bb_rect.combineExtentWith( x[i] >= 0.0 ? x[i] : x[i] + 360.0, y[i] );
591  }
592  else
593  {
594  bb_rect.combineExtentWith( x[i], y[i] );
595  }
596  }
597 
598  if ( handle180Crossover )
599  {
600  //subtract temporary addition of 360 degrees from longitudes
601  if ( bb_rect.xMinimum() > 180.0 )
602  bb_rect.setXMinimum( bb_rect.xMinimum() - 360.0 );
603  if ( bb_rect.xMaximum() > 180.0 )
604  bb_rect.setXMaximum( bb_rect.xMaximum() - 360.0 );
605  }
606 
607  QgsDebugMsg( "Projected extent: " + bb_rect.toString() );
608 
609  if ( bb_rect.isEmpty() )
610  {
611  QgsDebugMsg( "Original extent: " + rect.toString() );
612  }
613 
614  return bb_rect;
615 }
616 
617 void QgsCoordinateTransform::transformCoords( int numPoints, double *x, double *y, double *z, TransformDirection direction ) const
618 {
619  if ( mShortCircuit || !mInitialisedFlag )
620  return;
621  // Refuse to transform the points if the srs's are invalid
622  if ( !mSourceCRS.isValid() )
623  {
624  QgsMessageLog::logMessage( tr( "The source spatial reference system (CRS) is not valid. "
625  "The coordinates can not be reprojected. The CRS is: %1" )
626  .arg( mSourceCRS.toProj4() ), tr( "CRS" ) );
627  return;
628  }
629  if ( !mDestCRS.isValid() )
630  {
631  QgsMessageLog::logMessage( tr( "The destination spatial reference system (CRS) is not valid. "
632  "The coordinates can not be reprojected. The CRS is: %1" ).arg( mDestCRS.toProj4() ), tr( "CRS" ) );
633  return;
634  }
635 
636 #ifdef COORDINATE_TRANSFORM_VERBOSE
637  double xorg = *x;
638  double yorg = *y;
639  QgsDebugMsg( QString( "[[[[[[ Number of points to transform: %1 ]]]]]]" ).arg( numPoints ) );
640 #endif
641 
642  // use proj4 to do the transform
643  QString dir;
644  // if the source/destination projection is lat/long, convert the points to radians
645  // prior to transforming
646  QPair< projPJ, projPJ > projData = threadLocalProjData();
647  projPJ sourceProj = projData.first;
648  projPJ destProj = projData.second;
649 
650  if (( pj_is_latlong( destProj ) && ( direction == ReverseTransform ) )
651  || ( pj_is_latlong( sourceProj ) && ( direction == ForwardTransform ) ) )
652  {
653  for ( int i = 0; i < numPoints; ++i )
654  {
655  x[i] *= DEG_TO_RAD;
656  y[i] *= DEG_TO_RAD;
657  }
658 
659  }
660  int projResult;
661  if ( direction == ReverseTransform )
662  {
663  projResult = pj_transform( destProj, sourceProj, numPoints, 0, x, y, z );
664  }
665  else
666  {
667  Q_ASSERT( sourceProj );
668  Q_ASSERT( destProj );
669  projResult = pj_transform( sourceProj, destProj, numPoints, 0, x, y, z );
670  }
671 
672  if ( projResult != 0 )
673  {
674  //something bad happened....
675  QString points;
676 
677  for ( int i = 0; i < numPoints; ++i )
678  {
679  if ( direction == ForwardTransform )
680  {
681  points += QString( "(%1, %2)\n" ).arg( x[i], 0, 'f' ).arg( y[i], 0, 'f' );
682  }
683  else
684  {
685  points += QString( "(%1, %2)\n" ).arg( x[i] * RAD_TO_DEG, 0, 'f' ).arg( y[i] * RAD_TO_DEG, 0, 'f' );
686  }
687  }
688 
689  dir = ( direction == ForwardTransform ) ? tr( "forward transform" ) : tr( "inverse transform" );
690 
691  char *srcdef = pj_get_def( sourceProj, 0 );
692  char *dstdef = pj_get_def( destProj, 0 );
693 
694  QString msg = tr( "%1 of\n"
695  "%2"
696  "PROJ.4: %3 +to %4\n"
697  "Error: %5" )
698  .arg( dir,
699  points,
700  srcdef, dstdef,
701  QString::fromUtf8( pj_strerrno( projResult ) ) );
702 
703  pj_dalloc( srcdef );
704  pj_dalloc( dstdef );
705 
706  QgsDebugMsg( "Projection failed emitting invalid transform signal: " + msg );
707 
708  emit invalidTransformInput();
709 
710  QgsDebugMsg( "throwing exception" );
711 
712  throw QgsCsException( msg );
713  }
714 
715  // if the result is lat/long, convert the results from radians back
716  // to degrees
717  if (( pj_is_latlong( destProj ) && ( direction == ForwardTransform ) )
718  || ( pj_is_latlong( sourceProj ) && ( direction == ReverseTransform ) ) )
719  {
720  for ( int i = 0; i < numPoints; ++i )
721  {
722  x[i] *= RAD_TO_DEG;
723  y[i] *= RAD_TO_DEG;
724  }
725  }
726 #ifdef COORDINATE_TRANSFORM_VERBOSE
727  QgsDebugMsg( QString( "[[[[[[ Projected %1, %2 to %3, %4 ]]]]]]" )
728  .arg( xorg, 0, 'g', 15 ).arg( yorg, 0, 'g', 15 )
729  .arg( *x, 0, 'g', 15 ).arg( *y, 0, 'g', 15 ) );
730 #endif
731 }
732 
734 {
735 
736  QgsDebugMsg( "Reading Coordinate Transform from xml ------------------------!" );
737 
738  QDomNode mySrcNode = theNode.namedItem( "sourcesrs" );
739  mSourceCRS.readXML( mySrcNode );
740 
741  QDomNode myDestNode = theNode.namedItem( "destinationsrs" );
742  mDestCRS.readXML( myDestNode );
743 
744  mSourceDatumTransform = theNode.toElement().attribute( "sourceDatumTransform", "-1" ).toInt();
745  mDestinationDatumTransform = theNode.toElement().attribute( "destinationDatumTransform", "-1" ).toInt();
746 
747  initialise();
748 
749  return true;
750 }
751 
753 {
754  QDomElement myNodeElement = theNode.toElement();
755  QDomElement myTransformElement = theDoc.createElement( "coordinatetransform" );
756  myTransformElement.setAttribute( "sourceDatumTransform", QString::number( mSourceDatumTransform ) );
757  myTransformElement.setAttribute( "destinationDatumTransform", QString::number( mDestinationDatumTransform ) );
758 
759  QDomElement mySourceElement = theDoc.createElement( "sourcesrs" );
760  mSourceCRS.writeXML( mySourceElement, theDoc );
761  myTransformElement.appendChild( mySourceElement );
762 
763  QDomElement myDestElement = theDoc.createElement( "destinationsrs" );
764  mDestCRS.writeXML( myDestElement, theDoc );
765  myTransformElement.appendChild( myDestElement );
766 
767  myNodeElement.appendChild( myTransformElement );
768 
769  return true;
770 }
771 
772 const char *finder( const char *name )
773 {
774  QString proj;
775 #ifdef Q_OS_WIN
777  + "/share/proj/" + QString( name );
778 #else
779  Q_UNUSED( name );
780 #endif
781  return proj.toUtf8();
782 }
783 
784 void QgsCoordinateTransform::setFinder()
785 {
786 #if 0
787  // Attention! It should be possible to set PROJ_LIB
788  // but it can happen that it was previously set by installer
789  // (version 0.7) and the old installation was deleted
790 
791  // Another problem: PROJ checks if pj_finder was set before
792  // PROJ_LIB environment variable. pj_finder is probably set in
793  // GRASS gproj library when plugin is loaded, consequently
794  // PROJ_LIB is ignored
795 
796  pj_set_finder( finder );
797 #endif
798 }
799 
801 {
802  QList< QList< int > > transformations;
803 
804  QString srcGeoId = srcCRS.geographicCRSAuthId();
805  QString destGeoId = destCRS.geographicCRSAuthId();
806 
807  if ( srcGeoId.isEmpty() || destGeoId.isEmpty() )
808  {
809  return transformations;
810  }
811 
812  QStringList srcSplit = srcGeoId.split( ':' );
813  QStringList destSplit = destGeoId.split( ':' );
814 
815  if ( srcSplit.size() < 2 || destSplit.size() < 2 )
816  {
817  return transformations;
818  }
819 
820  int srcAuthCode = srcSplit.at( 1 ).toInt();
821  int destAuthCode = destSplit.at( 1 ).toInt();
822 
823  if ( srcAuthCode == destAuthCode )
824  {
825  return transformations; //crs have the same datum
826  }
827 
828  QList<int> directTransforms;
829  searchDatumTransform( QString( "SELECT coord_op_code FROM tbl_datum_transform WHERE source_crs_code=%1 AND target_crs_code=%2 ORDER BY deprecated ASC,preferred DESC" ).arg( srcAuthCode ).arg( destAuthCode ),
830  directTransforms );
831  QList<int> reverseDirectTransforms;
832  searchDatumTransform( QString( "SELECT coord_op_code FROM tbl_datum_transform WHERE source_crs_code = %1 AND target_crs_code=%2 ORDER BY deprecated ASC,preferred DESC" ).arg( destAuthCode ).arg( srcAuthCode ),
833  reverseDirectTransforms );
834  QList<int> srcToWgs84;
835  searchDatumTransform( QString( "SELECT coord_op_code FROM tbl_datum_transform WHERE (source_crs_code=%1 AND target_crs_code=%2) OR (source_crs_code=%2 AND target_crs_code=%1) ORDER BY deprecated ASC,preferred DESC" ).arg( srcAuthCode ).arg( 4326 ),
836  srcToWgs84 );
837  QList<int> destToWgs84;
838  searchDatumTransform( QString( "SELECT coord_op_code FROM tbl_datum_transform WHERE (source_crs_code=%1 AND target_crs_code=%2) OR (source_crs_code=%2 AND target_crs_code=%1) ORDER BY deprecated ASC,preferred DESC" ).arg( destAuthCode ).arg( 4326 ),
839  destToWgs84 );
840 
841  //add direct datum transformations
842  QList<int>::const_iterator directIt = directTransforms.constBegin();
843  for ( ; directIt != directTransforms.constEnd(); ++directIt )
844  {
845  transformations.push_back( QList<int>() << *directIt << -1 );
846  }
847 
848  //add direct datum transformations
849  directIt = reverseDirectTransforms.constBegin();
850  for ( ; directIt != reverseDirectTransforms.constEnd(); ++directIt )
851  {
852  transformations.push_back( QList<int>() << -1 << *directIt );
853  }
854 
855  QList<int>::const_iterator srcWgsIt = srcToWgs84.constBegin();
856  for ( ; srcWgsIt != srcToWgs84.constEnd(); ++srcWgsIt )
857  {
858  QList<int>::const_iterator dstWgsIt = destToWgs84.constBegin();
859  for ( ; dstWgsIt != destToWgs84.constEnd(); ++dstWgsIt )
860  {
861  transformations.push_back( QList<int>() << *srcWgsIt << *dstWgsIt );
862  }
863  }
864 
865  return transformations;
866 }
867 
868 QString QgsCoordinateTransform::stripDatumTransform( const QString& proj4 )
869 {
870  QStringList parameterSplit = proj4.split( '+', QString::SkipEmptyParts );
871  QString currentParameter;
872  QString newProjString;
873 
874  for ( int i = 0; i < parameterSplit.size(); ++i )
875  {
876  currentParameter = parameterSplit.at( i );
877  if ( !currentParameter.startsWith( "towgs84", Qt::CaseInsensitive )
878  && !currentParameter.startsWith( "nadgrids", Qt::CaseInsensitive ) )
879  {
880  newProjString.append( '+' );
881  newProjString.append( currentParameter );
882  newProjString.append( ' ' );
883  }
884  }
885  return newProjString;
886 }
887 
888 void QgsCoordinateTransform::searchDatumTransform( const QString& sql, QList< int >& transforms )
889 {
890  sqlite3* db;
891  int openResult = sqlite3_open_v2( QgsApplication::srsDbFilePath().toUtf8().constData(), &db, SQLITE_OPEN_READONLY, 0 );
892  if ( openResult != SQLITE_OK )
893  {
894  sqlite3_close( db );
895  return;
896  }
897 
898  sqlite3_stmt* stmt;
899  int prepareRes = sqlite3_prepare( db, sql.toAscii(), sql.size(), &stmt, nullptr );
900  if ( prepareRes != SQLITE_OK )
901  {
902  sqlite3_finalize( stmt );
903  sqlite3_close( db );
904  return;
905  }
906 
907  QString cOpCode;
908  while ( sqlite3_step( stmt ) == SQLITE_ROW )
909  {
910  cOpCode = reinterpret_cast< const char * >( sqlite3_column_text( stmt, 0 ) );
911  transforms.push_back( cOpCode.toInt() );
912  }
913  sqlite3_finalize( stmt );
914  sqlite3_close( db );
915 }
916 
918 {
919  QString transformString;
920 
921  sqlite3* db;
922  int openResult = sqlite3_open_v2( QgsApplication::srsDbFilePath().toUtf8().constData(), &db, SQLITE_OPEN_READONLY, 0 );
923  if ( openResult != SQLITE_OK )
924  {
925  sqlite3_close( db );
926  return transformString;
927  }
928 
929  sqlite3_stmt* stmt;
930  QString sql = QString( "SELECT coord_op_method_code,p1,p2,p3,p4,p5,p6,p7 FROM tbl_datum_transform WHERE coord_op_code=%1" ).arg( datumTransform );
931  int prepareRes = sqlite3_prepare( db, sql.toAscii(), sql.size(), &stmt, nullptr );
932  if ( prepareRes != SQLITE_OK )
933  {
934  sqlite3_finalize( stmt );
935  sqlite3_close( db );
936  return transformString;
937  }
938 
939  if ( sqlite3_step( stmt ) == SQLITE_ROW )
940  {
941  //coord_op_methode_code
942  int methodCode = sqlite3_column_int( stmt, 0 );
943  if ( methodCode == 9615 ) //ntv2
944  {
945  transformString = "+nadgrids=" + QString( reinterpret_cast< const char * >( sqlite3_column_text( stmt, 1 ) ) );
946  }
947  else if ( methodCode == 9603 || methodCode == 9606 || methodCode == 9607 )
948  {
949  transformString += "+towgs84=";
950  double p1 = sqlite3_column_double( stmt, 1 );
951  double p2 = sqlite3_column_double( stmt, 2 );
952  double p3 = sqlite3_column_double( stmt, 3 );
953  double p4 = sqlite3_column_double( stmt, 4 );
954  double p5 = sqlite3_column_double( stmt, 5 );
955  double p6 = sqlite3_column_double( stmt, 6 );
956  double p7 = sqlite3_column_double( stmt, 7 );
957  if ( methodCode == 9603 ) //3 parameter transformation
958  {
959  transformString += QString( "%1,%2,%3" ).arg( p1 ).arg( p2 ).arg( p3 );
960  }
961  else //7 parameter transformation
962  {
963  transformString += QString( "%1,%2,%3,%4,%5,%6,%7" ).arg( p1 ).arg( p2 ).arg( p3 ).arg( p4 ).arg( p5 ).arg( p6 ).arg( p7 );
964  }
965  }
966  }
967 
968  sqlite3_finalize( stmt );
969  sqlite3_close( db );
970  return transformString;
971 }
972 
973 bool QgsCoordinateTransform::datumTransformCrsInfo( int datumTransform, int& epsgNr, QString& srcProjection, QString& dstProjection, QString &remarks, QString &scope, bool &preferred, bool &deprecated )
974 {
975  sqlite3* db;
976  int openResult = sqlite3_open_v2( QgsApplication::srsDbFilePath().toUtf8().constData(), &db, SQLITE_OPEN_READONLY, 0 );
977  if ( openResult != SQLITE_OK )
978  {
979  sqlite3_close( db );
980  return false;
981  }
982 
983  sqlite3_stmt* stmt;
984  QString sql = QString( "SELECT epsg_nr,source_crs_code,target_crs_code,remarks,scope,preferred,deprecated FROM tbl_datum_transform WHERE coord_op_code=%1" ).arg( datumTransform );
985  int prepareRes = sqlite3_prepare( db, sql.toAscii(), sql.size(), &stmt, nullptr );
986  if ( prepareRes != SQLITE_OK )
987  {
988  sqlite3_finalize( stmt );
989  sqlite3_close( db );
990  return false;
991  }
992 
993  int srcCrsId, destCrsId;
994  if ( sqlite3_step( stmt ) != SQLITE_ROW )
995  {
996  sqlite3_finalize( stmt );
997  sqlite3_close( db );
998  return false;
999  }
1000 
1001  epsgNr = sqlite3_column_int( stmt, 0 );
1002  srcCrsId = sqlite3_column_int( stmt, 1 );
1003  destCrsId = sqlite3_column_int( stmt, 2 );
1004  remarks = QString::fromUtf8( reinterpret_cast< const char * >( sqlite3_column_text( stmt, 3 ) ) );
1005  scope = QString::fromUtf8( reinterpret_cast< const char * >( sqlite3_column_text( stmt, 4 ) ) );
1006  preferred = sqlite3_column_int( stmt, 5 ) != 0;
1007  deprecated = sqlite3_column_int( stmt, 6 ) != 0;
1008 
1009  QgsCoordinateReferenceSystem srcCrs = QgsCRSCache::instance()->crsByOgcWmsCrs( QString( "EPSG:%1" ).arg( srcCrsId ) );
1010  srcProjection = srcCrs.description();
1011  QgsCoordinateReferenceSystem destCrs = QgsCRSCache::instance()->crsByOgcWmsCrs( QString( "EPSG:%1" ).arg( destCrsId ) );
1012  dstProjection = destCrs.description();
1013 
1014  sqlite3_finalize( stmt );
1015  sqlite3_close( db );
1016  return true;
1017 }
1018 
1019 void QgsCoordinateTransform::addNullGridShifts( QString& srcProjString, QString& destProjString )
1020 {
1021  //if one transformation uses ntv2, the other one needs to be null grid shift
1022  if ( mDestinationDatumTransform == -1 && srcProjString.contains( "+nadgrids" ) ) //add null grid if source transformation is ntv2
1023  {
1024  destProjString += " [email protected]";
1025  return;
1026  }
1027  if ( mSourceDatumTransform == -1 && destProjString.contains( "+nadgrids" ) )
1028  {
1029  srcProjString += " [email protected]";
1030  return;
1031  }
1032 
1033  //add null shift grid for google mercator
1034  //(see e.g. http://trac.osgeo.org/proj/wiki/FAQ#ChangingEllipsoidWhycantIconvertfromWGS84toGoogleEarthVirtualGlobeMercator)
1035  if ( mSourceCRS.authid().compare( "EPSG:3857", Qt::CaseInsensitive ) == 0 && mSourceDatumTransform == -1 )
1036  {
1037  srcProjString += " [email protected]";
1038  }
1039  if ( mDestCRS.authid().compare( "EPSG:3857", Qt::CaseInsensitive ) == 0 && mDestinationDatumTransform == -1 )
1040  {
1041  destProjString += " [email protected]";
1042  }
1043 }
1044 
1045 QPair<projPJ, projPJ> QgsCoordinateTransform::threadLocalProjData() const
1046 {
1047  mProjLock.lockForRead();
1048  projCtx pContext = nullptr;
1049  if ( mProjContext.hasLocalData() )
1050  pContext = mProjContext.localData()->get();
1051  else
1052  {
1053  mProjContext.setLocalData( new QgsProjContextStore() );
1054  pContext = mProjContext.localData()->get();
1055  }
1056 
1057  QMap< uintptr_t, QPair< projPJ, projPJ > >::const_iterator it = mProjProjections.constFind( reinterpret_cast< uintptr_t >( pContext ) );
1058  if ( it != mProjProjections.constEnd() )
1059  {
1060  QPair< projPJ, projPJ > res = it.value();
1061  mProjLock.unlock();
1062  return res;
1063  }
1064 
1065  // proj projections don't exist yet, so we need to create
1066  mProjLock.unlock();
1067  mProjLock.lockForWrite();
1068  QPair< projPJ, projPJ > res = qMakePair( pj_init_plus_ctx( pContext, mSourceProjString.toUtf8() ),
1069  pj_init_plus_ctx( pContext, mDestProjString.toUtf8() ) );
1070  mProjProjections.insert( reinterpret_cast< uintptr_t >( pContext ), res );
1071  mProjLock.unlock();
1072  return res;
1073 }
1074 
1075 void QgsCoordinateTransform::freeProj()
1076 {
1077  mProjLock.lockForWrite();
1078  QMap< uintptr_t, QPair< projPJ, projPJ > >::const_iterator it = mProjProjections.constBegin();
1079  for ( ; it != mProjProjections.constEnd(); ++it )
1080  {
1081  pj_free( it.value().first );
1082  pj_free( it.value().second );
1083  }
1084  mProjProjections.clear();
1085  mProjLock.unlock();
1086 }
1087 
1088 QgsCoordinateTransform::QgsProjContextStore::QgsProjContextStore()
1089 {
1090  context = pj_ctx_alloc();
1091 }
1092 
1093 QgsCoordinateTransform::QgsProjContextStore::~QgsProjContextStore()
1094 {
1095  pj_ctx_free( context );
1096 }
QgsPoint transform(const QgsPoint &p, TransformDirection direction=ForwardTransform) const
Transform the point from Source Coordinate System to Destination Coordinate System If the direction i...
void lockForWrite()
void transformCoords(int numPoint, double *x, double *y, double *z, TransformDirection direction=ForwardTransform) const
Transform an array of coordinates to a different Coordinate System If the direction is ForwardTransfo...
A rectangle specified with double values.
Definition: qgsrectangle.h:35
QString & append(QChar ch)
void lockForRead()
QgsCoordinateReferenceSystem crsByOgcWmsCrs(const QString &ogcCrs) const
Returns the CRS from a given OGC WMS-format Coordinate Reference System string.
void setMinimal()
Set a rectangle so that min corner is at max and max corner is at min.
QDomNode appendChild(const QDomNode &newChild)
void setXMaximum(double x)
Set the maximum x value.
Definition: qgsrectangle.h:172
void push_back(const T &value)
QString attribute(const QString &name, const QString &defValue) const
static QList< QList< int > > datumTransformations(const QgsCoordinateReferenceSystem &srcCRS, const QgsCoordinateReferenceSystem &destCRS)
Returns list of datum transformations for the given src and dest CRS.
#define QgsDebugMsg(str)
Definition: qgslogger.h:33
void setSourceCrs(const QgsCoordinateReferenceSystem &theCRS)
QStringList split(const QString &sep, SplitBehavior behavior, Qt::CaseSensitivity cs) const
const_iterator constBegin() const
TransformDirection
Enum used to indicate the direction (forward or inverse) of the transform.
const T & at(int i) const
int size() const
bool readXML(QDomNode &theNode)
Restores state from the given Dom node.
QString toProj4() const
Returns a Proj4 string representation of this CRS.
bool createFromId(const long theId, CrsType theType=PostgisCrsId)
void initialise()
initialize is used to actually create the Transformer instance
QgsCoordinateReferenceSystem crsByWkt(const QString &wkt) const
Returns the CRS from a WKT spatial ref sys definition string.
const_iterator constFind(const Key &key) const
void clear()
QString tr(const char *sourceText, const char *disambiguation, int n)
QgsCoordinateReferenceSystem crsBySrsId(long srsId) const
Returns the CRS from a specified QGIS SRS ID.
int size() const
double y() const
Get the y value of the point.
Definition: qgspoint.h:193
void setLocalData(T data)
QDomElement toElement() const
QgsRectangle transformBoundingBox(const QgsRectangle &theRect, TransformDirection direction=ForwardTransform, const bool handle180Crossover=false) const
Transform a QgsRectangle to the dest Coordinate system If the direction is ForwardTransform then coor...
QString toString(bool automaticPrecision=false) const
returns string representation of form xmin,ymin xmax,ymax
T * data()
const char * name() const
QgsCoordinateTransform()
Default constructor.
QString number(int n, int base)
qreal x() const
qreal y() const
bool hasLocalData() const
QString fromUtf8(const char *str, int size)
const QgsCoordinateReferenceSystem & sourceCrs() const
bool writeXML(QDomNode &theNode, QDomDocument &theDoc)
Stores state to the given Dom node in the given document.
QString geographicCRSAuthId() const
Returns auth id of related geographic CRS.
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:34
bool isEmpty() const
test if rectangle is empty.
bool readXML(const QDomNode &theNode)
Restores state from the given Dom node.
void setAttribute(const QString &name, const QString &value)
static QString datumTransformString(int datumTransform)
double width() const
Width of the rectangle.
Definition: qgsrectangle.h:207
int toInt(bool *ok, int base) const
bool isEmpty() const
const_iterator constEnd() const
bool startsWith(const QString &s, Qt::CaseSensitivity cs) const
static void logMessage(const QString &message, const QString &tag=QString::null, MessageLevel level=WARNING)
add a message to the instance (and create it if necessary)
QString description() const
Returns the descriptive name of the CRS, eg "WGS 84" or "GDA 94 / Vicgrid94".
void * projCtx
A class to represent a point.
Definition: qgspoint.h:117
Caches QgsCoordinateReferenceSystem construction, which may be expensive.
Definition: qgscrscache.h:67
QDomNode namedItem(const QString &name) const
struct sqlite3 sqlite3
bool contains(QChar ch, Qt::CaseSensitivity cs) const
double yMinimum() const
Get the y minimum value (bottom side of rectangle)
Definition: qgsrectangle.h:202
double xMaximum() const
Get the x maximum value (right side of rectangle)
Definition: qgsrectangle.h:187
void combineExtentWith(const QgsRectangle &rect)
expand the rectangle so that covers both the original rectangle and the given rectangle ...
void setDestCRSID(long theCRSID)
Change the destination coordinate system by passing it a qgis srsid A QGIS srsid is a unique key valu...
void setDestCRS(const QgsCoordinateReferenceSystem &theCRS)
void * projPJ
const T & at(int i) const
Class for storing a coordinate reference system (CRS)
QgsCoordinateTransform * clone() const
const QgsCoordinateReferenceSystem & destCRS() const
qreal & rx()
qreal & ry()
Class for doing transforms between two map coordinate systems.
double xMinimum() const
Get the x minimum value (left side of rectangle)
Definition: qgsrectangle.h:192
double yMaximum() const
Get the y maximum value (top side of rectangle)
Definition: qgsrectangle.h:197
iterator insert(const Key &key, const T &value)
static QString srsDbFilePath()
Returns the path to the srs.db file.
void invalidTransformInput() const
Signal when an invalid pj_transform() has occurred.
Custom exception class for Coordinate Reference System related exceptions.
void transformPolygon(QPolygonF &poly, TransformDirection direction=ForwardTransform) const
const_iterator constEnd() const
bool writeXML(QDomNode &theNode, QDomDocument &theDoc) const
Stores state to the given Dom node in the given document.
QDomElement createElement(const QString &tagName)
const_iterator constBegin() const
static bool datumTransformCrsInfo(int datumTransform, int &epsgNr, QString &srcProjection, QString &dstProjection, QString &remarks, QString &scope, bool &preferred, bool &deprecated)
Gets name of source and dest geographical CRS (to show in a tooltip)
QString applicationDirPath()
void transformInPlace(double &x, double &y, double &z, TransformDirection direction=ForwardTransform) const
int size() const
int compare(const QString &other) const
QString arg(qlonglong a, int fieldWidth, int base, const QChar &fillChar) const
static QgsCRSCache * instance()
Returns a pointer to the QgsCRSCache singleton.
Definition: qgscrscache.cpp:91
QString authid() const
Returns the authority identifier for the CRS, which includes both the authority (eg EPSG) and the CRS...
double x() const
Get the x value of the point.
Definition: qgspoint.h:185
void setXMinimum(double x)
Set the minimum x value.
Definition: qgsrectangle.h:167
QByteArray toAscii() const
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:212
const char * finder(const char *name)
const T value(const Key &key) const
bool isValid() const
Returns whether this CRS is correctly initialized and usable.
QByteArray toUtf8() const