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