Quantum GIS API Documentation  1.7.5-Wroclaw
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 /* $Id$ */
18 #include "qgscoordinatetransform.h"
19 #include "qgsmessageoutput.h"
20 #include "qgslogger.h"
21 
22 //qt includes
23 #include <QDomNode>
24 #include <QDomElement>
25 #include <QApplication>
26 
27 extern "C"
28 {
29 #include <proj_api.h>
30 }
31 
32 // if defined shows all information about transform to stdout
33 // #define COORDINATE_TRANSFORM_VERBOSE
34 
36  : QObject()
37  , mInitialisedFlag( false )
38  , mSourceProjection( 0 )
39  , mDestinationProjection( 0 )
40 {
41  setFinder();
42 }
43 
45  : QObject()
46  , mInitialisedFlag( false )
47  , mSourceProjection( 0 )
48  , mDestinationProjection( 0 )
49 {
50  setFinder();
51  mSourceCRS = source;
52  mDestCRS = dest;
53  initialise();
54 }
55 
56 QgsCoordinateTransform::QgsCoordinateTransform( long theSourceSrsId, long theDestSrsId )
57  : QObject()
58  , mInitialisedFlag( false )
59  , mSourceCRS( theSourceSrsId, QgsCoordinateReferenceSystem::InternalCrsId )
60  , mDestCRS( theDestSrsId, QgsCoordinateReferenceSystem::InternalCrsId )
61  , mSourceProjection( 0 )
62  , mDestinationProjection( 0 )
63 {
64  initialise();
65 }
66 
67 QgsCoordinateTransform::QgsCoordinateTransform( QString theSourceCRS, QString theDestCRS )
68  : QObject()
69  , mInitialisedFlag( false )
70  , mSourceProjection( 0 )
71  , mDestinationProjection( 0 )
72 {
73  setFinder();
74  mSourceCRS.createFromWkt( theSourceCRS );
75  mDestCRS.createFromWkt( theDestCRS );
76  // initialize the coordinate system data structures
77  //XXX Who spells initialize initialise?
78  //XXX A: Its the queen's english....
79  //XXX : Long live the queen! Lets get on with the initialisation...
80  initialise();
81 }
82 
84  QString theDestWkt,
85  QgsCoordinateReferenceSystem::CrsType theSourceCRSType )
86  : QObject()
87  , mInitialisedFlag( false )
88  , mSourceProjection( 0 )
89  , mDestinationProjection( 0 )
90 {
91  setFinder();
92 
93  mSourceCRS.createFromId( theSourceSrid, theSourceCRSType );
94  mDestCRS.createFromWkt( theDestWkt );
95  // initialize the coordinate system data structures
96  //XXX Who spells initialize initialise?
97  //XXX A: Its the queen's english....
98  //XXX : Long live the queen! Lets get on with the initialisation...
99  initialise();
100 }
101 
103 {
104  // free the proj objects
105  if ( mSourceProjection )
106  {
107  pj_free( mSourceProjection );
108  }
110  {
111  pj_free( mDestinationProjection );
112  }
113 }
114 
116 {
117  mSourceCRS = theCRS;
118  initialise();
119 }
121 {
122  mDestCRS = theCRS;
123  initialise();
124 }
125 
127 {
129  mDestCRS.createFromSrsId( theCRSID );
130  initialise();
131 }
132 
133 // XXX This whole function is full of multiple return statements!!!
134 // And probably shouldn't be a void
136 {
137  // XXX Warning - multiple return paths in this block!!
138  if ( !mSourceCRS.isValid() )
139  {
140  //mSourceCRS = defaultWkt;
141  // Pass through with no projection since we have no idea what the layer
142  // coordinates are and projecting them may not be appropriate
143  mShortCircuit = true;
144  QgsDebugMsg( "SourceCRS seemed invalid!" );
145  return;
146  }
147 
148  if ( !mDestCRS.isValid() )
149  {
150  //No destination projection is set so we set the default output projection to
151  //be the same as input proj. This only happens on the first layer loaded
152  //whatever that may be...
154  }
155 
156  // init the projections (destination and source)
157  mDestinationProjection = pj_init_plus( mDestCRS.toProj4().toUtf8() );
158  mSourceProjection = pj_init_plus( mSourceCRS.toProj4().toUtf8() );
159 
160 #ifdef COORDINATE_TRANSFORM_VERBOSE
161  QgsDebugMsg( "From proj : " + mSourceCRS.toProj4() );
162  QgsDebugMsg( "To proj : " + mDestCRS.toProj4() );
163 #endif
164 
165  mInitialisedFlag = true;
166  if ( !mDestinationProjection )
167  {
168  mInitialisedFlag = false;
169  }
170  if ( !mSourceProjection )
171  {
172  mInitialisedFlag = false;
173  }
174 #ifdef COORDINATE_TRANSFORM_VERBOSE
175  if ( mInitialisedFlag )
176  {
177  QgsDebugMsg( "------------------------------------------------------------" );
178  QgsDebugMsg( "The OGR Coordinate transformation for this layer was set to" );
179  QgsLogger::debug<QgsCoordinateReferenceSystem>( "Input", mSourceCRS, __FILE__, __FUNCTION__, __LINE__ );
180  QgsLogger::debug<QgsCoordinateReferenceSystem>( "Output", mDestCRS, __FILE__, __FUNCTION__, __LINE__ );
181  QgsDebugMsg( "------------------------------------------------------------" );
182  }
183  else
184  {
185  QgsDebugMsg( "------------------------------------------------------------" );
186  QgsDebugMsg( "The OGR Coordinate transformation FAILED TO INITIALISE!" );
187  QgsDebugMsg( "------------------------------------------------------------" );
188  }
189 #else
190  if ( !mInitialisedFlag )
191  {
192  QgsDebugMsg( "Coordinate transformation failed to initialize!" );
193  }
194 #endif
195 
196  //XXX todo overload == operator for QgsCoordinateReferenceSystem
197  //at the moment srs.parameters contains the whole proj def...soon it wont...
198  //if (mSourceCRS->toProj4() == mDestCRS->toProj4())
199  if ( mSourceCRS == mDestCRS )
200  {
201  // If the source and destination projection are the same, set the short
202  // circuit flag (no transform takes place)
203  mShortCircuit = true;
204  QgsDebugMsgLevel( "Source/Dest CRS equal, shortcircuit is set.", 3 );
205  }
206  else
207  {
208  // Transform must take place
209  mShortCircuit = false;
210  QgsDebugMsgLevel( "Source/Dest CRS UNequal, shortcircuit is NOt set.", 3 );
211  }
212 
213 }
214 
215 //
216 //
217 // TRANSFORMERS BELOW THIS POINT .........
218 //
219 //
220 //
221 
222 
224 {
225  if ( mShortCircuit || !mInitialisedFlag ) return thePoint;
226  // transform x
227  double x = thePoint.x();
228  double y = thePoint.y();
229  double z = 0.0;
230  try
231  {
232  transformCoords( 1, &x, &y, &z, direction );
233  }
234  catch ( QgsCsException &cse )
235  {
236  // rethrow the exception
237  QgsDebugMsg( "rethrowing exception" );
238  throw cse;
239  }
240 
241  return QgsPoint( x, y );
242 }
243 
244 
245 QgsPoint QgsCoordinateTransform::transform( const double theX, const double theY = 0, TransformDirection direction ) const
246 {
247  try
248  {
249  return transform( QgsPoint( theX, theY ), direction );
250  }
251  catch ( QgsCsException &cse )
252  {
253  // rethrow the exception
254  QgsDebugMsg( "rethrowing exception" );
255  throw cse;
256  }
257 }
258 
260 {
261  if ( mShortCircuit || !mInitialisedFlag ) return theRect;
262  // transform x
263  double x1 = theRect.xMinimum();
264  double y1 = theRect.yMinimum();
265  double x2 = theRect.xMaximum();
266  double y2 = theRect.yMaximum();
267 
268  // Number of points to reproject------+
269  // |
270  // V
271  try
272  {
273  double z = 0.0;
274  transformCoords( 1, &x1, &y1, &z, direction );
275  transformCoords( 1, &x2, &y2, &z, direction );
276  }
277  catch ( QgsCsException &cse )
278  {
279  // rethrow the exception
280  QgsDebugMsg( "rethrowing exception" );
281  throw cse;
282  }
283 
284 #ifdef COORDINATE_TRANSFORM_VERBOSE
285  QgsDebugMsg( "Rect projection..." );
286  QgsLogger::debug( "Xmin : ", theRect.xMinimum(), 1, __FILE__, __FUNCTION__, __LINE__ );
287  QgsLogger::debug( "-->", x1, 1, __FILE__, __FUNCTION__, __LINE__ );
288  QgsLogger::debug( "Ymin : ", theRect.yMinimum(), 1, __FILE__, __FUNCTION__, __LINE__ );
289  QgsLogger::debug( "-->", y1, 1, __FILE__, __FUNCTION__, __LINE__ );
290  QgsLogger::debug( "Xmax : ", theRect.xMaximum(), 1, __FILE__, __FUNCTION__, __LINE__ );
291  QgsLogger::debug( "-->", x2, 1, __FILE__, __FUNCTION__, __LINE__ );
292  QgsLogger::debug( "Ymax : ", theRect.yMaximum(), 1, __FILE__, __FUNCTION__, __LINE__ );
293  QgsLogger::debug( "-->", y2, 1, __FILE__, __FUNCTION__, __LINE__ );
294 #endif
295  return QgsRectangle( x1, y1, x2, y2 );
296 }
297 
298 void QgsCoordinateTransform::transformInPlace( double& x, double& y, double& z,
299  TransformDirection direction ) const
300 {
302  return;
303 #ifdef QGISDEBUG
304 // QgsDebugMsg(QString("Using transform in place %1 %2").arg(__FILE__).arg(__LINE__));
305 #endif
306  // transform x
307  try
308  {
309  transformCoords( 1, &x, &y, &z, direction );
310  }
311  catch ( QgsCsException &cse )
312  {
313  // rethrow the exception
314  QgsDebugMsg( "rethrowing exception" );
315  throw cse;
316  }
317 }
318 
319 void QgsCoordinateTransform::transformInPlace( std::vector<double>& x,
320  std::vector<double>& y, std::vector<double>& z,
321  TransformDirection direction ) const
322 {
324  return;
325 
326  Q_ASSERT( x.size() == y.size() );
327 
328  // Apparently, if one has a std::vector, it is valid to use the
329  // address of the first element in the vector as a pointer to an
330  // array of the vectors data, and hence easily interface with code
331  // that wants C-style arrays.
332 
333  try
334  {
335  transformCoords( x.size(), &x[0], &y[0], &z[0], direction );
336  }
337  catch ( QgsCsException &cse )
338  {
339  // rethrow the exception
340  QgsDebugMsg( "rethrowing exception" );
341  throw cse;
342  }
343 }
344 
345 
347 {
348  // Calculate the bounding box of a QgsRectangle in the source CRS
349  // when projected to the destination CRS (or the inverse).
350  // This is done by looking at a number of points spread evenly
351  // across the rectangle
352 
354  return rect;
355 
356  if ( rect.isEmpty() )
357  {
358  QgsPoint p = transform( rect.xMinimum(), rect.yMinimum(), direction );
359  return QgsRectangle( p, p );
360  }
361 
362  static const int numP = 8;
363 
364  QgsRectangle bb_rect;
365  bb_rect.setMinimal();
366 
367  // We're interfacing with C-style vectors in the
368  // end, so let's do C-style vectors here too.
369 
370  double x[numP * numP];
371  double y[numP * numP];
372  double z[numP * numP];
373 
374  QgsDebugMsg( "Entering transformBoundingBox..." );
375 
376  // Populate the vectors
377 
378  double dx = rect.width() / ( double )( numP - 1 );
379  double dy = rect.height() / ( double )( numP - 1 );
380 
381  double pointY = rect.yMinimum();
382 
383  for ( int i = 0; i < numP ; i++ )
384  {
385 
386  // Start at right edge
387  double pointX = rect.xMinimum();
388 
389  for ( int j = 0; j < numP; j++ )
390  {
391  x[( i*numP ) + j] = pointX;
392  y[( i*numP ) + j] = pointY;
393  // and the height...
394  z[( i*numP ) + j] = 0.0;
395  // QgsDebugMsg(QString("BBox coord: (%1, %2)").arg(x[(i*numP) + j]).arg(y[(i*numP) + j]));
396  pointX += dx;
397  }
398  pointY += dy;
399  }
400 
401  // Do transformation. Any exception generated must
402  // be handled in above layers.
403  try
404  {
405  transformCoords( numP * numP, x, y, z, direction );
406  }
407  catch ( QgsCsException &cse )
408  {
409  // rethrow the exception
410  QgsDebugMsg( "rethrowing exception" );
411  throw cse;
412  }
413 
414  // Calculate the bounding box and use that for the extent
415 
416  for ( int i = 0; i < numP * numP; i++ )
417  {
418  if ( qIsFinite( x[i] ) && qIsFinite( y[i] ) )
419  bb_rect.combineExtentWith( x[i], y[i] );
420  }
421 
422  QgsDebugMsg( "Projected extent: " + QString(( bb_rect.toString() ).toLocal8Bit().data() ) );
423 
424  return bb_rect;
425 }
426 
427 void QgsCoordinateTransform::transformCoords( const int& numPoints, double *x, double *y, double *z, TransformDirection direction ) const
428 {
429  // Refuse to transform the points if the srs's are invalid
430  if ( !mSourceCRS.isValid() )
431  {
432  QgsLogger::critical( tr( "The source spatial reference system (CRS) is not valid. " )
433  + tr( "The coordinates can not be reprojected."
434  " The CRS is: %1" )
435  .arg( mSourceCRS.toProj4() ) );
436  return;
437  }
438  if ( !mDestCRS.isValid() )
439  {
440  QgsLogger::critical( tr( "The destination spatial reference system (CRS) is not valid. " )
441  + tr( "The coordinates can not be reprojected."
442  " The CRS is: %1" ).arg( mDestCRS.toProj4() ) );
443  return;
444  }
445 
446 #ifdef COORDINATE_TRANSFORM_VERBOSE
447  double xorg = *x;
448  double yorg = *y;
449  QgsDebugMsg( QString( "[[[[[[ Number of points to transform: %1 ]]]]]]" ).arg( numPoints ) );
450 #endif
451 
452  // use proj4 to do the transform
453  QString dir;
454  // if the source/destination projection is lat/long, convert the points to radians
455  // prior to transforming
456  if (( pj_is_latlong( mDestinationProjection ) && ( direction == ReverseTransform ) )
457  || ( pj_is_latlong( mSourceProjection ) && ( direction == ForwardTransform ) ) )
458  {
459  for ( int i = 0; i < numPoints; ++i )
460  {
461  x[i] *= DEG_TO_RAD;
462  y[i] *= DEG_TO_RAD;
463  z[i] *= DEG_TO_RAD;
464  }
465 
466  }
467  int projResult;
468  if ( direction == ReverseTransform )
469  {
470  projResult = pj_transform( mDestinationProjection, mSourceProjection, numPoints, 0, x, y, z );
471  dir = tr( "inverse transform" );
472  }
473  else
474  {
475  Q_ASSERT( mSourceProjection != 0 );
476  Q_ASSERT( mDestinationProjection != 0 );
477  projResult = pj_transform( mSourceProjection, mDestinationProjection, numPoints, 0, x, y, z );
478  dir = tr( "forward transform" );
479  }
480 
481  if ( projResult != 0 )
482  {
483  //something bad happened....
484  QString points;
485 
486  for ( int i = 0; i < numPoints; ++i )
487  {
488  if ( direction == ForwardTransform )
489  {
490  points += QString( "(%1, %2)\n" ).arg( x[i] ).arg( y[i] );
491  }
492  else
493  {
494  points += QString( "(%1, %2)\n" ).arg( x[i] * RAD_TO_DEG ).arg( y[i] * RAD_TO_DEG );
495  }
496  }
497 
498  QString msg = tr( "%1 of\n%2\nfailed with error: %3\n" )
499  .arg( dir )
500  .arg( points )
501  .arg( QString::fromUtf8( pj_strerrno( projResult ) ) );
502 
503  QgsDebugMsg( "Projection failed emitting invalid transform signal: " + msg );
504 
505  emit invalidTransformInput();
506 
507  QgsDebugMsg( "throwing exception" );
508 
509  throw QgsCsException( msg );
510  }
511 
512  // if the result is lat/long, convert the results from radians back
513  // to degrees
514  if (( pj_is_latlong( mDestinationProjection ) && ( direction == ForwardTransform ) )
515  || ( pj_is_latlong( mSourceProjection ) && ( direction == ReverseTransform ) ) )
516  {
517  for ( int i = 0; i < numPoints; ++i )
518  {
519  x[i] *= RAD_TO_DEG;
520  y[i] *= RAD_TO_DEG;
521  z[i] *= RAD_TO_DEG;
522  }
523  }
524 #ifdef COORDINATE_TRANSFORM_VERBOSE
525  QgsDebugMsg( QString( "[[[[[[ Projected %1, %2 to %3, %4 ]]]]]]" )
526  .arg( xorg, 0, 'g', 15 ).arg( yorg, 0, 'g', 15 )
527  .arg( *x, 0, 'g', 15 ).arg( *y, 0, 'g', 15 ) );
528 #endif
529 }
530 
531 bool QgsCoordinateTransform::readXML( QDomNode & theNode )
532 {
533 
534  QgsDebugMsg( "Reading Coordinate Transform from xml ------------------------!" );
535 
536  QDomNode mySrcNode = theNode.namedItem( "sourcesrs" );
537  mSourceCRS.readXML( mySrcNode );
538 
539  QDomNode myDestNode = theNode.namedItem( "destinationsrs" );
540  mDestCRS.readXML( myDestNode );
541 
542  initialise();
543 
544  return true;
545 }
546 
547 bool QgsCoordinateTransform::writeXML( QDomNode & theNode, QDomDocument & theDoc )
548 {
549  QDomElement myNodeElement = theNode.toElement();
550  QDomElement myTransformElement = theDoc.createElement( "coordinatetransform" );
551 
552  QDomElement mySourceElement = theDoc.createElement( "sourcesrs" );
553  mSourceCRS.writeXML( mySourceElement, theDoc );
554  myTransformElement.appendChild( mySourceElement );
555 
556  QDomElement myDestElement = theDoc.createElement( "destinationsrs" );
557  mDestCRS.writeXML( myDestElement, theDoc );
558  myTransformElement.appendChild( myDestElement );
559 
560  myNodeElement.appendChild( myTransformElement );
561 
562  return true;
563 }
564 
565 const char *finder( const char *name )
566 {
567  QString proj;
568 #ifdef WIN32
569  proj = QApplication::applicationDirPath()
570  + "/share/proj/" + QString( name );
571 #endif
572  return proj.toUtf8();
573 }
574 
576 {
577 #if 0
578  // Attention! It should be possible to set PROJ_LIB
579  // but it can happen that it was previously set by installer
580  // (version 0.7) and the old installation was deleted
581 
582  // Another problem: PROJ checks if pj_finder was set before
583  // PROJ_LIB environment variable. pj_finder is probably set in
584  // GRASS gproj library when plugin is loaded, consequently
585  // PROJ_LIB is ignored
586 
587  pj_set_finder( finder );
588 #endif
589 }