Quantum GIS API Documentation  1.7.4
src/core/qgscoordinatetransform.cpp
Go to the documentation of this file.
00001 /***************************************************************************
00002                QgsCoordinateTransform.cpp  - Coordinate Transforms
00003                              -------------------
00004     begin                : Dec 2004
00005     copyright            : (C) 2004 Tim Sutton
00006     email                : tim at linfiniti.com
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  *                                                                         *
00011  *   This program is free software; you can redistribute it and/or modify  *
00012  *   it under the terms of the GNU General Public License as published by  *
00013  *   the Free Software Foundation; either version 2 of the License, or     *
00014  *   (at your option) any later version.                                   *
00015  *                                                                         *
00016  ***************************************************************************/
00017 /* $Id$ */
00018 #include "qgscoordinatetransform.h"
00019 #include "qgsmessageoutput.h"
00020 #include "qgslogger.h"
00021 
00022 //qt includes
00023 #include <QDomNode>
00024 #include <QDomElement>
00025 #include <QApplication>
00026 
00027 extern "C"
00028 {
00029 #include <proj_api.h>
00030 }
00031 
00032 // if defined shows all information about transform to stdout
00033 // #define COORDINATE_TRANSFORM_VERBOSE
00034 
00035 QgsCoordinateTransform::QgsCoordinateTransform()
00036     : QObject()
00037     , mInitialisedFlag( false )
00038     , mSourceProjection( 0 )
00039     , mDestinationProjection( 0 )
00040 {
00041   setFinder();
00042 }
00043 
00044 QgsCoordinateTransform::QgsCoordinateTransform( const QgsCoordinateReferenceSystem& source, const QgsCoordinateReferenceSystem& dest )
00045     : QObject()
00046     , mInitialisedFlag( false )
00047     , mSourceProjection( 0 )
00048     , mDestinationProjection( 0 )
00049 {
00050   setFinder();
00051   mSourceCRS = source;
00052   mDestCRS = dest;
00053   initialise();
00054 }
00055 
00056 QgsCoordinateTransform::QgsCoordinateTransform( long theSourceSrsId, long theDestSrsId )
00057     : QObject()
00058     , mInitialisedFlag( false )
00059     , mSourceCRS( theSourceSrsId, QgsCoordinateReferenceSystem::InternalCrsId )
00060     , mDestCRS( theDestSrsId, QgsCoordinateReferenceSystem::InternalCrsId )
00061     , mSourceProjection( 0 )
00062     , mDestinationProjection( 0 )
00063 {
00064   initialise();
00065 }
00066 
00067 QgsCoordinateTransform::QgsCoordinateTransform( QString theSourceCRS, QString theDestCRS )
00068     : QObject()
00069     , mInitialisedFlag( false )
00070     , mSourceProjection( 0 )
00071     , mDestinationProjection( 0 )
00072 {
00073   setFinder();
00074   mSourceCRS.createFromWkt( theSourceCRS );
00075   mDestCRS.createFromWkt( theDestCRS );
00076   // initialize the coordinate system data structures
00077   //XXX Who spells initialize initialise?
00078   //XXX A: Its the queen's english....
00079   //XXX  : Long live the queen! Lets get on with the initialisation...
00080   initialise();
00081 }
00082 
00083 QgsCoordinateTransform::QgsCoordinateTransform( long theSourceSrid,
00084     QString theDestWkt,
00085     QgsCoordinateReferenceSystem::CrsType theSourceCRSType )
00086     : QObject()
00087     , mInitialisedFlag( false )
00088     , mSourceProjection( 0 )
00089     , mDestinationProjection( 0 )
00090 {
00091   setFinder();
00092 
00093   mSourceCRS.createFromId( theSourceSrid, theSourceCRSType );
00094   mDestCRS.createFromWkt( theDestWkt );
00095   // initialize the coordinate system data structures
00096   //XXX Who spells initialize initialise?
00097   //XXX A: Its the queen's english....
00098   //XXX  : Long live the queen! Lets get on with the initialisation...
00099   initialise();
00100 }
00101 
00102 QgsCoordinateTransform::~QgsCoordinateTransform()
00103 {
00104   // free the proj objects
00105   if ( mSourceProjection )
00106   {
00107     pj_free( mSourceProjection );
00108   }
00109   if ( mDestinationProjection )
00110   {
00111     pj_free( mDestinationProjection );
00112   }
00113 }
00114 
00115 void QgsCoordinateTransform::setSourceCrs( const QgsCoordinateReferenceSystem& theCRS )
00116 {
00117   mSourceCRS = theCRS;
00118   initialise();
00119 }
00120 void QgsCoordinateTransform::setDestCRS( const QgsCoordinateReferenceSystem& theCRS )
00121 {
00122   mDestCRS = theCRS;
00123   initialise();
00124 }
00125 
00126 void QgsCoordinateTransform::setDestCRSID( long theCRSID )
00127 {
00129   mDestCRS.createFromSrsId( theCRSID );
00130   initialise();
00131 }
00132 
00133 // XXX This whole function is full of multiple return statements!!!
00134 // And probably shouldn't be a void
00135 void QgsCoordinateTransform::initialise()
00136 {
00137   // XXX Warning - multiple return paths in this block!!
00138   if ( !mSourceCRS.isValid() )
00139   {
00140     //mSourceCRS = defaultWkt;
00141     // Pass through with no projection since we have no idea what the layer
00142     // coordinates are and projecting them may not be appropriate
00143     mShortCircuit = true;
00144     QgsDebugMsg( "SourceCRS seemed invalid!" );
00145     return;
00146   }
00147 
00148   if ( !mDestCRS.isValid() )
00149   {
00150     //No destination projection is set so we set the default output projection to
00151     //be the same as input proj. This only happens on the first layer loaded
00152     //whatever that may be...
00153     mDestCRS.createFromOgcWmsCrs( mSourceCRS.authid() );
00154   }
00155 
00156   // init the projections (destination and source)
00157   mDestinationProjection = pj_init_plus( mDestCRS.toProj4().toUtf8() );
00158   mSourceProjection = pj_init_plus( mSourceCRS.toProj4().toUtf8() );
00159 
00160 #ifdef COORDINATE_TRANSFORM_VERBOSE
00161   QgsDebugMsg( "From proj : " + mSourceCRS.toProj4() );
00162   QgsDebugMsg( "To proj   : " + mDestCRS.toProj4() );
00163 #endif
00164 
00165   mInitialisedFlag = true;
00166   if ( !mDestinationProjection )
00167   {
00168     mInitialisedFlag = false;
00169   }
00170   if ( !mSourceProjection )
00171   {
00172     mInitialisedFlag = false;
00173   }
00174 #ifdef COORDINATE_TRANSFORM_VERBOSE
00175   if ( mInitialisedFlag )
00176   {
00177     QgsDebugMsg( "------------------------------------------------------------" );
00178     QgsDebugMsg( "The OGR Coordinate transformation for this layer was set to" );
00179     QgsLogger::debug<QgsCoordinateReferenceSystem>( "Input", mSourceCRS, __FILE__, __FUNCTION__, __LINE__ );
00180     QgsLogger::debug<QgsCoordinateReferenceSystem>( "Output", mDestCRS, __FILE__, __FUNCTION__, __LINE__ );
00181     QgsDebugMsg( "------------------------------------------------------------" );
00182   }
00183   else
00184   {
00185     QgsDebugMsg( "------------------------------------------------------------" );
00186     QgsDebugMsg( "The OGR Coordinate transformation FAILED TO INITIALISE!" );
00187     QgsDebugMsg( "------------------------------------------------------------" );
00188   }
00189 #else
00190   if ( !mInitialisedFlag )
00191   {
00192     QgsDebugMsg( "Coordinate transformation failed to initialize!" );
00193   }
00194 #endif
00195 
00196   //XXX todo overload == operator for QgsCoordinateReferenceSystem
00197   //at the moment srs.parameters contains the whole proj def...soon it wont...
00198   //if (mSourceCRS->toProj4() == mDestCRS->toProj4())
00199   if ( mSourceCRS == mDestCRS )
00200   {
00201     // If the source and destination projection are the same, set the short
00202     // circuit flag (no transform takes place)
00203     mShortCircuit = true;
00204     QgsDebugMsgLevel( "Source/Dest CRS equal, shortcircuit is set.", 3 );
00205   }
00206   else
00207   {
00208     // Transform must take place
00209     mShortCircuit = false;
00210     QgsDebugMsgLevel( "Source/Dest CRS UNequal, shortcircuit is NOt set.", 3 );
00211   }
00212 
00213 }
00214 
00215 //
00216 //
00217 // TRANSFORMERS BELOW THIS POINT .........
00218 //
00219 //
00220 //
00221 
00222 
00223 QgsPoint QgsCoordinateTransform::transform( const QgsPoint thePoint, TransformDirection direction ) const
00224 {
00225   if ( mShortCircuit || !mInitialisedFlag ) return thePoint;
00226   // transform x
00227   double x = thePoint.x();
00228   double y = thePoint.y();
00229   double z = 0.0;
00230   try
00231   {
00232     transformCoords( 1, &x, &y, &z, direction );
00233   }
00234   catch ( QgsCsException &cse )
00235   {
00236     // rethrow the exception
00237     QgsDebugMsg( "rethrowing exception" );
00238     throw cse;
00239   }
00240 
00241   return QgsPoint( x, y );
00242 }
00243 
00244 
00245 QgsPoint QgsCoordinateTransform::transform( const double theX, const double theY = 0, TransformDirection direction ) const
00246 {
00247   try
00248   {
00249     return transform( QgsPoint( theX, theY ), direction );
00250   }
00251   catch ( QgsCsException &cse )
00252   {
00253     // rethrow the exception
00254     QgsDebugMsg( "rethrowing exception" );
00255     throw cse;
00256   }
00257 }
00258 
00259 QgsRectangle QgsCoordinateTransform::transform( const QgsRectangle theRect, TransformDirection direction ) const
00260 {
00261   if ( mShortCircuit || !mInitialisedFlag ) return theRect;
00262   // transform x
00263   double x1 = theRect.xMinimum();
00264   double y1 = theRect.yMinimum();
00265   double x2 = theRect.xMaximum();
00266   double y2 = theRect.yMaximum();
00267 
00268   // Number of points to reproject------+
00269   //                                    |
00270   //                                    V
00271   try
00272   {
00273     double z = 0.0;
00274     transformCoords( 1, &x1, &y1, &z, direction );
00275     transformCoords( 1, &x2, &y2, &z, direction );
00276   }
00277   catch ( QgsCsException &cse )
00278   {
00279     // rethrow the exception
00280     QgsDebugMsg( "rethrowing exception" );
00281     throw cse;
00282   }
00283 
00284 #ifdef COORDINATE_TRANSFORM_VERBOSE
00285   QgsDebugMsg( "Rect projection..." );
00286   QgsLogger::debug( "Xmin : ", theRect.xMinimum(), 1, __FILE__, __FUNCTION__, __LINE__ );
00287   QgsLogger::debug( "-->", x1, 1, __FILE__, __FUNCTION__, __LINE__ );
00288   QgsLogger::debug( "Ymin : ", theRect.yMinimum(), 1, __FILE__, __FUNCTION__, __LINE__ );
00289   QgsLogger::debug( "-->", y1, 1, __FILE__, __FUNCTION__, __LINE__ );
00290   QgsLogger::debug( "Xmax : ", theRect.xMaximum(), 1, __FILE__, __FUNCTION__, __LINE__ );
00291   QgsLogger::debug( "-->", x2, 1, __FILE__, __FUNCTION__, __LINE__ );
00292   QgsLogger::debug( "Ymax : ", theRect.yMaximum(), 1, __FILE__, __FUNCTION__, __LINE__ );
00293   QgsLogger::debug( "-->", y2, 1, __FILE__, __FUNCTION__, __LINE__ );
00294 #endif
00295   return QgsRectangle( x1, y1, x2, y2 );
00296 }
00297 
00298 void QgsCoordinateTransform::transformInPlace( double& x, double& y, double& z,
00299     TransformDirection direction ) const
00300 {
00301   if ( mShortCircuit || !mInitialisedFlag )
00302     return;
00303 #ifdef QGISDEBUG
00304 // QgsDebugMsg(QString("Using transform in place %1 %2").arg(__FILE__).arg(__LINE__));
00305 #endif
00306   // transform x
00307   try
00308   {
00309     transformCoords( 1, &x, &y, &z, direction );
00310   }
00311   catch ( QgsCsException &cse )
00312   {
00313     // rethrow the exception
00314     QgsDebugMsg( "rethrowing exception" );
00315     throw cse;
00316   }
00317 }
00318 
00319 void QgsCoordinateTransform::transformInPlace( std::vector<double>& x,
00320     std::vector<double>& y, std::vector<double>& z,
00321     TransformDirection direction ) const
00322 {
00323   if ( mShortCircuit || !mInitialisedFlag )
00324     return;
00325 
00326   Q_ASSERT( x.size() == y.size() );
00327 
00328   // Apparently, if one has a std::vector, it is valid to use the
00329   // address of the first element in the vector as a pointer to an
00330   // array of the vectors data, and hence easily interface with code
00331   // that wants C-style arrays.
00332 
00333   try
00334   {
00335     transformCoords( x.size(), &x[0], &y[0], &z[0], direction );
00336   }
00337   catch ( QgsCsException &cse )
00338   {
00339     // rethrow the exception
00340     QgsDebugMsg( "rethrowing exception" );
00341     throw cse;
00342   }
00343 }
00344 
00345 
00346 QgsRectangle QgsCoordinateTransform::transformBoundingBox( const QgsRectangle rect, TransformDirection direction ) const
00347 {
00348   // Calculate the bounding box of a QgsRectangle in the source CRS
00349   // when projected to the destination CRS (or the inverse).
00350   // This is done by looking at a number of points spread evenly
00351   // across the rectangle
00352 
00353   if ( mShortCircuit || !mInitialisedFlag )
00354     return rect;
00355 
00356   if ( rect.isEmpty() )
00357   {
00358     QgsPoint p = transform( rect.xMinimum(), rect.yMinimum(), direction );
00359     return QgsRectangle( p, p );
00360   }
00361 
00362   static const int numP = 8;
00363 
00364   QgsRectangle bb_rect;
00365   bb_rect.setMinimal();
00366 
00367   // We're interfacing with C-style vectors in the
00368   // end, so let's do C-style vectors here too.
00369 
00370   double x[numP * numP];
00371   double y[numP * numP];
00372   double z[numP * numP];
00373 
00374   QgsDebugMsg( "Entering transformBoundingBox..." );
00375 
00376   // Populate the vectors
00377 
00378   double dx = rect.width()  / ( double )( numP - 1 );
00379   double dy = rect.height() / ( double )( numP - 1 );
00380 
00381   double pointY = rect.yMinimum();
00382 
00383   for ( int i = 0; i < numP ; i++ )
00384   {
00385 
00386     // Start at right edge
00387     double pointX = rect.xMinimum();
00388 
00389     for ( int j = 0; j < numP; j++ )
00390     {
00391       x[( i*numP ) + j] = pointX;
00392       y[( i*numP ) + j] = pointY;
00393       // and the height...
00394       z[( i*numP ) + j] = 0.0;
00395       // QgsDebugMsg(QString("BBox coord: (%1, %2)").arg(x[(i*numP) + j]).arg(y[(i*numP) + j]));
00396       pointX += dx;
00397     }
00398     pointY += dy;
00399   }
00400 
00401   // Do transformation. Any exception generated must
00402   // be handled in above layers.
00403   try
00404   {
00405     transformCoords( numP * numP, x, y, z, direction );
00406   }
00407   catch ( QgsCsException &cse )
00408   {
00409     // rethrow the exception
00410     QgsDebugMsg( "rethrowing exception" );
00411     throw cse;
00412   }
00413 
00414   // Calculate the bounding box and use that for the extent
00415 
00416   for ( int i = 0; i < numP * numP; i++ )
00417   {
00418     if ( qIsFinite( x[i] ) && qIsFinite( y[i] ) )
00419       bb_rect.combineExtentWith( x[i], y[i] );
00420   }
00421 
00422   QgsDebugMsg( "Projected extent: " + QString(( bb_rect.toString() ).toLocal8Bit().data() ) );
00423 
00424   return bb_rect;
00425 }
00426 
00427 void QgsCoordinateTransform::transformCoords( const int& numPoints, double *x, double *y, double *z, TransformDirection direction ) const
00428 {
00429   // Refuse to transform the points if the srs's are invalid
00430   if ( !mSourceCRS.isValid() )
00431   {
00432     QgsLogger::critical( tr( "The source spatial reference system (CRS) is not valid. " )
00433                          + tr( "The coordinates can not be reprojected."
00434                                " The CRS is: %1" )
00435                          .arg( mSourceCRS.toProj4() ) );
00436     return;
00437   }
00438   if ( !mDestCRS.isValid() )
00439   {
00440     QgsLogger::critical( tr( "The destination spatial reference system (CRS) is not valid. " )
00441                          + tr( "The coordinates can not be reprojected."
00442                                " The CRS is: %1" ).arg( mDestCRS.toProj4() ) );
00443     return;
00444   }
00445 
00446 #ifdef COORDINATE_TRANSFORM_VERBOSE
00447   double xorg = *x;
00448   double yorg = *y;
00449   QgsDebugMsg( QString( "[[[[[[ Number of points to transform: %1 ]]]]]]" ).arg( numPoints ) );
00450 #endif
00451 
00452   // use proj4 to do the transform
00453   QString dir;
00454   // if the source/destination projection is lat/long, convert the points to radians
00455   // prior to transforming
00456   if (( pj_is_latlong( mDestinationProjection ) && ( direction == ReverseTransform ) )
00457       || ( pj_is_latlong( mSourceProjection ) && ( direction == ForwardTransform ) ) )
00458   {
00459     for ( int i = 0; i < numPoints; ++i )
00460     {
00461       x[i] *= DEG_TO_RAD;
00462       y[i] *= DEG_TO_RAD;
00463       z[i] *= DEG_TO_RAD;
00464     }
00465 
00466   }
00467   int projResult;
00468   if ( direction == ReverseTransform )
00469   {
00470     projResult = pj_transform( mDestinationProjection, mSourceProjection, numPoints, 0, x, y, z );
00471     dir = tr( "inverse transform" );
00472   }
00473   else
00474   {
00475     Q_ASSERT( mSourceProjection != 0 );
00476     Q_ASSERT( mDestinationProjection != 0 );
00477     projResult = pj_transform( mSourceProjection, mDestinationProjection, numPoints, 0, x, y, z );
00478     dir = tr( "forward transform" );
00479   }
00480 
00481   if ( projResult != 0 )
00482   {
00483     //something bad happened....
00484     QString points;
00485 
00486     for ( int i = 0; i < numPoints; ++i )
00487     {
00488       if ( direction == ForwardTransform )
00489       {
00490         points += QString( "(%1, %2)\n" ).arg( x[i] ).arg( y[i] );
00491       }
00492       else
00493       {
00494         points += QString( "(%1, %2)\n" ).arg( x[i] * RAD_TO_DEG ).arg( y[i] * RAD_TO_DEG );
00495       }
00496     }
00497 
00498     QString msg = tr( "%1 of\n%2\nfailed with error: %3\n" )
00499                   .arg( dir )
00500                   .arg( points )
00501                   .arg( QString::fromUtf8( pj_strerrno( projResult ) ) );
00502 
00503     QgsDebugMsg( "Projection failed emitting invalid transform signal: " + msg );
00504 
00505     emit invalidTransformInput();
00506 
00507     QgsDebugMsg( "throwing exception" );
00508 
00509     throw QgsCsException( msg );
00510   }
00511 
00512   // if the result is lat/long, convert the results from radians back
00513   // to degrees
00514   if (( pj_is_latlong( mDestinationProjection ) && ( direction == ForwardTransform ) )
00515       || ( pj_is_latlong( mSourceProjection ) && ( direction == ReverseTransform ) ) )
00516   {
00517     for ( int i = 0; i < numPoints; ++i )
00518     {
00519       x[i] *= RAD_TO_DEG;
00520       y[i] *= RAD_TO_DEG;
00521       z[i] *= RAD_TO_DEG;
00522     }
00523   }
00524 #ifdef COORDINATE_TRANSFORM_VERBOSE
00525   QgsDebugMsg( QString( "[[[[[[ Projected %1, %2 to %3, %4 ]]]]]]" )
00526                .arg( xorg, 0, 'g', 15 ).arg( yorg, 0, 'g', 15 )
00527                .arg( *x, 0, 'g', 15 ).arg( *y, 0, 'g', 15 ) );
00528 #endif
00529 }
00530 
00531 bool QgsCoordinateTransform::readXML( QDomNode & theNode )
00532 {
00533 
00534   QgsDebugMsg( "Reading Coordinate Transform from xml ------------------------!" );
00535 
00536   QDomNode mySrcNode = theNode.namedItem( "sourcesrs" );
00537   mSourceCRS.readXML( mySrcNode );
00538 
00539   QDomNode myDestNode = theNode.namedItem( "destinationsrs" );
00540   mDestCRS.readXML( myDestNode );
00541 
00542   initialise();
00543 
00544   return true;
00545 }
00546 
00547 bool QgsCoordinateTransform::writeXML( QDomNode & theNode, QDomDocument & theDoc )
00548 {
00549   QDomElement myNodeElement = theNode.toElement();
00550   QDomElement myTransformElement = theDoc.createElement( "coordinatetransform" );
00551 
00552   QDomElement mySourceElement = theDoc.createElement( "sourcesrs" );
00553   mSourceCRS.writeXML( mySourceElement, theDoc );
00554   myTransformElement.appendChild( mySourceElement );
00555 
00556   QDomElement myDestElement = theDoc.createElement( "destinationsrs" );
00557   mDestCRS.writeXML( myDestElement, theDoc );
00558   myTransformElement.appendChild( myDestElement );
00559 
00560   myNodeElement.appendChild( myTransformElement );
00561 
00562   return true;
00563 }
00564 
00565 const char *finder( const char *name )
00566 {
00567   QString proj;
00568 #ifdef WIN32
00569   proj = QApplication::applicationDirPath()
00570          + "/share/proj/" + QString( name );
00571 #endif
00572   return proj.toUtf8();
00573 }
00574 
00575 void QgsCoordinateTransform::setFinder()
00576 {
00577 #if 0
00578   // Attention! It should be possible to set PROJ_LIB
00579   // but it can happen that it was previously set by installer
00580   // (version 0.7) and the old installation was deleted
00581 
00582   // Another problem: PROJ checks if pj_finder was set before
00583   // PROJ_LIB environment variable. pj_finder is probably set in
00584   // GRASS gproj library when plugin is loaded, consequently
00585   // PROJ_LIB is ignored
00586 
00587   pj_set_finder( finder );
00588 #endif
00589 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines