QGIS API Documentation  2.17.0-Master (872e6d2)
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 }
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...
const QgsCoordinateReferenceSystem & sourceCrs() const
A rectangle specified with double values.
Definition: qgsrectangle.h:35
QString & append(QChar ch)
bool isEmpty() const
test if rectangle is empty.
void setMinimal()
Set a rectangle so that min corner is at max and max corner is at min.
void invalidTransformInput() const
Signal when an invalid pj_transform() has occurred.
QgsCoordinateReferenceSystem crsBySrsId(long srsId) const
Returns the CRS from a specified QGIS SRS ID.
QDomNode appendChild(const QDomNode &newChild)
void setXMaximum(double x)
Set the maximum x value.
Definition: qgsrectangle.h:172
void push_back(const T &value)
QgsCoordinateReferenceSystem crsByWkt(const QString &wkt) const
Returns the CRS from a WKT spatial ref sys definition string.
void transformPolygon(QPolygonF &poly, TransformDirection direction=ForwardTransform) const
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.
double yMaximum() const
Get the y maximum value (top side of rectangle)
Definition: qgsrectangle.h:197
#define QgsDebugMsg(str)
Definition: qgslogger.h:33
QString geographicCRSAuthId() const
Returns auth id of related geographic CRS.
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.
bool createFromId(const long theId, CrsType theType=PostgisCrsId)
void initialise()
initialize is used to actually create the Transformer instance
QgsPoint transform(const QgsPoint &p, TransformDirection direction=ForwardTransform) const
Transform the point from Source Coordinate System to Destination Coordinate System If the direction i...
QString tr(const char *sourceText, const char *disambiguation, int n)
double x() const
Get the x value of the point.
Definition: qgspoint.h:185
int size() const
QDomElement toElement() const
T * data()
const char * name() const
void transformInPlace(double &x, double &y, double &z, TransformDirection direction=ForwardTransform) const
QgsCoordinateTransform()
Default constructor.
QString number(int n, int base)
qreal x() const
qreal y() const
QString fromUtf8(const char *str, int size)
bool writeXML(QDomNode &theNode, QDomDocument &theDoc)
Stores state to the given Dom node in the given document.
QgsCoordinateReferenceSystem crsByOgcWmsCrs(const QString &ogcCrs) const
Returns the CRS from a given OGC WMS-format Coordinate Reference System string.
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
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:34
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)
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)
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
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)
bool isValid() const
Returns whether this CRS is correctly initialized and usable.
const T & at(int i) const
bool writeXML(QDomNode &theNode, QDomDocument &theDoc) const
Stores state to the given Dom node in the given document.
Class for storing a coordinate reference system (CRS)
qreal & rx()
qreal & ry()
QString authid() const
Returns the authority identifier for the CRS, which includes both the authority (eg EPSG) and the CRS...
Class for doing transforms between two map coordinate systems.
double y() const
Get the y value of the point.
Definition: qgspoint.h:193
static QString srsDbFilePath()
Returns the path to the srs.db file.
QString description() const
Returns the descriptive name of the CRS, eg "WGS 84" or "GDA 94 / Vicgrid94".
Custom exception class for Coordinate Reference System related exceptions.
const_iterator constEnd() const
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)
double width() const
Width of the rectangle.
Definition: qgsrectangle.h:207
QString applicationDirPath()
int size() const
const QgsCoordinateReferenceSystem & destCRS() const
int compare(const QString &other) const
QString toString(bool automaticPrecision=false) const
returns string representation of form xmin,ymin xmax,ymax
QString arg(qlonglong a, int fieldWidth, int base, const QChar &fillChar) const
double xMinimum() const
Get the x minimum value (left side of rectangle)
Definition: qgsrectangle.h:192
static QgsCRSCache * instance()
Returns a pointer to the QgsCRSCache singleton.
Definition: qgscrscache.cpp:91
void setXMinimum(double x)
Set the minimum x value.
Definition: qgsrectangle.h:167
QByteArray toAscii() 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...
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:212
const char * finder(const char *name)
QgsCoordinateTransform * clone() const
QString toProj4() const
Returns a Proj4 string representation of this CRS.
QByteArray toUtf8() const