Quantum GIS API Documentation
1.7.4
|
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 }