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