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