QGIS API Documentation  2.4.0-Chugiak
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups 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  }
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.
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 {
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, 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 {
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 {
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 {
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 {
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 {
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 {
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 
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 
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  bb_rect.combineExtentWith( x[i], y[i] );
597  }
598 
599  QgsDebugMsg( "Projected extent: " + bb_rect.toString() );
600 
601  if ( bb_rect.isEmpty() )
602  {
603  QgsDebugMsg( "Original extent: " + rect.toString() );
604  }
605 
606  return bb_rect;
607 }
608 
609 void QgsCoordinateTransform::transformCoords( const int& numPoints, double *x, double *y, double *z, TransformDirection direction ) const
610 {
611  // Refuse to transform the points if the srs's are invalid
612  if ( !mSourceCRS.isValid() )
613  {
614  QgsMessageLog::logMessage( tr( "The source spatial reference system (CRS) is not valid. "
615  "The coordinates can not be reprojected. The CRS is: %1" )
616  .arg( mSourceCRS.toProj4() ), tr( "CRS" ) );
617  return;
618  }
619  if ( !mDestCRS.isValid() )
620  {
621  QgsMessageLog::logMessage( tr( "The destination spatial reference system (CRS) is not valid. "
622  "The coordinates can not be reprojected. The CRS is: %1" ).arg( mDestCRS.toProj4() ), tr( "CRS" ) );
623  return;
624  }
625 
626 #ifdef COORDINATE_TRANSFORM_VERBOSE
627  double xorg = *x;
628  double yorg = *y;
629  QgsDebugMsg( QString( "[[[[[[ Number of points to transform: %1 ]]]]]]" ).arg( numPoints ) );
630 #endif
631 
632  // use proj4 to do the transform
633  QString dir;
634  // if the source/destination projection is lat/long, convert the points to radians
635  // prior to transforming
636  if (( pj_is_latlong( mDestinationProjection ) && ( direction == ReverseTransform ) )
637  || ( pj_is_latlong( mSourceProjection ) && ( direction == ForwardTransform ) ) )
638  {
639  for ( int i = 0; i < numPoints; ++i )
640  {
641  x[i] *= DEG_TO_RAD;
642  y[i] *= DEG_TO_RAD;
643  z[i] *= DEG_TO_RAD;
644  }
645 
646  }
647  int projResult;
648  if ( direction == ReverseTransform )
649  {
650  projResult = pj_transform( mDestinationProjection, mSourceProjection, numPoints, 0, x, y, z );
651  }
652  else
653  {
654  Q_ASSERT( mSourceProjection != 0 );
655  Q_ASSERT( mDestinationProjection != 0 );
656  projResult = pj_transform( mSourceProjection, mDestinationProjection, numPoints, 0, x, y, z );
657  }
658 
659  if ( projResult != 0 )
660  {
661  //something bad happened....
662  QString points;
663 
664  for ( int i = 0; i < numPoints; ++i )
665  {
666  if ( direction == ForwardTransform )
667  {
668  points += QString( "(%1, %2)\n" ).arg( x[i], 0, 'f' ).arg( y[i], 0, 'f' );
669  }
670  else
671  {
672  points += QString( "(%1, %2)\n" ).arg( x[i] * RAD_TO_DEG, 0, 'f' ).arg( y[i] * RAD_TO_DEG, 0, 'f' );
673  }
674  }
675 
676  dir = ( direction == ForwardTransform ) ? tr( "forward transform" ) : tr( "inverse transform" );
677 
678  QString msg = tr( "%1 of\n"
679  "%2"
680  "PROJ.4: %3 +to %4\n"
681  "Error: %5" )
682  .arg( dir )
683  .arg( points )
684  .arg( mSourceCRS.toProj4() ).arg( mDestCRS.toProj4() )
685  .arg( QString::fromUtf8( pj_strerrno( projResult ) ) );
686 
687  QgsDebugMsg( "Projection failed emitting invalid transform signal: " + msg );
688 
689  emit invalidTransformInput();
690 
691  QgsDebugMsg( "throwing exception" );
692 
693  throw QgsCsException( msg );
694  }
695 
696  // if the result is lat/long, convert the results from radians back
697  // to degrees
698  if (( pj_is_latlong( mDestinationProjection ) && ( direction == ForwardTransform ) )
699  || ( pj_is_latlong( mSourceProjection ) && ( direction == ReverseTransform ) ) )
700  {
701  for ( int i = 0; i < numPoints; ++i )
702  {
703  x[i] *= RAD_TO_DEG;
704  y[i] *= RAD_TO_DEG;
705  z[i] *= RAD_TO_DEG;
706  }
707  }
708 #ifdef COORDINATE_TRANSFORM_VERBOSE
709  QgsDebugMsg( QString( "[[[[[[ Projected %1, %2 to %3, %4 ]]]]]]" )
710  .arg( xorg, 0, 'g', 15 ).arg( yorg, 0, 'g', 15 )
711  .arg( *x, 0, 'g', 15 ).arg( *y, 0, 'g', 15 ) );
712 #endif
713 }
714 
715 bool QgsCoordinateTransform::readXML( QDomNode & theNode )
716 {
717 
718  QgsDebugMsg( "Reading Coordinate Transform from xml ------------------------!" );
719 
720  QDomNode mySrcNode = theNode.namedItem( "sourcesrs" );
721  mSourceCRS.readXML( mySrcNode );
722 
723  QDomNode myDestNode = theNode.namedItem( "destinationsrs" );
724  mDestCRS.readXML( myDestNode );
725 
726  mSourceDatumTransform = theNode.toElement().attribute( "sourceDatumTransform", "-1" ).toInt();
727  mDestinationDatumTransform = theNode.toElement().attribute( "destinationDatumTransform", "-1" ).toInt();
728 
729  initialise();
730 
731  return true;
732 }
733 
734 bool QgsCoordinateTransform::writeXML( QDomNode & theNode, QDomDocument & theDoc )
735 {
736  QDomElement myNodeElement = theNode.toElement();
737  QDomElement myTransformElement = theDoc.createElement( "coordinatetransform" );
738  myTransformElement.setAttribute( "sourceDatumTransform", QString::number( mSourceDatumTransform ) );
739  myTransformElement.setAttribute( "destinationDatumTransform", QString::number( mDestinationDatumTransform ) );
740 
741  QDomElement mySourceElement = theDoc.createElement( "sourcesrs" );
742  mSourceCRS.writeXML( mySourceElement, theDoc );
743  myTransformElement.appendChild( mySourceElement );
744 
745  QDomElement myDestElement = theDoc.createElement( "destinationsrs" );
746  mDestCRS.writeXML( myDestElement, theDoc );
747  myTransformElement.appendChild( myDestElement );
748 
749  myNodeElement.appendChild( myTransformElement );
750 
751  return true;
752 }
753 
754 const char *finder( const char *name )
755 {
756  QString proj;
757 #ifdef WIN32
758  proj = QApplication::applicationDirPath()
759  + "/share/proj/" + QString( name );
760 #else
761  Q_UNUSED( name );
762 #endif
763  return proj.toUtf8();
764 }
765 
767 {
768 #if 0
769  // Attention! It should be possible to set PROJ_LIB
770  // but it can happen that it was previously set by installer
771  // (version 0.7) and the old installation was deleted
772 
773  // Another problem: PROJ checks if pj_finder was set before
774  // PROJ_LIB environment variable. pj_finder is probably set in
775  // GRASS gproj library when plugin is loaded, consequently
776  // PROJ_LIB is ignored
777 
778  pj_set_finder( finder );
779 #endif
780 }
781 
783 {
784  QList< QList< int > > transformations;
785 
786  QString srcGeoId = srcCRS.geographicCRSAuthId();
787  QString destGeoId = destCRS.geographicCRSAuthId();
788 
789  if ( srcGeoId.isEmpty() || destGeoId.isEmpty() )
790  {
791  return transformations;
792  }
793 
794  QStringList srcSplit = srcGeoId.split( ":" );
795  QStringList destSplit = destGeoId.split( ":" );
796 
797  if ( srcSplit.size() < 2 || destSplit.size() < 2 )
798  {
799  return transformations;
800  }
801 
802  int srcAuthCode = srcSplit.at( 1 ).toInt();
803  int destAuthCode = destSplit.at( 1 ).toInt();
804 
805  if ( srcAuthCode == destAuthCode )
806  {
807  return transformations; //crs have the same datum
808  }
809 
810  QList<int> directTransforms;
811  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 ),
812  directTransforms );
813  QList<int> reverseDirectTransforms;
814  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 ),
815  reverseDirectTransforms );
816  QList<int> srcToWgs84;
817  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 ),
818  srcToWgs84 );
819  QList<int> destToWgs84;
820  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 ),
821  destToWgs84 );
822 
823  //add direct datum transformations
824  QList<int>::const_iterator directIt = directTransforms.constBegin();
825  for ( ; directIt != directTransforms.constEnd(); ++directIt )
826  {
827  transformations.push_back( QList<int>() << *directIt << -1 );
828  }
829 
830  //add direct datum transformations
831  directIt = reverseDirectTransforms.constBegin();
832  for ( ; directIt != reverseDirectTransforms.constEnd(); ++directIt )
833  {
834  transformations.push_back( QList<int>() << -1 << *directIt );
835  }
836 
837  QList<int>::const_iterator srcWgsIt = srcToWgs84.constBegin();
838  for ( ; srcWgsIt != srcToWgs84.constEnd(); ++srcWgsIt )
839  {
840  QList<int>::const_iterator dstWgsIt = destToWgs84.constBegin();
841  for ( ; dstWgsIt != destToWgs84.constEnd(); ++dstWgsIt )
842  {
843  transformations.push_back( QList<int>() << *srcWgsIt << *dstWgsIt );
844  }
845  }
846 
847  return transformations;
848 }
849 
850 QString QgsCoordinateTransform::stripDatumTransform( const QString& proj4 )
851 {
852  QStringList parameterSplit = proj4.split( "+", QString::SkipEmptyParts );
853  QString currentParameter;
854  QString newProjString;
855 
856  for ( int i = 0; i < parameterSplit.size(); ++i )
857  {
858  currentParameter = parameterSplit.at( i );
859  if ( !currentParameter.startsWith( "towgs84", Qt::CaseInsensitive )
860  && !currentParameter.startsWith( "nadgrids", Qt::CaseInsensitive ) )
861  {
862  newProjString.append( "+" );
863  newProjString.append( currentParameter );
864  newProjString.append( " " );
865  }
866  }
867  return newProjString;
868 }
869 
870 void QgsCoordinateTransform::searchDatumTransform( const QString& sql, QList< int >& transforms )
871 {
872  sqlite3* db;
873  int openResult = sqlite3_open( QgsApplication::srsDbFilePath().toUtf8().constData(), &db );
874  if ( openResult != SQLITE_OK )
875  {
876  sqlite3_close( db );
877  return;
878  }
879 
880  sqlite3_stmt* stmt;
881  int prepareRes = sqlite3_prepare( db, sql.toAscii(), sql.size(), &stmt, NULL );
882  if ( prepareRes != SQLITE_OK )
883  {
884  sqlite3_finalize( stmt ); sqlite3_close( db );
885  return;
886  }
887 
888  QString cOpCode;
889  while ( sqlite3_step( stmt ) == SQLITE_ROW )
890  {
891  cOpCode = ( const char * ) sqlite3_column_text( stmt, 0 );
892  transforms.push_back( cOpCode.toInt() );
893  }
894  sqlite3_finalize( stmt ); sqlite3_close( db );
895 }
896 
897 QString QgsCoordinateTransform::datumTransformString( int datumTransform )
898 {
899  QString transformString;
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 transformString;
907  }
908 
909  sqlite3_stmt* stmt;
910  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 );
911  int prepareRes = sqlite3_prepare( db, sql.toAscii(), sql.size(), &stmt, NULL );
912  if ( prepareRes != SQLITE_OK )
913  {
914  sqlite3_finalize( stmt ); sqlite3_close( db );
915  return transformString;
916  }
917 
918  if ( sqlite3_step( stmt ) == SQLITE_ROW )
919  {
920  //coord_op_methode_code
921  int methodCode = sqlite3_column_int( stmt, 0 );
922  if ( methodCode == 9615 ) //ntv2
923  {
924  transformString = "+nadgrids=" + QString(( const char * )sqlite3_column_text( stmt, 1 ) );
925  }
926  else if ( methodCode == 9603 || methodCode == 9606 || methodCode == 9607 )
927  {
928  transformString += "+towgs84=";
929  double p1 = sqlite3_column_double( stmt, 1 );
930  double p2 = sqlite3_column_double( stmt, 2 );
931  double p3 = sqlite3_column_double( stmt, 3 );
932  double p4 = sqlite3_column_double( stmt, 4 );
933  double p5 = sqlite3_column_double( stmt, 5 );
934  double p6 = sqlite3_column_double( stmt, 6 );
935  double p7 = sqlite3_column_double( stmt, 7 );
936  if ( methodCode == 9603 ) //3 parameter transformation
937  {
938  transformString += QString( "%1,%2,%3" ).arg( p1 ).arg( p2 ).arg( p3 );
939  }
940  else //7 parameter transformation
941  {
942  transformString += QString( "%1,%2,%3,%4,%5,%6,%7" ).arg( p1 ).arg( p2 ).arg( p3 ).arg( p4 ).arg( p5 ).arg( p6 ).arg( p7 );
943  }
944  }
945  }
946 
947  sqlite3_finalize( stmt ); sqlite3_close( db );
948  return transformString;
949 }
950 
951 bool QgsCoordinateTransform::datumTransformCrsInfo( int datumTransform, int& epsgNr, QString& srcProjection, QString& dstProjection, QString &remarks, QString &scope, bool &preferred, bool &deprecated )
952 {
953  sqlite3* db;
954  int openResult = sqlite3_open( QgsApplication::srsDbFilePath().toUtf8().constData(), &db );
955  if ( openResult != SQLITE_OK )
956  {
957  sqlite3_close( db );
958  return false;
959  }
960 
961  sqlite3_stmt* stmt;
962  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 );
963  int prepareRes = sqlite3_prepare( db, sql.toAscii(), sql.size(), &stmt, NULL );
964  if ( prepareRes != SQLITE_OK )
965  {
966  sqlite3_finalize( stmt ); sqlite3_close( db );
967  return false;
968  }
969 
970  int srcCrsId, destCrsId;
971  if ( sqlite3_step( stmt ) != SQLITE_ROW )
972  {
973  sqlite3_finalize( stmt );
974  sqlite3_close( db );
975  return false;
976  }
977 
978  epsgNr = sqlite3_column_int( stmt, 0 );
979  srcCrsId = sqlite3_column_int( stmt, 1 );
980  destCrsId = sqlite3_column_int( stmt, 2 );
981  remarks = QString::fromUtf8(( const char * ) sqlite3_column_text( stmt, 3 ) );
982  scope = QString::fromUtf8(( const char * ) sqlite3_column_text( stmt, 4 ) );
983  preferred = sqlite3_column_int( stmt, 5 ) != 0;
984  deprecated = sqlite3_column_int( stmt, 6 ) != 0;
985 
987  srcCrs.createFromOgcWmsCrs( QString( "EPSG:%1" ).arg( srcCrsId ) );
988  srcProjection = srcCrs.description();
990  destCrs.createFromOgcWmsCrs( QString( "EPSG:%1" ).arg( destCrsId ) );
991  dstProjection = destCrs.description();
992 
993  sqlite3_finalize( stmt );
994  sqlite3_close( db );
995  return true;
996 }
997 
998 void QgsCoordinateTransform::addNullGridShifts( QString& srcProjString, QString& destProjString )
999 {
1000  //if one transformation uses ntv2, the other one needs to be null grid shift
1001  if ( mDestinationDatumTransform == -1 && srcProjString.contains( "+nadgrids" ) ) //add null grid if source transformation is ntv2
1002  {
1003  destProjString += " +nadgrids=@null";
1004  return;
1005  }
1006  if ( mSourceDatumTransform == -1 && destProjString.contains( "+nadgrids" ) )
1007  {
1008  srcProjString += " +nadgrids=@null";
1009  return;
1010  }
1011 
1012  //add null shift grid for google mercator
1013  //(see e.g. http://trac.osgeo.org/proj/wiki/FAQ#ChangingEllipsoidWhycantIconvertfromWGS84toGoogleEarthVirtualGlobeMercator)
1014  if ( mSourceCRS.authid().compare( "EPSG:3857", Qt::CaseInsensitive ) == 0 && mSourceDatumTransform == -1 )
1015  {
1016  srcProjString += " +nadgrids=@null";
1017  }
1018  if ( mDestCRS.authid().compare( "EPSG:3857", Qt::CaseInsensitive ) == 0 && mDestinationDatumTransform == -1 )
1019  {
1020  destProjString += " +nadgrids=@null";
1021  }
1022 }
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 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:194
#define QgsDebugMsg(str)
Definition: qgslogger.h:36
QString geographicCRSAuthId() const
Returns auth id of related geographic CRS.
void setSourceCrs(const QgsCoordinateReferenceSystem &theCRS)
TransformDirection
Enum used to indicate the direction (forward or inverse) of the transform.
bool readXML(QDomNode &theNode)
bool createFromId(const long theId, CrsType theType=PostgisCrsId)
void initialise()
initialise is used to actually create the Transformer instance
double x() const
Definition: qgspoint.h:110
static QString stripDatumTransform(const QString &proj4)
Removes +nadgrids and +towgs84 from proj4 string.
void addNullGridShifts(QString &srcProjString, QString &destProjString)
In certain situations, null grid shifts have to be added to src / dst proj string.
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.
QgsCoordinateReferenceSystem mSourceCRS
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:199
double xMaximum() const
Get the x maximum value (right side of rectangle)
Definition: qgsrectangle.h:184
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:37
static void searchDatumTransform(const QString &sql, QList< int > &transforms)
static QString datumTransformString(int datumTransform)
QgsPoint transform(const QgsPoint p, TransformDirection direction=ForwardTransform) const
A class to represent a point geometry.
Definition: qgspoint.h:63
struct sqlite3 sqlite3
void setDestCRSID(long theCRSID)
QgsCoordinateReferenceSystem mDestCRS
void setDestCRS(const QgsCoordinateReferenceSystem &theCRS)
bool writeXML(QDomNode &theNode, QDomDocument &theDoc) const
QgsRectangle transformBoundingBox(const QgsRectangle theRect, TransformDirection direction=ForwardTransform) 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:118
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:204
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:189
static QgsCRSCache * instance()
Definition: qgscrscache.cpp:85
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:209
const char * finder(const char *name)
QgsCoordinateTransform * clone() const
QString toProj4() const
Get the Proj Proj4 string representation of this srs.
#define tr(sourceText)