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