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