QGIS API Documentation  2.9.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  , mShortCircuit( false )
43  , mInitialisedFlag( false )
44  , mSourceProjection( 0 )
45  , mDestinationProjection( 0 )
46  , mSourceDatumTransform( -1 )
47  , mDestinationDatumTransform( -1 )
48 {
49  setFinder();
50 }
51 
53  : QObject()
54  , mShortCircuit( false )
55  , mInitialisedFlag( false )
56  , mSourceProjection( 0 )
57  , mDestinationProjection( 0 )
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( 0 )
73  , mDestinationProjection( 0 )
74  , mSourceDatumTransform( -1 )
75  , mDestinationDatumTransform( -1 )
76 {
77  initialise();
78 }
79 
80 QgsCoordinateTransform::QgsCoordinateTransform( QString theSourceCRS, QString theDestCRS )
81  : QObject()
82  , mInitialisedFlag( false )
83  , mSourceProjection( 0 )
84  , mDestinationProjection( 0 )
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 initialisation...
95  initialise();
96 }
97 
99  QString theDestWkt,
100  QgsCoordinateReferenceSystem::CrsType theSourceCRSType )
101  : QObject()
102  , mInitialisedFlag( false )
103  , mSourceProjection( 0 )
104  , mDestinationProjection( 0 )
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 initialisation...
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::transformPolygon( QPolygonF& poly, TransformDirection direction ) const
373 {
374  if ( mShortCircuit || !mInitialisedFlag )
375  {
376  return;
377  }
378 
379  //create x, y arrays
380  int nVertices = poly.size();
381 
382  QVector<double> x( nVertices );
383  QVector<double> y( nVertices );
384  QVector<double> z( nVertices );
385 
386  for ( int i = 0; i < nVertices; ++i )
387  {
388  const QPointF& pt = poly.at( i );
389  x[i] = pt.x();
390  y[i] = pt.y();
391  z[i] = 0;
392  }
393 
394  try
395  {
396  transformCoords( nVertices, x.data(), y.data(), z.data(), direction );
397  }
398  catch ( const QgsCsException & )
399  {
400  // rethrow the exception
401  QgsDebugMsg( "rethrowing exception" );
402  throw;
403  }
404 
405  for ( int i = 0; i < nVertices; ++i )
406  {
407  QPointF& pt = poly[i];
408  pt.rx() = x[i];
409  pt.ry() = y[i];
410  }
411 }
412 
414  QVector<double>& x, QVector<double>& y, QVector<double>& z,
415  TransformDirection direction ) const
416 {
417  if ( mShortCircuit || !mInitialisedFlag )
418  return;
419 
420  Q_ASSERT( x.size() == y.size() );
421 
422  // Apparently, if one has a std::vector, it is valid to use the
423  // address of the first element in the vector as a pointer to an
424  // array of the vectors data, and hence easily interface with code
425  // that wants C-style arrays.
426 
427  try
428  {
429  transformCoords( x.size(), &x[0], &y[0], &z[0], direction );
430  }
431  catch ( const QgsCsException & )
432  {
433  // rethrow the exception
434  QgsDebugMsg( "rethrowing exception" );
435  throw;
436  }
437 }
438 
439 #ifdef QT_ARCH_ARM
440 void QgsCoordinateTransform::transformInPlace( qreal& x, qreal& y, double& z,
441  TransformDirection direction ) const
442 {
443  double xd = ( double ) x, yd = ( double ) y;
444  transformInPlace( xd, yd, z, direction );
445  x = xd;
446  y = yd;
447 }
448 #endif
449 
450 #ifdef ANDROID
451 void QgsCoordinateTransform::transformInPlace( float& x, float& y, float& z,
452  TransformDirection direction ) const
453 {
454  if ( mShortCircuit || !mInitialisedFlag )
455  return;
456 #ifdef QGISDEBUG
457 // QgsDebugMsg(QString("Using transform in place %1 %2").arg(__FILE__).arg(__LINE__));
458 #endif
459  // transform x
460  try
461  {
462  double xd = x;
463  double yd = y;
464  double zd = z;
465  transformCoords( 1, &xd, &yd, &zd, direction );
466  x = xd;
467  y = yd;
468  z = zd;
469  }
470  catch ( QgsCsException & )
471  {
472  // rethrow the exception
473  QgsDebugMsg( "rethrowing exception" );
474  throw;
475  }
476 }
477 
479  QVector<float>& x, QVector<float>& y, QVector<float>& z,
480  TransformDirection direction ) const
481 {
482  if ( mShortCircuit || !mInitialisedFlag )
483  return;
484 
485  Q_ASSERT( x.size() == y.size() );
486 
487  // Apparently, if one has a std::vector, it is valid to use the
488  // address of the first element in the vector as a pointer to an
489  // array of the vectors data, and hence easily interface with code
490  // that wants C-style arrays.
491 
492  try
493  {
494  //copy everything to double vectors since proj needs double
495  int vectorSize = x.size();
496  QVector<double> xd( x.size() );
497  QVector<double> yd( y.size() );
498  QVector<double> zd( z.size() );
499  for ( int i = 0; i < vectorSize; ++i )
500  {
501  xd[i] = x[i];
502  yd[i] = y[i];
503  zd[i] = z[i];
504  }
505  transformCoords( x.size(), &xd[0], &yd[0], &zd[0], direction );
506 
507  //copy back
508  for ( int i = 0; i < vectorSize; ++i )
509  {
510  x[i] = xd[i];
511  y[i] = yd[i];
512  z[i] = zd[i];
513  }
514  }
515  catch ( QgsCsException & )
516  {
517  // rethrow the exception
518  QgsDebugMsg( "rethrowing exception" );
519  throw;
520  }
521 }
522 #endif //ANDROID
523 
524 
525 QgsRectangle QgsCoordinateTransform::transformBoundingBox( const QgsRectangle &rect, TransformDirection direction, const bool handle180Crossover ) const
526 {
527  // Calculate the bounding box of a QgsRectangle in the source CRS
528  // when projected to the destination CRS (or the inverse).
529  // This is done by looking at a number of points spread evenly
530  // across the rectangle
531 
532  if ( mShortCircuit || !mInitialisedFlag )
533  return rect;
534 
535  if ( rect.isEmpty() )
536  {
537  QgsPoint p = transform( rect.xMinimum(), rect.yMinimum(), direction );
538  return QgsRectangle( p, p );
539  }
540 
541  static const int numP = 8;
542 
543  QgsRectangle bb_rect;
544  bb_rect.setMinimal();
545 
546  // We're interfacing with C-style vectors in the
547  // end, so let's do C-style vectors here too.
548 
549  double x[numP * numP];
550  double y[numP * numP];
551  double z[numP * numP];
552 
553  QgsDebugMsg( "Entering transformBoundingBox..." );
554 
555  // Populate the vectors
556 
557  double dx = rect.width() / ( double )( numP - 1 );
558  double dy = rect.height() / ( double )( numP - 1 );
559 
560  double pointY = rect.yMinimum();
561 
562  for ( int i = 0; i < numP ; i++ )
563  {
564 
565  // Start at right edge
566  double pointX = rect.xMinimum();
567 
568  for ( int j = 0; j < numP; j++ )
569  {
570  x[( i*numP ) + j] = pointX;
571  y[( i*numP ) + j] = pointY;
572  // and the height...
573  z[( i*numP ) + j] = 0.0;
574  // QgsDebugMsg(QString("BBox coord: (%1, %2)").arg(x[(i*numP) + j]).arg(y[(i*numP) + j]));
575  pointX += dx;
576  }
577  pointY += dy;
578  }
579 
580  // Do transformation. Any exception generated must
581  // be handled in above layers.
582  try
583  {
584  transformCoords( numP * numP, x, y, z, direction );
585  }
586  catch ( const QgsCsException & )
587  {
588  // rethrow the exception
589  QgsDebugMsg( "rethrowing exception" );
590  throw;
591  }
592 
593  // Calculate the bounding box and use that for the extent
594 
595  for ( int i = 0; i < numP * numP; i++ )
596  {
597  if ( !qIsFinite( x[i] ) || !qIsFinite( y[i] ) )
598  {
599  continue;
600  }
601 
602  if ( handle180Crossover )
603  {
604  //if crossing the date line, temporarily add 360 degrees to -ve longitudes
605  bb_rect.combineExtentWith( x[i] >= 0.0 ? x[i] : x[i] + 360.0, y[i] );
606  }
607  else
608  {
609  bb_rect.combineExtentWith( x[i], y[i] );
610  }
611  }
612 
613  if ( handle180Crossover )
614  {
615  //subtract temporary addition of 360 degrees from longitudes
616  if ( bb_rect.xMinimum() > 180.0 )
617  bb_rect.setXMinimum( bb_rect.xMinimum() - 360.0 );
618  if ( bb_rect.xMaximum() > 180.0 )
619  bb_rect.setXMaximum( bb_rect.xMaximum() - 360.0 );
620  }
621 
622  QgsDebugMsg( "Projected extent: " + bb_rect.toString() );
623 
624  if ( bb_rect.isEmpty() )
625  {
626  QgsDebugMsg( "Original extent: " + rect.toString() );
627  }
628 
629  return bb_rect;
630 }
631 
632 void QgsCoordinateTransform::transformCoords( const int& numPoints, double *x, double *y, double *z, TransformDirection direction ) const
633 {
634  // Refuse to transform the points if the srs's are invalid
635  if ( !mSourceCRS.isValid() )
636  {
637  QgsMessageLog::logMessage( tr( "The source spatial reference system (CRS) is not valid. "
638  "The coordinates can not be reprojected. The CRS is: %1" )
639  .arg( mSourceCRS.toProj4() ), tr( "CRS" ) );
640  return;
641  }
642  if ( !mDestCRS.isValid() )
643  {
644  QgsMessageLog::logMessage( tr( "The destination spatial reference system (CRS) is not valid. "
645  "The coordinates can not be reprojected. The CRS is: %1" ).arg( mDestCRS.toProj4() ), tr( "CRS" ) );
646  return;
647  }
648 
649 #ifdef COORDINATE_TRANSFORM_VERBOSE
650  double xorg = *x;
651  double yorg = *y;
652  QgsDebugMsg( QString( "[[[[[[ Number of points to transform: %1 ]]]]]]" ).arg( numPoints ) );
653 #endif
654 
655  // use proj4 to do the transform
656  QString dir;
657  // if the source/destination projection is lat/long, convert the points to radians
658  // prior to transforming
659  if (( pj_is_latlong( mDestinationProjection ) && ( direction == ReverseTransform ) )
660  || ( pj_is_latlong( mSourceProjection ) && ( direction == ForwardTransform ) ) )
661  {
662  for ( int i = 0; i < numPoints; ++i )
663  {
664  x[i] *= DEG_TO_RAD;
665  y[i] *= DEG_TO_RAD;
666  z[i] *= DEG_TO_RAD;
667  }
668 
669  }
670  int projResult;
671  if ( direction == ReverseTransform )
672  {
673  projResult = pj_transform( mDestinationProjection, mSourceProjection, numPoints, 0, x, y, z );
674  }
675  else
676  {
677  Q_ASSERT( mSourceProjection != 0 );
678  Q_ASSERT( mDestinationProjection != 0 );
679  projResult = pj_transform( mSourceProjection, mDestinationProjection, numPoints, 0, x, y, z );
680  }
681 
682  if ( projResult != 0 )
683  {
684  //something bad happened....
685  QString points;
686 
687  for ( int i = 0; i < numPoints; ++i )
688  {
689  if ( direction == ForwardTransform )
690  {
691  points += QString( "(%1, %2)\n" ).arg( x[i], 0, 'f' ).arg( y[i], 0, 'f' );
692  }
693  else
694  {
695  points += QString( "(%1, %2)\n" ).arg( x[i] * RAD_TO_DEG, 0, 'f' ).arg( y[i] * RAD_TO_DEG, 0, 'f' );
696  }
697  }
698 
699  dir = ( direction == ForwardTransform ) ? tr( "forward transform" ) : tr( "inverse transform" );
700 
701  char *srcdef = pj_get_def( mSourceProjection, 0 );
702  char *dstdef = pj_get_def( mDestinationProjection, 0 );
703 
704  QString msg = tr( "%1 of\n"
705  "%2"
706  "PROJ.4: %3 +to %4\n"
707  "Error: %5" )
708  .arg( dir )
709  .arg( points )
710  .arg( srcdef ).arg( dstdef )
711  .arg( QString::fromUtf8( pj_strerrno( projResult ) ) );
712 
713  pj_dalloc( srcdef );
714  pj_dalloc( dstdef );
715 
716  QgsDebugMsg( "Projection failed emitting invalid transform signal: " + msg );
717 
718  emit invalidTransformInput();
719 
720  QgsDebugMsg( "throwing exception" );
721 
722  throw QgsCsException( msg );
723  }
724 
725  // if the result is lat/long, convert the results from radians back
726  // to degrees
727  if (( pj_is_latlong( mDestinationProjection ) && ( direction == ForwardTransform ) )
728  || ( pj_is_latlong( mSourceProjection ) && ( direction == ReverseTransform ) ) )
729  {
730  for ( int i = 0; i < numPoints; ++i )
731  {
732  x[i] *= RAD_TO_DEG;
733  y[i] *= RAD_TO_DEG;
734  z[i] *= RAD_TO_DEG;
735  }
736  }
737 #ifdef COORDINATE_TRANSFORM_VERBOSE
738  QgsDebugMsg( QString( "[[[[[[ Projected %1, %2 to %3, %4 ]]]]]]" )
739  .arg( xorg, 0, 'g', 15 ).arg( yorg, 0, 'g', 15 )
740  .arg( *x, 0, 'g', 15 ).arg( *y, 0, 'g', 15 ) );
741 #endif
742 }
743 
744 bool QgsCoordinateTransform::readXML( QDomNode & theNode )
745 {
746 
747  QgsDebugMsg( "Reading Coordinate Transform from xml ------------------------!" );
748 
749  QDomNode mySrcNode = theNode.namedItem( "sourcesrs" );
750  mSourceCRS.readXML( mySrcNode );
751 
752  QDomNode myDestNode = theNode.namedItem( "destinationsrs" );
753  mDestCRS.readXML( myDestNode );
754 
755  mSourceDatumTransform = theNode.toElement().attribute( "sourceDatumTransform", "-1" ).toInt();
756  mDestinationDatumTransform = theNode.toElement().attribute( "destinationDatumTransform", "-1" ).toInt();
757 
758  initialise();
759 
760  return true;
761 }
762 
763 bool QgsCoordinateTransform::writeXML( QDomNode & theNode, QDomDocument & theDoc )
764 {
765  QDomElement myNodeElement = theNode.toElement();
766  QDomElement myTransformElement = theDoc.createElement( "coordinatetransform" );
767  myTransformElement.setAttribute( "sourceDatumTransform", QString::number( mSourceDatumTransform ) );
768  myTransformElement.setAttribute( "destinationDatumTransform", QString::number( mDestinationDatumTransform ) );
769 
770  QDomElement mySourceElement = theDoc.createElement( "sourcesrs" );
771  mSourceCRS.writeXML( mySourceElement, theDoc );
772  myTransformElement.appendChild( mySourceElement );
773 
774  QDomElement myDestElement = theDoc.createElement( "destinationsrs" );
775  mDestCRS.writeXML( myDestElement, theDoc );
776  myTransformElement.appendChild( myDestElement );
777 
778  myNodeElement.appendChild( myTransformElement );
779 
780  return true;
781 }
782 
783 const char *finder( const char *name )
784 {
785  QString proj;
786 #ifdef WIN32
787  proj = QApplication::applicationDirPath()
788  + "/share/proj/" + QString( name );
789 #else
790  Q_UNUSED( name );
791 #endif
792  return proj.toUtf8();
793 }
794 
795 void QgsCoordinateTransform::setFinder()
796 {
797 #if 0
798  // Attention! It should be possible to set PROJ_LIB
799  // but it can happen that it was previously set by installer
800  // (version 0.7) and the old installation was deleted
801 
802  // Another problem: PROJ checks if pj_finder was set before
803  // PROJ_LIB environment variable. pj_finder is probably set in
804  // GRASS gproj library when plugin is loaded, consequently
805  // PROJ_LIB is ignored
806 
807  pj_set_finder( finder );
808 #endif
809 }
810 
812 {
813  QList< QList< int > > transformations;
814 
815  QString srcGeoId = srcCRS.geographicCRSAuthId();
816  QString destGeoId = destCRS.geographicCRSAuthId();
817 
818  if ( srcGeoId.isEmpty() || destGeoId.isEmpty() )
819  {
820  return transformations;
821  }
822 
823  QStringList srcSplit = srcGeoId.split( ":" );
824  QStringList destSplit = destGeoId.split( ":" );
825 
826  if ( srcSplit.size() < 2 || destSplit.size() < 2 )
827  {
828  return transformations;
829  }
830 
831  int srcAuthCode = srcSplit.at( 1 ).toInt();
832  int destAuthCode = destSplit.at( 1 ).toInt();
833 
834  if ( srcAuthCode == destAuthCode )
835  {
836  return transformations; //crs have the same datum
837  }
838 
839  QList<int> directTransforms;
840  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 ),
841  directTransforms );
842  QList<int> reverseDirectTransforms;
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( destAuthCode ).arg( srcAuthCode ),
844  reverseDirectTransforms );
845  QList<int> srcToWgs84;
846  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 ),
847  srcToWgs84 );
848  QList<int> destToWgs84;
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( destAuthCode ).arg( 4326 ),
850  destToWgs84 );
851 
852  //add direct datum transformations
853  QList<int>::const_iterator directIt = directTransforms.constBegin();
854  for ( ; directIt != directTransforms.constEnd(); ++directIt )
855  {
856  transformations.push_back( QList<int>() << *directIt << -1 );
857  }
858 
859  //add direct datum transformations
860  directIt = reverseDirectTransforms.constBegin();
861  for ( ; directIt != reverseDirectTransforms.constEnd(); ++directIt )
862  {
863  transformations.push_back( QList<int>() << -1 << *directIt );
864  }
865 
866  QList<int>::const_iterator srcWgsIt = srcToWgs84.constBegin();
867  for ( ; srcWgsIt != srcToWgs84.constEnd(); ++srcWgsIt )
868  {
869  QList<int>::const_iterator dstWgsIt = destToWgs84.constBegin();
870  for ( ; dstWgsIt != destToWgs84.constEnd(); ++dstWgsIt )
871  {
872  transformations.push_back( QList<int>() << *srcWgsIt << *dstWgsIt );
873  }
874  }
875 
876  return transformations;
877 }
878 
879 QString QgsCoordinateTransform::stripDatumTransform( const QString& proj4 )
880 {
881  QStringList parameterSplit = proj4.split( "+", QString::SkipEmptyParts );
882  QString currentParameter;
883  QString newProjString;
884 
885  for ( int i = 0; i < parameterSplit.size(); ++i )
886  {
887  currentParameter = parameterSplit.at( i );
888  if ( !currentParameter.startsWith( "towgs84", Qt::CaseInsensitive )
889  && !currentParameter.startsWith( "nadgrids", Qt::CaseInsensitive ) )
890  {
891  newProjString.append( "+" );
892  newProjString.append( currentParameter );
893  newProjString.append( " " );
894  }
895  }
896  return newProjString;
897 }
898 
899 void QgsCoordinateTransform::searchDatumTransform( const QString& sql, QList< int >& transforms )
900 {
901  sqlite3* db;
902  int openResult = sqlite3_open( QgsApplication::srsDbFilePath().toUtf8().constData(), &db );
903  if ( openResult != SQLITE_OK )
904  {
905  sqlite3_close( db );
906  return;
907  }
908 
909  sqlite3_stmt* stmt;
910  int prepareRes = sqlite3_prepare( db, sql.toAscii(), sql.size(), &stmt, NULL );
911  if ( prepareRes != SQLITE_OK )
912  {
913  sqlite3_finalize( stmt ); sqlite3_close( db );
914  return;
915  }
916 
917  QString cOpCode;
918  while ( sqlite3_step( stmt ) == SQLITE_ROW )
919  {
920  cOpCode = ( const char * ) sqlite3_column_text( stmt, 0 );
921  transforms.push_back( cOpCode.toInt() );
922  }
923  sqlite3_finalize( stmt ); sqlite3_close( db );
924 }
925 
926 QString QgsCoordinateTransform::datumTransformString( int datumTransform )
927 {
928  QString transformString;
929 
930  sqlite3* db;
931  int openResult = sqlite3_open( QgsApplication::srsDbFilePath().toUtf8().constData(), &db );
932  if ( openResult != SQLITE_OK )
933  {
934  sqlite3_close( db );
935  return transformString;
936  }
937 
938  sqlite3_stmt* stmt;
939  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 );
940  int prepareRes = sqlite3_prepare( db, sql.toAscii(), sql.size(), &stmt, NULL );
941  if ( prepareRes != SQLITE_OK )
942  {
943  sqlite3_finalize( stmt ); sqlite3_close( db );
944  return transformString;
945  }
946 
947  if ( sqlite3_step( stmt ) == SQLITE_ROW )
948  {
949  //coord_op_methode_code
950  int methodCode = sqlite3_column_int( stmt, 0 );
951  if ( methodCode == 9615 ) //ntv2
952  {
953  transformString = "+nadgrids=" + QString(( const char * )sqlite3_column_text( stmt, 1 ) );
954  }
955  else if ( methodCode == 9603 || methodCode == 9606 || methodCode == 9607 )
956  {
957  transformString += "+towgs84=";
958  double p1 = sqlite3_column_double( stmt, 1 );
959  double p2 = sqlite3_column_double( stmt, 2 );
960  double p3 = sqlite3_column_double( stmt, 3 );
961  double p4 = sqlite3_column_double( stmt, 4 );
962  double p5 = sqlite3_column_double( stmt, 5 );
963  double p6 = sqlite3_column_double( stmt, 6 );
964  double p7 = sqlite3_column_double( stmt, 7 );
965  if ( methodCode == 9603 ) //3 parameter transformation
966  {
967  transformString += QString( "%1,%2,%3" ).arg( p1 ).arg( p2 ).arg( p3 );
968  }
969  else //7 parameter transformation
970  {
971  transformString += QString( "%1,%2,%3,%4,%5,%6,%7" ).arg( p1 ).arg( p2 ).arg( p3 ).arg( p4 ).arg( p5 ).arg( p6 ).arg( p7 );
972  }
973  }
974  }
975 
976  sqlite3_finalize( stmt ); sqlite3_close( db );
977  return transformString;
978 }
979 
980 bool QgsCoordinateTransform::datumTransformCrsInfo( int datumTransform, int& epsgNr, QString& srcProjection, QString& dstProjection, QString &remarks, QString &scope, bool &preferred, bool &deprecated )
981 {
982  sqlite3* db;
983  int openResult = sqlite3_open( QgsApplication::srsDbFilePath().toUtf8().constData(), &db );
984  if ( openResult != SQLITE_OK )
985  {
986  sqlite3_close( db );
987  return false;
988  }
989 
990  sqlite3_stmt* stmt;
991  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 );
992  int prepareRes = sqlite3_prepare( db, sql.toAscii(), sql.size(), &stmt, NULL );
993  if ( prepareRes != SQLITE_OK )
994  {
995  sqlite3_finalize( stmt ); sqlite3_close( db );
996  return false;
997  }
998 
999  int srcCrsId, destCrsId;
1000  if ( sqlite3_step( stmt ) != SQLITE_ROW )
1001  {
1002  sqlite3_finalize( stmt );
1003  sqlite3_close( db );
1004  return false;
1005  }
1006 
1007  epsgNr = sqlite3_column_int( stmt, 0 );
1008  srcCrsId = sqlite3_column_int( stmt, 1 );
1009  destCrsId = sqlite3_column_int( stmt, 2 );
1010  remarks = QString::fromUtf8(( const char * ) sqlite3_column_text( stmt, 3 ) );
1011  scope = QString::fromUtf8(( const char * ) sqlite3_column_text( stmt, 4 ) );
1012  preferred = sqlite3_column_int( stmt, 5 ) != 0;
1013  deprecated = sqlite3_column_int( stmt, 6 ) != 0;
1014 
1016  srcCrs.createFromOgcWmsCrs( QString( "EPSG:%1" ).arg( srcCrsId ) );
1017  srcProjection = srcCrs.description();
1019  destCrs.createFromOgcWmsCrs( QString( "EPSG:%1" ).arg( destCrsId ) );
1020  dstProjection = destCrs.description();
1021 
1022  sqlite3_finalize( stmt );
1023  sqlite3_close( db );
1024  return true;
1025 }
1026 
1027 void QgsCoordinateTransform::addNullGridShifts( QString& srcProjString, QString& destProjString )
1028 {
1029  //if one transformation uses ntv2, the other one needs to be null grid shift
1030  if ( mDestinationDatumTransform == -1 && srcProjString.contains( "+nadgrids" ) ) //add null grid if source transformation is ntv2
1031  {
1032  destProjString += " +nadgrids=@null";
1033  return;
1034  }
1035  if ( mSourceDatumTransform == -1 && destProjString.contains( "+nadgrids" ) )
1036  {
1037  srcProjString += " +nadgrids=@null";
1038  return;
1039  }
1040 
1041  //add null shift grid for google mercator
1042  //(see e.g. http://trac.osgeo.org/proj/wiki/FAQ#ChangingEllipsoidWhycantIconvertfromWGS84toGoogleEarthVirtualGlobeMercator)
1043  if ( mSourceCRS.authid().compare( "EPSG:3857", Qt::CaseInsensitive ) == 0 && mSourceDatumTransform == -1 )
1044  {
1045  srcProjString += " +nadgrids=@null";
1046  }
1047  if ( mDestCRS.authid().compare( "EPSG:3857", Qt::CaseInsensitive ) == 0 && mDestinationDatumTransform == -1 )
1048  {
1049  destProjString += " +nadgrids=@null";
1050  }
1051 }
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:82
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)