00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include "qgscoordinatetransform.h"
00019 #include "qgsmessageoutput.h"
00020 #include "qgslogger.h"
00021
00022
00023 #include <QDomNode>
00024 #include <QDomElement>
00025 #include <QApplication>
00026 #include "qgslogger.h"
00027
00028 extern "C"
00029 {
00030 #include <proj_api.h>
00031 }
00032
00033
00034
00035
00036 QgsCoordinateTransform::QgsCoordinateTransform( ) : QObject(), mSourceCRS(), mDestCRS()
00037 {
00038 setFinder();
00039 }
00040
00041 QgsCoordinateTransform::QgsCoordinateTransform( const QgsCoordinateReferenceSystem& source,
00042 const QgsCoordinateReferenceSystem& dest )
00043 {
00044 setFinder();
00045
00046 mSourceCRS = source;
00047 mDestCRS = dest;
00048 initialise();
00049 }
00050
00051 QgsCoordinateTransform::QgsCoordinateTransform( long theSourceSrsId, long theDestSrsId )
00052 : mSourceCRS( theSourceSrsId, QgsCoordinateReferenceSystem::InternalCrsId ),
00053 mDestCRS( theDestSrsId, QgsCoordinateReferenceSystem::InternalCrsId )
00054 {
00055 initialise();
00056 }
00057
00058 QgsCoordinateTransform::QgsCoordinateTransform( QString theSourceCRS, QString theDestCRS ) : QObject()
00059
00060 {
00061 setFinder();
00062 mSourceCRS.createFromWkt( theSourceCRS );
00063 mDestCRS.createFromWkt( theDestCRS );
00064
00065
00066
00067
00068 initialise();
00069 }
00070
00071 QgsCoordinateTransform::QgsCoordinateTransform( long theSourceSrid,
00072 QString theDestWkt,
00073 QgsCoordinateReferenceSystem::CrsType theSourceCRSType ): QObject()
00074 {
00075 setFinder();
00076
00077 mSourceCRS.createFromId( theSourceSrid, theSourceCRSType );
00078 mDestCRS.createFromWkt( theDestWkt );
00079
00080
00081
00082
00083 initialise();
00084 }
00085
00086 QgsCoordinateTransform::~QgsCoordinateTransform()
00087 {
00088
00089 if ( mSourceProjection != 0 )
00090 {
00091 pj_free( mSourceProjection );
00092 }
00093 if ( mDestinationProjection != 0 )
00094 {
00095 pj_free( mDestinationProjection );
00096 }
00097 }
00098
00099 void QgsCoordinateTransform::setSourceCrs( const QgsCoordinateReferenceSystem& theCRS )
00100 {
00101 mSourceCRS = theCRS;
00102 initialise();
00103 }
00104 void QgsCoordinateTransform::setDestCRS( const QgsCoordinateReferenceSystem& theCRS )
00105 {
00106 mDestCRS = theCRS;
00107 initialise();
00108 }
00109
00110
00111 void QgsCoordinateTransform::setDestCRSID( long theCRSID )
00112 {
00114 mDestCRS.createFromSrsId( theCRSID );
00115 initialise();
00116 }
00117
00118
00119
00120 void QgsCoordinateTransform::initialise()
00121 {
00122
00123 mInitialisedFlag = false;
00124 mSourceProjection = NULL;
00125 mDestinationProjection = NULL;
00126
00127
00128 if ( !mSourceCRS.isValid() )
00129 {
00130
00131
00132
00133 mShortCircuit = true;
00134 QgsDebugMsg( "SourceCRS seemed invalid!" );
00135 return;
00136 }
00137
00138 if ( !mDestCRS.isValid() )
00139 {
00140
00141
00142
00143 mDestCRS.createFromProj4( mSourceCRS.toProj4() );
00144 }
00145
00146
00147 mDestinationProjection = pj_init_plus( mDestCRS.toProj4().toUtf8() );
00148 mSourceProjection = pj_init_plus( mSourceCRS.toProj4().toUtf8() );
00149
00150 #ifdef COORDINATE_TRANSFORM_VERBOSE
00151 QgsDebugMsg( "From proj : " + mSourceCRS.toProj4() );
00152 QgsDebugMsg( "To proj : " + mDestCRS.toProj4() );
00153 #endif
00154
00155 mInitialisedFlag = true;
00156 if ( mDestinationProjection == NULL )
00157 {
00158 mInitialisedFlag = false;
00159 }
00160 if ( mSourceProjection == NULL )
00161 {
00162 mInitialisedFlag = false;
00163 }
00164 #ifdef COORDINATE_TRANSFORM_VERBOSE
00165 if ( mInitialisedFlag )
00166 {
00167 QgsDebugMsg( "------------------------------------------------------------" );
00168 QgsDebugMsg( "The OGR Coordinate transformation for this layer was set to" );
00169 QgsLogger::debug<QgsCoordinateReferenceSystem>( "Input", mSourceCRS, __FILE__, __FUNCTION__, __LINE__ );
00170 QgsLogger::debug<QgsCoordinateReferenceSystem>( "Output", mDestCRS, __FILE__, __FUNCTION__, __LINE__ );
00171 QgsDebugMsg( "------------------------------------------------------------" );
00172 }
00173 else
00174 {
00175 QgsDebugMsg( "------------------------------------------------------------" );
00176 QgsDebugMsg( "The OGR Coordinate transformation FAILED TO INITIALISE!" );
00177 QgsDebugMsg( "------------------------------------------------------------" );
00178 }
00179 #else
00180 if ( !mInitialisedFlag )
00181 {
00182 QgsDebugMsg( "Coordinate transformation failed to initialize!" );
00183 }
00184 #endif
00185
00186
00187
00188
00189 if ( mSourceCRS == mDestCRS )
00190 {
00191
00192
00193 mShortCircuit = true;
00194 QgsDebugMsgLevel( "Source/Dest CRS equal, shortcircuit is set.", 3 );
00195 }
00196 else
00197 {
00198
00199 mShortCircuit = false;
00200 QgsDebugMsgLevel( "Source/Dest CRS UNequal, shortcircuit is NOt set.", 3 );
00201 }
00202
00203 }
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213 QgsPoint QgsCoordinateTransform::transform( const QgsPoint thePoint, TransformDirection direction ) const
00214 {
00215 if ( mShortCircuit || !mInitialisedFlag ) return thePoint;
00216
00217 double x = thePoint.x();
00218 double y = thePoint.y();
00219 double z = 0.0;
00220 try
00221 {
00222 transformCoords( 1, &x, &y, &z, direction );
00223 }
00224 catch ( QgsCsException &cse )
00225 {
00226
00227 QgsDebugMsg( "rethrowing exception" );
00228 throw cse;
00229 }
00230
00231 return QgsPoint( x, y );
00232 }
00233
00234
00235 QgsPoint QgsCoordinateTransform::transform( const double theX, const double theY = 0, TransformDirection direction ) const
00236 {
00237 try
00238 {
00239 return transform( QgsPoint( theX, theY ), direction );
00240 }
00241 catch ( QgsCsException &cse )
00242 {
00243
00244 QgsDebugMsg( "rethrowing exception" );
00245 throw cse;
00246 }
00247 }
00248
00249 QgsRectangle QgsCoordinateTransform::transform( const QgsRectangle theRect, TransformDirection direction ) const
00250 {
00251 if ( mShortCircuit || !mInitialisedFlag ) return theRect;
00252
00253 double x1 = theRect.xMinimum();
00254 double y1 = theRect.yMinimum();
00255 double x2 = theRect.xMaximum();
00256 double y2 = theRect.yMaximum();
00257
00258
00259
00260
00261 try
00262 {
00263 double z = 0.0;
00264 transformCoords( 1, &x1, &y1, &z, direction );
00265 transformCoords( 1, &x2, &y2, &z, direction );
00266 }
00267 catch ( QgsCsException &cse )
00268 {
00269
00270 QgsDebugMsg( "rethrowing exception" );
00271 throw cse;
00272 }
00273
00274 #ifdef COORDINATE_TRANSFORM_VERBOSE
00275 QgsDebugMsg( "Rect projection..." );
00276 QgsLogger::debug( "Xmin : ", theRect.xMinimum(), 1, __FILE__, __FUNCTION__, __LINE__ );
00277 QgsLogger::debug( "-->", x1, 1, __FILE__, __FUNCTION__, __LINE__ );
00278 QgsLogger::debug( "Ymin : ", theRect.yMinimum(), 1, __FILE__, __FUNCTION__, __LINE__ );
00279 QgsLogger::debug( "-->", y1, 1, __FILE__, __FUNCTION__, __LINE__ );
00280 QgsLogger::debug( "Xmax : ", theRect.xMaximum(), 1, __FILE__, __FUNCTION__, __LINE__ );
00281 QgsLogger::debug( "-->", x2, 1, __FILE__, __FUNCTION__, __LINE__ );
00282 QgsLogger::debug( "Ymax : ", theRect.yMaximum(), 1, __FILE__, __FUNCTION__, __LINE__ );
00283 QgsLogger::debug( "-->", y2, 1, __FILE__, __FUNCTION__, __LINE__ );
00284 #endif
00285 return QgsRectangle( x1, y1, x2, y2 );
00286 }
00287
00288 void QgsCoordinateTransform::transformInPlace( double& x, double& y, double& z,
00289 TransformDirection direction ) const
00290 {
00291 if ( mShortCircuit || !mInitialisedFlag )
00292 return;
00293 #ifdef QGISDEBUG
00294
00295 #endif
00296
00297 try
00298 {
00299 transformCoords( 1, &x, &y, &z, direction );
00300 }
00301 catch ( QgsCsException &cse )
00302 {
00303
00304 QgsDebugMsg( "rethrowing exception" );
00305 throw cse;
00306 }
00307 }
00308
00309 void QgsCoordinateTransform::transformInPlace( std::vector<double>& x,
00310 std::vector<double>& y, std::vector<double>& z,
00311 TransformDirection direction ) const
00312 {
00313 if ( mShortCircuit || !mInitialisedFlag )
00314 return;
00315
00316 Q_ASSERT( x.size() == y.size() );
00317
00318
00319
00320
00321
00322
00323 try
00324 {
00325 transformCoords( x.size(), &x[0], &y[0], &z[0], direction );
00326 }
00327 catch ( QgsCsException &cse )
00328 {
00329
00330 QgsDebugMsg( "rethrowing exception" );
00331 throw cse;
00332 }
00333 }
00334
00335
00336 QgsRectangle QgsCoordinateTransform::transformBoundingBox( const QgsRectangle rect, TransformDirection direction ) const
00337 {
00338
00339
00340
00341
00342
00343 if ( mShortCircuit || !mInitialisedFlag )
00344 return rect;
00345
00346 static const int numP = 8;
00347
00348 QgsRectangle bb_rect;
00349 bb_rect.setMinimal();
00350
00351
00352
00353
00354 double x[numP * numP];
00355 double y[numP * numP];
00356 double z[numP * numP];
00357
00358 QgsDebugMsg( "Entering transformBoundingBox..." );
00359
00360
00361
00362 double dx = rect.width() / ( double )( numP - 1 );
00363 double dy = rect.height() / ( double )( numP - 1 );
00364
00365 double pointY = rect.yMinimum();
00366
00367 for ( int i = 0; i < numP ; i++ )
00368 {
00369
00370
00371 double pointX = rect.xMinimum();
00372
00373 for ( int j = 0; j < numP; j++ )
00374 {
00375 x[( i*numP ) + j] = pointX;
00376 y[( i*numP ) + j] = pointY;
00377
00378 z[( i*numP ) + j] = 0.0;
00379
00380 pointX += dx;
00381 }
00382 pointY += dy;
00383 }
00384
00385
00386
00387 try
00388 {
00389 transformCoords( numP * numP, x, y, z, direction );
00390 }
00391 catch ( QgsCsException &cse )
00392 {
00393
00394 QgsDebugMsg( "rethrowing exception" );
00395 throw cse;
00396 }
00397
00398
00399
00400 for ( int i = 0; i < numP * numP; i++ )
00401 {
00402 bb_rect.combineExtentWith( x[i], y[i] );
00403 }
00404
00405 QgsDebugMsg( "Projected extent: " + QString(( bb_rect.toString() ).toLocal8Bit().data() ) );
00406
00407 return bb_rect;
00408 }
00409
00410 void QgsCoordinateTransform::transformCoords( const int& numPoints, double *x, double *y, double *z, TransformDirection direction ) const
00411 {
00412
00413 if ( !mSourceCRS.isValid() )
00414 {
00415 QgsLogger::critical( tr( "The source spatial reference system (CRS) is not valid. " )
00416 + tr( "The coordinates can not be reprojected."
00417 " The CRS is: %1" )
00418 .arg( mSourceCRS.toProj4() ) );
00419 return;
00420 }
00421 if ( !mDestCRS.isValid() )
00422 {
00423 QgsLogger::critical( tr( "The destination spatial reference system (CRS) is not valid. " )
00424 + tr( "The coordinates can not be reprojected."
00425 " The CRS is: %1" ).arg( mDestCRS.toProj4() ) );
00426 return;
00427 }
00428
00429 #ifdef COORDINATE_TRANSFORM_VERBOSE
00430 double xorg = *x;
00431 double yorg = *y;
00432 QgsDebugMsg( QString( "[[[[[[ Number of points to transform: %1 ]]]]]]" ).arg( numPoints ) );
00433 #endif
00434
00435
00436 QString dir;
00437
00438
00439 if (( pj_is_latlong( mDestinationProjection ) && ( direction == ReverseTransform ) )
00440 || ( pj_is_latlong( mSourceProjection ) && ( direction == ForwardTransform ) ) )
00441 {
00442 for ( int i = 0; i < numPoints; ++i )
00443 {
00444 x[i] *= DEG_TO_RAD;
00445 y[i] *= DEG_TO_RAD;
00446 z[i] *= DEG_TO_RAD;
00447 }
00448
00449 }
00450 int projResult;
00451 if ( direction == ReverseTransform )
00452 {
00453 projResult = pj_transform( mDestinationProjection, mSourceProjection, numPoints, 0, x, y, z );
00454 dir = tr( "inverse transform" );
00455 }
00456 else
00457 {
00458 Q_ASSERT( mSourceProjection != 0 );
00459 Q_ASSERT( mDestinationProjection != 0 );
00460 projResult = pj_transform( mSourceProjection, mDestinationProjection, numPoints, 0, x, y, z );
00461 dir = tr( "forward transform" );
00462 }
00463
00464 if ( projResult != 0 )
00465 {
00466
00467 QString points;
00468
00469 for ( int i = 0; i < numPoints; ++i )
00470 {
00471 if ( direction == ForwardTransform )
00472 {
00473 points += QString( "(%1, %2)\n" ).arg( x[i] ).arg( y[i] );
00474 }
00475 else
00476 {
00477 points += QString( "(%1, %2)\n" ).arg( x[i] * RAD_TO_DEG ).arg( y[i] * RAD_TO_DEG );
00478 }
00479 }
00480
00481 QString msg = tr( "%1 of\n%2\nfailed with error: %3\n" )
00482 .arg( dir )
00483 .arg( points )
00484 .arg( QString::fromUtf8( pj_strerrno( projResult ) ) );
00485
00486 QgsDebugMsg( "Projection failed emitting invalid transform signal: " + msg );
00487
00488 emit invalidTransformInput();
00489
00490 QgsDebugMsg( "throwing exception" );
00491
00492 throw QgsCsException( msg );
00493 }
00494
00495
00496
00497 if (( pj_is_latlong( mDestinationProjection ) && ( direction == ForwardTransform ) )
00498 || ( pj_is_latlong( mSourceProjection ) && ( direction == ReverseTransform ) ) )
00499 {
00500 for ( int i = 0; i < numPoints; ++i )
00501 {
00502 x[i] *= RAD_TO_DEG;
00503 y[i] *= RAD_TO_DEG;
00504 z[i] *= RAD_TO_DEG;
00505 }
00506 }
00507 #ifdef COORDINATE_TRANSFORM_VERBOSE
00508 QgsDebugMsg( QString( "[[[[[[ Projected %1, %2 to %3, %4 ]]]]]]" )
00509 .arg( xorg, 0, 'g', 15 ).arg( yorg, 0, 'g', 15 )
00510 .arg( *x, 0, 'g', 15 ).arg( *y, 0, 'g', 15 ) );
00511 #endif
00512 }
00513
00514 bool QgsCoordinateTransform::readXML( QDomNode & theNode )
00515 {
00516
00517 QgsDebugMsg( "Reading Coordinate Transform from xml ------------------------!" );
00518
00519 QDomNode mySrcNode = theNode.namedItem( "sourcesrs" );
00520 mSourceCRS.readXML( mySrcNode );
00521
00522 QDomNode myDestNode = theNode.namedItem( "destinationsrs" );
00523 mDestCRS.readXML( myDestNode );
00524
00525 initialise();
00526
00527 return true;
00528 }
00529
00530 bool QgsCoordinateTransform::writeXML( QDomNode & theNode, QDomDocument & theDoc )
00531 {
00532 QDomElement myNodeElement = theNode.toElement();
00533 QDomElement myTransformElement = theDoc.createElement( "coordinatetransform" );
00534
00535 QDomElement mySourceElement = theDoc.createElement( "sourcesrs" );
00536 mSourceCRS.writeXML( mySourceElement, theDoc );
00537 myTransformElement.appendChild( mySourceElement );
00538
00539 QDomElement myDestElement = theDoc.createElement( "destinationsrs" );
00540 mDestCRS.writeXML( myDestElement, theDoc );
00541 myTransformElement.appendChild( myDestElement );
00542
00543 myNodeElement.appendChild( myTransformElement );
00544
00545 return true;
00546 }
00547
00548 const char *finder( const char *name )
00549 {
00550 QString proj;
00551 #ifdef WIN32
00552 proj = QApplication::applicationDirPath()
00553 + "/share/proj/" + QString( name );
00554 #endif
00555 return proj.toUtf8();
00556 }
00557
00558 void QgsCoordinateTransform::setFinder()
00559 {
00560 #ifdef WIN32
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570 pj_set_finder( finder );
00571 #endif
00572 }