QGIS API Documentation  2.15.0-Master (af20121)
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( theSourceSrsId, QgsCoordinateReferenceSystem::InternalCrsId )
71  , mDestCRS( theDestSrsId, QgsCoordinateReferenceSystem::InternalCrsId )
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.createFromWkt( theSourceCRS );
90  mDestCRS.createFromWkt( 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.createFromWkt( 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.createFromSrsId( 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()->crsByAuthId( 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  z[i] *= DEG_TO_RAD;
672  }
673 
674  }
675  int projResult;
676  if ( direction == ReverseTransform )
677  {
678  projResult = pj_transform( mDestinationProjection, mSourceProjection, numPoints, 0, x, y, z );
679  }
680  else
681  {
682  Q_ASSERT( mSourceProjection );
683  Q_ASSERT( mDestinationProjection );
684  projResult = pj_transform( mSourceProjection, mDestinationProjection, numPoints, 0, x, y, z );
685  }
686 
687  if ( projResult != 0 )
688  {
689  //something bad happened....
690  QString points;
691 
692  for ( int i = 0; i < numPoints; ++i )
693  {
694  if ( direction == ForwardTransform )
695  {
696  points += QString( "(%1, %2)\n" ).arg( x[i], 0, 'f' ).arg( y[i], 0, 'f' );
697  }
698  else
699  {
700  points += QString( "(%1, %2)\n" ).arg( x[i] * RAD_TO_DEG, 0, 'f' ).arg( y[i] * RAD_TO_DEG, 0, 'f' );
701  }
702  }
703 
704  dir = ( direction == ForwardTransform ) ? tr( "forward transform" ) : tr( "inverse transform" );
705 
706  char *srcdef = pj_get_def( mSourceProjection, 0 );
707  char *dstdef = pj_get_def( mDestinationProjection, 0 );
708 
709  QString msg = tr( "%1 of\n"
710  "%2"
711  "PROJ.4: %3 +to %4\n"
712  "Error: %5" )
713  .arg( dir,
714  points,
715  srcdef, dstdef,
716  QString::fromUtf8( pj_strerrno( projResult ) ) );
717 
718  pj_dalloc( srcdef );
719  pj_dalloc( dstdef );
720 
721  QgsDebugMsg( "Projection failed emitting invalid transform signal: " + msg );
722 
723  emit invalidTransformInput();
724 
725  QgsDebugMsg( "throwing exception" );
726 
727  throw QgsCsException( msg );
728  }
729 
730  // if the result is lat/long, convert the results from radians back
731  // to degrees
732  if (( pj_is_latlong( mDestinationProjection ) && ( direction == ForwardTransform ) )
733  || ( pj_is_latlong( mSourceProjection ) && ( direction == ReverseTransform ) ) )
734  {
735  for ( int i = 0; i < numPoints; ++i )
736  {
737  x[i] *= RAD_TO_DEG;
738  y[i] *= RAD_TO_DEG;
739  z[i] *= RAD_TO_DEG;
740  }
741  }
742 #ifdef COORDINATE_TRANSFORM_VERBOSE
743  QgsDebugMsg( QString( "[[[[[[ Projected %1, %2 to %3, %4 ]]]]]]" )
744  .arg( xorg, 0, 'g', 15 ).arg( yorg, 0, 'g', 15 )
745  .arg( *x, 0, 'g', 15 ).arg( *y, 0, 'g', 15 ) );
746 #endif
747 }
748 
750 {
751 
752  QgsDebugMsg( "Reading Coordinate Transform from xml ------------------------!" );
753 
754  QDomNode mySrcNode = theNode.namedItem( "sourcesrs" );
755  mSourceCRS.readXML( mySrcNode );
756 
757  QDomNode myDestNode = theNode.namedItem( "destinationsrs" );
758  mDestCRS.readXML( myDestNode );
759 
760  mSourceDatumTransform = theNode.toElement().attribute( "sourceDatumTransform", "-1" ).toInt();
761  mDestinationDatumTransform = theNode.toElement().attribute( "destinationDatumTransform", "-1" ).toInt();
762 
763  initialise();
764 
765  return true;
766 }
767 
769 {
770  QDomElement myNodeElement = theNode.toElement();
771  QDomElement myTransformElement = theDoc.createElement( "coordinatetransform" );
772  myTransformElement.setAttribute( "sourceDatumTransform", QString::number( mSourceDatumTransform ) );
773  myTransformElement.setAttribute( "destinationDatumTransform", QString::number( mDestinationDatumTransform ) );
774 
775  QDomElement mySourceElement = theDoc.createElement( "sourcesrs" );
776  mSourceCRS.writeXML( mySourceElement, theDoc );
777  myTransformElement.appendChild( mySourceElement );
778 
779  QDomElement myDestElement = theDoc.createElement( "destinationsrs" );
780  mDestCRS.writeXML( myDestElement, theDoc );
781  myTransformElement.appendChild( myDestElement );
782 
783  myNodeElement.appendChild( myTransformElement );
784 
785  return true;
786 }
787 
788 const char *finder( const char *name )
789 {
790  QString proj;
791 #ifdef Q_OS_WIN
793  + "/share/proj/" + QString( name );
794 #else
795  Q_UNUSED( name );
796 #endif
797  return proj.toUtf8();
798 }
799 
800 void QgsCoordinateTransform::setFinder()
801 {
802 #if 0
803  // Attention! It should be possible to set PROJ_LIB
804  // but it can happen that it was previously set by installer
805  // (version 0.7) and the old installation was deleted
806 
807  // Another problem: PROJ checks if pj_finder was set before
808  // PROJ_LIB environment variable. pj_finder is probably set in
809  // GRASS gproj library when plugin is loaded, consequently
810  // PROJ_LIB is ignored
811 
812  pj_set_finder( finder );
813 #endif
814 }
815 
817 {
818  QList< QList< int > > transformations;
819 
820  QString srcGeoId = srcCRS.geographicCRSAuthId();
821  QString destGeoId = destCRS.geographicCRSAuthId();
822 
823  if ( srcGeoId.isEmpty() || destGeoId.isEmpty() )
824  {
825  return transformations;
826  }
827 
828  QStringList srcSplit = srcGeoId.split( ':' );
829  QStringList destSplit = destGeoId.split( ':' );
830 
831  if ( srcSplit.size() < 2 || destSplit.size() < 2 )
832  {
833  return transformations;
834  }
835 
836  int srcAuthCode = srcSplit.at( 1 ).toInt();
837  int destAuthCode = destSplit.at( 1 ).toInt();
838 
839  if ( srcAuthCode == destAuthCode )
840  {
841  return transformations; //crs have the same datum
842  }
843 
844  QList<int> directTransforms;
845  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 ),
846  directTransforms );
847  QList<int> reverseDirectTransforms;
848  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 ),
849  reverseDirectTransforms );
850  QList<int> srcToWgs84;
851  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 ),
852  srcToWgs84 );
853  QList<int> destToWgs84;
854  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 ),
855  destToWgs84 );
856 
857  //add direct datum transformations
858  QList<int>::const_iterator directIt = directTransforms.constBegin();
859  for ( ; directIt != directTransforms.constEnd(); ++directIt )
860  {
861  transformations.push_back( QList<int>() << *directIt << -1 );
862  }
863 
864  //add direct datum transformations
865  directIt = reverseDirectTransforms.constBegin();
866  for ( ; directIt != reverseDirectTransforms.constEnd(); ++directIt )
867  {
868  transformations.push_back( QList<int>() << -1 << *directIt );
869  }
870 
871  QList<int>::const_iterator srcWgsIt = srcToWgs84.constBegin();
872  for ( ; srcWgsIt != srcToWgs84.constEnd(); ++srcWgsIt )
873  {
874  QList<int>::const_iterator dstWgsIt = destToWgs84.constBegin();
875  for ( ; dstWgsIt != destToWgs84.constEnd(); ++dstWgsIt )
876  {
877  transformations.push_back( QList<int>() << *srcWgsIt << *dstWgsIt );
878  }
879  }
880 
881  return transformations;
882 }
883 
884 QString QgsCoordinateTransform::stripDatumTransform( const QString& proj4 )
885 {
886  QStringList parameterSplit = proj4.split( '+', QString::SkipEmptyParts );
887  QString currentParameter;
888  QString newProjString;
889 
890  for ( int i = 0; i < parameterSplit.size(); ++i )
891  {
892  currentParameter = parameterSplit.at( i );
893  if ( !currentParameter.startsWith( "towgs84", Qt::CaseInsensitive )
894  && !currentParameter.startsWith( "nadgrids", Qt::CaseInsensitive ) )
895  {
896  newProjString.append( '+' );
897  newProjString.append( currentParameter );
898  newProjString.append( ' ' );
899  }
900  }
901  return newProjString;
902 }
903 
904 void QgsCoordinateTransform::searchDatumTransform( const QString& sql, QList< int >& transforms )
905 {
906  sqlite3* db;
907  int openResult = sqlite3_open_v2( QgsApplication::srsDbFilePath().toUtf8().constData(), &db, SQLITE_OPEN_READONLY, 0 );
908  if ( openResult != SQLITE_OK )
909  {
910  sqlite3_close( db );
911  return;
912  }
913 
914  sqlite3_stmt* stmt;
915  int prepareRes = sqlite3_prepare( db, sql.toAscii(), sql.size(), &stmt, nullptr );
916  if ( prepareRes != SQLITE_OK )
917  {
918  sqlite3_finalize( stmt );
919  sqlite3_close( db );
920  return;
921  }
922 
923  QString cOpCode;
924  while ( sqlite3_step( stmt ) == SQLITE_ROW )
925  {
926  cOpCode = reinterpret_cast< const char * >( sqlite3_column_text( stmt, 0 ) );
927  transforms.push_back( cOpCode.toInt() );
928  }
929  sqlite3_finalize( stmt );
930  sqlite3_close( db );
931 }
932 
934 {
935  QString transformString;
936 
937  sqlite3* db;
938  int openResult = sqlite3_open_v2( QgsApplication::srsDbFilePath().toUtf8().constData(), &db, SQLITE_OPEN_READONLY, 0 );
939  if ( openResult != SQLITE_OK )
940  {
941  sqlite3_close( db );
942  return transformString;
943  }
944 
945  sqlite3_stmt* stmt;
946  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 );
947  int prepareRes = sqlite3_prepare( db, sql.toAscii(), sql.size(), &stmt, nullptr );
948  if ( prepareRes != SQLITE_OK )
949  {
950  sqlite3_finalize( stmt );
951  sqlite3_close( db );
952  return transformString;
953  }
954 
955  if ( sqlite3_step( stmt ) == SQLITE_ROW )
956  {
957  //coord_op_methode_code
958  int methodCode = sqlite3_column_int( stmt, 0 );
959  if ( methodCode == 9615 ) //ntv2
960  {
961  transformString = "+nadgrids=" + QString( reinterpret_cast< const char * >( sqlite3_column_text( stmt, 1 ) ) );
962  }
963  else if ( methodCode == 9603 || methodCode == 9606 || methodCode == 9607 )
964  {
965  transformString += "+towgs84=";
966  double p1 = sqlite3_column_double( stmt, 1 );
967  double p2 = sqlite3_column_double( stmt, 2 );
968  double p3 = sqlite3_column_double( stmt, 3 );
969  double p4 = sqlite3_column_double( stmt, 4 );
970  double p5 = sqlite3_column_double( stmt, 5 );
971  double p6 = sqlite3_column_double( stmt, 6 );
972  double p7 = sqlite3_column_double( stmt, 7 );
973  if ( methodCode == 9603 ) //3 parameter transformation
974  {
975  transformString += QString( "%1,%2,%3" ).arg( p1 ).arg( p2 ).arg( p3 );
976  }
977  else //7 parameter transformation
978  {
979  transformString += QString( "%1,%2,%3,%4,%5,%6,%7" ).arg( p1 ).arg( p2 ).arg( p3 ).arg( p4 ).arg( p5 ).arg( p6 ).arg( p7 );
980  }
981  }
982  }
983 
984  sqlite3_finalize( stmt );
985  sqlite3_close( db );
986  return transformString;
987 }
988 
989 bool QgsCoordinateTransform::datumTransformCrsInfo( int datumTransform, int& epsgNr, QString& srcProjection, QString& dstProjection, QString &remarks, QString &scope, bool &preferred, bool &deprecated )
990 {
991  sqlite3* db;
992  int openResult = sqlite3_open_v2( QgsApplication::srsDbFilePath().toUtf8().constData(), &db, SQLITE_OPEN_READONLY, 0 );
993  if ( openResult != SQLITE_OK )
994  {
995  sqlite3_close( db );
996  return false;
997  }
998 
999  sqlite3_stmt* stmt;
1000  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 );
1001  int prepareRes = sqlite3_prepare( db, sql.toAscii(), sql.size(), &stmt, nullptr );
1002  if ( prepareRes != SQLITE_OK )
1003  {
1004  sqlite3_finalize( stmt );
1005  sqlite3_close( db );
1006  return false;
1007  }
1008 
1009  int srcCrsId, destCrsId;
1010  if ( sqlite3_step( stmt ) != SQLITE_ROW )
1011  {
1012  sqlite3_finalize( stmt );
1013  sqlite3_close( db );
1014  return false;
1015  }
1016 
1017  epsgNr = sqlite3_column_int( stmt, 0 );
1018  srcCrsId = sqlite3_column_int( stmt, 1 );
1019  destCrsId = sqlite3_column_int( stmt, 2 );
1020  remarks = QString::fromUtf8( reinterpret_cast< const char * >( sqlite3_column_text( stmt, 3 ) ) );
1021  scope = QString::fromUtf8( reinterpret_cast< const char * >( sqlite3_column_text( stmt, 4 ) ) );
1022  preferred = sqlite3_column_int( stmt, 5 ) != 0;
1023  deprecated = sqlite3_column_int( stmt, 6 ) != 0;
1024 
1026  srcCrs.createFromOgcWmsCrs( QString( "EPSG:%1" ).arg( srcCrsId ) );
1027  srcProjection = srcCrs.description();
1029  destCrs.createFromOgcWmsCrs( QString( "EPSG:%1" ).arg( destCrsId ) );
1030  dstProjection = destCrs.description();
1031 
1032  sqlite3_finalize( stmt );
1033  sqlite3_close( db );
1034  return true;
1035 }
1036 
1037 void QgsCoordinateTransform::addNullGridShifts( QString& srcProjString, QString& destProjString )
1038 {
1039  //if one transformation uses ntv2, the other one needs to be null grid shift
1040  if ( mDestinationDatumTransform == -1 && srcProjString.contains( "+nadgrids" ) ) //add null grid if source transformation is ntv2
1041  {
1042  destProjString += " [email protected]";
1043  return;
1044  }
1045  if ( mSourceDatumTransform == -1 && destProjString.contains( "+nadgrids" ) )
1046  {
1047  srcProjString += " [email protected]";
1048  return;
1049  }
1050 
1051  //add null shift grid for google mercator
1052  //(see e.g. http://trac.osgeo.org/proj/wiki/FAQ#ChangingEllipsoidWhycantIconvertfromWGS84toGoogleEarthVirtualGlobeMercator)
1053  if ( mSourceCRS.authid().compare( "EPSG:3857", Qt::CaseInsensitive ) == 0 && mSourceDatumTransform == -1 )
1054  {
1055  srcProjString += " [email protected]";
1056  }
1057  if ( mDestCRS.authid().compare( "EPSG:3857", Qt::CaseInsensitive ) == 0 && mDestinationDatumTransform == -1 )
1058  {
1059  destProjString += " [email protected]";
1060  }
1061 }
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)
const QgsCoordinateReferenceSystem & crsByAuthId(const QString &authid)
Returns the CRS for authid, e.g.
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.
bool createFromWkt(const QString &theWkt)
Set up this CRS using a WKT spatial ref sys definition.
QDomNode appendChild(const QDomNode &newChild)
void setXMaximum(double x)
Set the maximum x value.
Definition: qgsrectangle.h:172
void push_back(const T &value)
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)
bool createFromOgcWmsCrs(const QString &theCrs)
Set up this CRS from the given OGC CRS.
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)
void combineExtentWith(QgsRectangle *rect)
expand the rectangle so that covers both the original rectangle and the given rectangle ...
bool createFromSrsId(const long theSrsId)
Set up this CRS by fetching the appropriate information from the sqlite backend.
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.
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
QDomNode namedItem(const QString &name) const
struct sqlite3 sqlite3
bool contains(QChar ch, Qt::CaseSensitivity cs) const
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()
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