00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132 #include <math.h>
00133 #include <stdlib.h>
00134 #include <ctype.h>
00135 #include "DatumLibraryImplementation.h"
00136 #include "EllipsoidLibraryImplementation.h"
00137 #include "EllipsoidParameters.h"
00138 #include "SevenParameterDatum.h"
00139 #include "ThreeParameterDatum.h"
00140 #include "Geocentric.h"
00141 #include "Datum.h"
00142 #include "CartesianCoordinates.h"
00143 #include "GeodeticCoordinates.h"
00144 #include "CoordinateConversionException.h"
00145 #include "ErrorMessages.h"
00146 #include "Accuracy.h"
00147 #include "CCSThreadMutex.h"
00148 #include "CCSThreadLock.h"
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167 using namespace MSP::CCS;
00168 using MSP::CCSThreadMutex;
00169 using MSP::CCSThreadLock;
00170
00171
00172
00173
00174
00175
00176
00177 const double SECONDS_PER_RADIAN = 206264.8062471;
00178 const double PI = 3.14159265358979323e0;
00179 const double PI_OVER_2 = (PI / 2.0);
00180 const double PI_OVER_180 = (PI / 180.0);
00181 const double _180_OVER_PI = (180.0 / PI);
00182 const double TWO_PI = (2.0 * PI);
00183 const double MIN_LAT = (-PI/2.0);
00184 const double MAX_LAT = (+PI/2.0);
00185 const double MIN_LON = -PI;
00186 const double MAX_LON = (2.0 * PI);
00187 const int DATUM_CODE_LENGTH = 7;
00188 const int DATUM_NAME_LENGTH = 33;
00189 const int ELLIPSOID_CODE_LENGTH = 3;
00190 const int MAX_7PARAM = 25;
00191 const int MAX_3PARAM = 250;
00192 const int MAX_WGS = 2;
00193 const double MOLODENSKY_MAX = (89.75 * PI_OVER_180);
00194 const int FILENAME_LENGTH = 128;
00195 const char *WGS84_Datum_Code = "WGE";
00196 const char *WGS72_Datum_Code = "WGC";
00197
00198
00199
00200
00201
00202
00203
00204 GeodeticCoordinates* molodenskyShift( const double a, const double da, const double f, const double df,
00205 const double dx, const double dy, const double dz,
00206 const double sourceLongitude, const double sourceLatitude, const double sourceHeight )
00207
00208 {
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227 double tLon_in;
00228 double e2;
00229 double ep2;
00230 double sin_Lat;
00231 double sin2_Lat;
00232 double sin_Lon;
00233 double cos_Lat;
00234 double cos_Lon;
00235 double w2;
00236 double w;
00237 double w3;
00238 double m;
00239 double n;
00240 double dp;
00241 double dp1;
00242 double dp2;
00243 double dp3;
00244 double dl;
00245 double dh;
00246 double dh1;
00247 double dh2;
00248
00249 if (sourceLongitude > PI)
00250 tLon_in = sourceLongitude - TWO_PI;
00251 else
00252 tLon_in = sourceLongitude;
00253
00254 e2 = 2 * f - f * f;
00255 ep2 = e2 / (1 - e2);
00256 sin_Lat = sin(sourceLatitude);
00257 cos_Lat = cos(sourceLatitude);
00258 sin_Lon = sin(tLon_in);
00259 cos_Lon = cos(tLon_in);
00260 sin2_Lat = sin_Lat * sin_Lat;
00261 w2 = 1.0 - e2 * sin2_Lat;
00262 w = sqrt(w2);
00263 w3 = w * w2;
00264 m = (a * (1.0 - e2)) / w3;
00265 n = a / w;
00266 dp1 = cos_Lat * dz - sin_Lat * cos_Lon * dx - sin_Lat * sin_Lon * dy;
00267 dp2 = ((e2 * sin_Lat * cos_Lat) / w) * da;
00268 dp3 = sin_Lat * cos_Lat * (2.0 * n + ep2 * m * sin2_Lat) * (1.0 - f) * df;
00269 dp = (dp1 + dp2 + dp3) / (m + sourceHeight);
00270 dl = (-sin_Lon * dx + cos_Lon * dy) / ((n + sourceHeight) * cos_Lat);
00271 dh1 = (cos_Lat * cos_Lon * dx) + (cos_Lat * sin_Lon * dy) + (sin_Lat * dz);
00272 dh2 = -(w * da) + ((a * (1 - f)) / w) * sin2_Lat * df;
00273 dh = dh1 + dh2;
00274
00275 double targetLatitude = sourceLatitude + dp;
00276 double targetLongitude = sourceLongitude + dl;
00277 double targetHeight = sourceHeight + dh;
00278
00279 if (targetLongitude > TWO_PI)
00280 targetLongitude -= TWO_PI;
00281 if (targetLongitude < (- PI))
00282 targetLongitude += TWO_PI;
00283
00284 return new GeodeticCoordinates( CoordinateType::geodetic, targetLongitude, targetLatitude, targetHeight );
00285 }
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296 class MSP::CCS::DatumLibraryImplementationCleaner
00297 {
00298 public:
00299
00300 ~DatumLibraryImplementationCleaner()
00301 {
00302 CCSThreadLock lock(&DatumLibraryImplementation::mutex);
00303 DatumLibraryImplementation::deleteInstance();
00304 }
00305
00306 } datumLibraryImplementationCleanerInstance;
00307
00308
00309
00310 CCSThreadMutex DatumLibraryImplementation::mutex;
00311 DatumLibraryImplementation* DatumLibraryImplementation::instance = 0;
00312 int DatumLibraryImplementation::instanceCount = 0;
00313
00314
00315 DatumLibraryImplementation* DatumLibraryImplementation::getInstance()
00316 {
00317 CCSThreadLock lock(&mutex);
00318 if( instance == 0 )
00319 instance = new DatumLibraryImplementation;
00320
00321 instanceCount++;
00322
00323 return instance;
00324 }
00325
00326
00327 void DatumLibraryImplementation::removeInstance()
00328 {
00329
00330
00331
00332
00333 CCSThreadLock lock(&mutex);
00334 if( --instanceCount < 1 )
00335 {
00336 deleteInstance();
00337 }
00338 }
00339
00340
00341 void DatumLibraryImplementation::deleteInstance()
00342 {
00343
00344
00345
00346
00347 if( instance != 0 )
00348 {
00349 delete instance;
00350 instance = 0;
00351 }
00352 }
00353
00354
00355 DatumLibraryImplementation::DatumLibraryImplementation():
00356 _ellipsoidLibraryImplementation( 0 ),
00357 datum3ParamCount( 0 ),
00358 datum7ParamCount( 0 )
00359 {
00360
00361
00362
00363
00364
00365 datumList.reserve( 2 + MAX_3PARAM + MAX_7PARAM );
00366
00367 try
00368 {
00369 loadDatums();
00370 }
00371 catch(CoordinateConversionException e)
00372 {
00373 throw e;
00374 }
00375 }
00376
00377
00378 DatumLibraryImplementation::DatumLibraryImplementation( const DatumLibraryImplementation &dl )
00379 {
00380 int size = dl.datumList.size();
00381 for( int i = 0; i < size; i++ )
00382 {
00383 switch( dl.datumList[i]->datumType() )
00384 {
00385 case DatumType::threeParamDatum:
00386 datumList.push_back( new ThreeParameterDatum( *( ( ThreeParameterDatum* )( dl.datumList[i] ) ) ) );
00387 break;
00388 case DatumType::sevenParamDatum:
00389 datumList.push_back( new SevenParameterDatum( *( ( SevenParameterDatum* )( dl.datumList[i] ) ) ) );
00390 break;
00391 case DatumType::wgs84Datum:
00392 case DatumType::wgs72Datum:
00393 datumList.push_back( new Datum( *( dl.datumList[i] ) ) );
00394 break;
00395 }
00396 }
00397
00398 _ellipsoidLibraryImplementation = dl._ellipsoidLibraryImplementation;
00399 datum3ParamCount = dl.datum3ParamCount;
00400 datum7ParamCount = dl.datum7ParamCount;
00401 }
00402
00403
00404 DatumLibraryImplementation::~DatumLibraryImplementation()
00405 {
00406 std::vector<Datum*>::iterator iter = datumList.begin();
00407 while( iter != datumList.end() )
00408 {
00409 delete( *iter );
00410 iter++;
00411 }
00412 datumList.clear();
00413
00414 _ellipsoidLibraryImplementation = 0;
00415 }
00416
00417
00418 DatumLibraryImplementation& DatumLibraryImplementation::operator=( const DatumLibraryImplementation &dl )
00419 {
00420 if ( &dl == this )
00421 return *this;
00422
00423 int size = dl.datumList.size();
00424 for( int i = 0; i < size; i++ )
00425 {
00426 switch( dl.datumList[i]->datumType() )
00427 {
00428 case DatumType::threeParamDatum:
00429 datumList.push_back( new ThreeParameterDatum( *( ( ThreeParameterDatum* )( dl.datumList[i] ) ) ) );
00430 break;
00431 case DatumType::sevenParamDatum:
00432 datumList.push_back( new SevenParameterDatum( *( ( SevenParameterDatum* )( dl.datumList[i] ) ) ) );
00433 break;
00434 case DatumType::wgs84Datum:
00435 case DatumType::wgs72Datum:
00436 datumList.push_back( new Datum( *( dl.datumList[i] ) ) );
00437 break;
00438 }
00439 }
00440
00441 _ellipsoidLibraryImplementation = dl._ellipsoidLibraryImplementation;
00442 datum3ParamCount = dl.datum3ParamCount;
00443 datum7ParamCount = dl.datum7ParamCount;
00444
00445 return *this;
00446 }
00447
00448
00449 void DatumLibraryImplementation::define3ParamDatum( const char *code, const char *name, const char *ellipsoidCode,
00450 double deltaX, double deltaY, double deltaZ,
00451 double sigmaX, double sigmaY, double sigmaZ,
00452 double westLongitude, double eastLongitude, double southLatitude, double northLatitude )
00453 {
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477 char datum_Code[DATUM_CODE_LENGTH];
00478 long index = 0;
00479 long ellipsoid_index = 0;
00480 long code_length = 0;
00481 char *PathName = NULL;
00482 FILE *fp_3param = NULL;
00483 char errorStatus[256] = "";
00484
00485 if (!(datum3ParamCount < MAX_3PARAM))
00486 strcat( errorStatus, ErrorMessages::datumOverflow );
00487 if (!(((sigmaX > 0.0) || (sigmaX == -1.0)) &&
00488 ((sigmaY > 0.0) || (sigmaY == -1.0)) &&
00489 ((sigmaZ > 0.0) || (sigmaZ == -1.0))))
00490 strcat( errorStatus, ErrorMessages::datumSigma );
00491
00492 if ((southLatitude < MIN_LAT) || (southLatitude > MAX_LAT))
00493 strcat( errorStatus, ErrorMessages::latitude );
00494 if ((westLongitude < MIN_LON) || (westLongitude > MAX_LON))
00495 strcat( errorStatus, ErrorMessages::longitude );
00496 if ((northLatitude < MIN_LAT) || (northLatitude > MAX_LAT))
00497 strcat( errorStatus, ErrorMessages::latitude );
00498 if ((eastLongitude < MIN_LON) || (eastLongitude > MAX_LON))
00499 strcat( errorStatus, ErrorMessages::longitude );
00500 if (southLatitude >= northLatitude)
00501 strcat( errorStatus, ErrorMessages::datumDomain );
00502 if ((westLongitude >= eastLongitude) && (westLongitude >= 0 && westLongitude < 180) && (eastLongitude >= 0 && eastLongitude < 180))
00503 strcat( errorStatus, ErrorMessages::datumDomain );
00504
00505 try
00506 {
00507 datumIndex( code, &index );
00508 strcat( errorStatus, ErrorMessages::invalidDatumCode );
00509 }
00510 catch(CoordinateConversionException e)
00511 {
00512 }
00513
00514 code_length = strlen( code );
00515
00516 if( code_length > ( DATUM_CODE_LENGTH-1 ) )
00517 strcat( errorStatus, ErrorMessages::invalidDatumCode );
00518
00519 if( _ellipsoidLibraryImplementation )
00520 {
00521 _ellipsoidLibraryImplementation->ellipsoidIndex( ellipsoidCode, &ellipsoid_index );
00522
00523
00524 }
00525 else
00526 strcat( errorStatus, ErrorMessages::ellipse );
00527
00528 if( strlen( errorStatus ) > 0)
00529 throw CoordinateConversionException( errorStatus );
00530
00531 strcpy( datum_Code, code );
00532
00533 for( long i = 0; i < code_length; i++ )
00534 datum_Code[i] = ( char )toupper( datum_Code[i] );
00535
00536 int numDatums = datumList.size();
00537 datumList.push_back( new ThreeParameterDatum( numDatums, ( char* )datum_Code, ( char* )ellipsoidCode, ( char* )name, DatumType::threeParamDatum, deltaX, deltaY, deltaZ,
00538 westLongitude, eastLongitude, southLatitude, northLatitude, sigmaX, sigmaY, sigmaZ, true ) );
00539 datum3ParamCount++;
00540
00541 write3ParamFile();
00542 }
00543
00544
00545 void DatumLibraryImplementation::define7ParamDatum( const char *code, const char *name, const char *ellipsoidCode,
00546 double deltaX, double deltaY, double deltaZ,
00547 double rotationX, double rotationY, double rotationZ,
00548 double scale )
00549 {
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570 char datum_Code[DATUM_CODE_LENGTH];
00571 long index = 0;
00572 long ellipsoid_index = 0;
00573 long code_length = 0;
00574 char *PathName = NULL;
00575 FILE *fp_7param = NULL;
00576 char errorStatus[256] = "";
00577
00578 if (!(datum7ParamCount < MAX_7PARAM))
00579 strcat( errorStatus, ErrorMessages::datumOverflow );
00580
00581 if ((rotationX < -60.0) || (rotationX > 60.0))
00582 strcat( errorStatus, ErrorMessages::datumRotation );
00583 if ((rotationY < -60.0) || (rotationY > 60.0))
00584 strcat( errorStatus, ErrorMessages::datumRotation );
00585 if ((rotationZ < -60.0) || (rotationZ > 60.0))
00586 strcat( errorStatus, ErrorMessages::datumRotation );
00587
00588 if ((scale < -0.001) || (scale > 0.001))
00589 strcat( errorStatus, ErrorMessages::scaleFactor );
00590
00591 try
00592 {
00593 datumIndex( code, &index );
00594 strcat( errorStatus, ErrorMessages::invalidDatumCode );
00595 }
00596 catch(CoordinateConversionException e)
00597 {
00598 }
00599
00600 code_length = strlen( code );
00601 if( code_length > ( DATUM_CODE_LENGTH-1 ) )
00602 strcat( errorStatus, ErrorMessages::invalidDatumCode );
00603
00604 if( _ellipsoidLibraryImplementation )
00605 {
00606 _ellipsoidLibraryImplementation->ellipsoidIndex( ellipsoidCode, &ellipsoid_index );
00607
00608
00609 }
00610 else
00611 strcat( errorStatus, ErrorMessages::ellipse );
00612
00613 if( strlen( errorStatus ) > 0)
00614 throw CoordinateConversionException( errorStatus );
00615
00616 long i;
00617 strcpy( datum_Code, code );
00618
00619 for( i = 0; i < code_length; i++ )
00620 datum_Code[i] = ( char )toupper( datum_Code[i] );
00621
00622 datumList.insert( datumList.begin() + MAX_WGS + datum7ParamCount, new SevenParameterDatum( datum7ParamCount, ( char* )datum_Code, ( char* )ellipsoidCode, ( char* )name, DatumType::sevenParamDatum, deltaX, deltaY, deltaZ,
00623 0, 0, 0, 0, rotationX / SECONDS_PER_RADIAN, rotationY / SECONDS_PER_RADIAN, rotationZ / SECONDS_PER_RADIAN, scale, true ) );
00624 datum7ParamCount++;
00625
00626 write7ParamFile();
00627 }
00628
00629
00630 void DatumLibraryImplementation::removeDatum( const char* code )
00631 {
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643 char *PathName = NULL;
00644 FILE *fp_3param = NULL;
00645 FILE *fp_7param = NULL;
00646 long index = 0;
00647 bool delete_3param_datum = true;
00648
00649 datumIndex( code, &index );
00650
00651 if( datumList[index]->datumType() == DatumType::threeParamDatum )
00652 {
00653 if( !( ( ThreeParameterDatum* )datumList[index] )->userDefined() )
00654 throw CoordinateConversionException( ErrorMessages::notUserDefined );
00655 }
00656 else if( datumList[index]->datumType() == DatumType::sevenParamDatum )
00657 {
00658 delete_3param_datum = false;
00659 if( !( ( SevenParameterDatum* )datumList[index] )->userDefined() )
00660 throw CoordinateConversionException( ErrorMessages::notUserDefined );
00661 }
00662 else
00663 throw CoordinateConversionException( ErrorMessages::notUserDefined );
00664
00665 long i = 0;
00666 long count = 0;
00667
00668 int numDatums = datumList.size();
00669 for(i = index; i < numDatums; i++)
00670 datumList[i] = datumList[i+1];
00671
00672 if( numDatums != ( 2 + MAX_3PARAM + MAX_7PARAM ) )
00673 datumList[i] = datumList[i+1];
00674 else
00675 datumList.erase( datumList.end() - 1 );
00676
00677 numDatums--;
00678 datumList.resize( numDatums );
00679
00680 if( !delete_3param_datum )
00681 {
00682 datum7ParamCount--;
00683
00684 write7ParamFile();
00685 }
00686 else if( delete_3param_datum )
00687 {
00688 datum3ParamCount--;
00689
00690 write3ParamFile();
00691 }
00692 }
00693
00694
00695 void DatumLibraryImplementation::datumCount( long *count )
00696 {
00697
00698
00699
00700
00701
00702
00703
00704 *count = datumList.size();
00705 }
00706
00707
00708 void DatumLibraryImplementation::datumIndex( const char *code, long *index )
00709 {
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719 char temp_code[DATUM_CODE_LENGTH];
00720 long length;
00721 long pos = 0;
00722 long i = 0;
00723
00724 *index = 0;
00725
00726 if( !code )
00727 throw CoordinateConversionException( ErrorMessages::invalidDatumCode );
00728
00729 length = strlen( code );
00730 if ( length > ( DATUM_CODE_LENGTH-1 ) )
00731 throw CoordinateConversionException( ErrorMessages::invalidDatumCode );
00732 else
00733 {
00734 strcpy( temp_code, code );
00735
00736
00737 for( i=0; i < length; i++ )
00738 temp_code[i] = ( char )toupper( temp_code[i] );
00739
00740
00741 while( pos < length )
00742 {
00743 if( isspace( temp_code[pos] ) )
00744 {
00745 for( i=pos; i <= length; i++ )
00746 temp_code[i] = temp_code[i+1];
00747 length -= 1;
00748 }
00749 else
00750 pos += 1;
00751 }
00752
00753 int numDatums = datumList.size();
00754
00755 i = 0;
00756 while( i < numDatums && strcmp( temp_code, datumList[i]->code() ) )
00757 {
00758 i++;
00759 }
00760 if( i == numDatums || strcmp( temp_code, datumList[i]->code() ) )
00761 throw CoordinateConversionException( ErrorMessages::invalidDatumCode );
00762 else
00763 *index = i;
00764 }
00765 }
00766
00767
00768 void DatumLibraryImplementation::datumCode( const long index, char *code )
00769 {
00770
00771
00772
00773
00774
00775
00776
00777
00778 if( index < 0 || index >= datumList.size() )
00779 throw CoordinateConversionException( ErrorMessages::invalidIndex );
00780 else
00781 strcpy( code, datumList[index]->code() );
00782 }
00783
00784
00785 void DatumLibraryImplementation::datumName( const long index, char *name )
00786 {
00787
00788
00789
00790
00791
00792
00793
00794
00795 if( index < 0 || index >= datumList.size() )
00796 throw CoordinateConversionException( ErrorMessages::invalidIndex );
00797 else
00798 strcpy( name, datumList[index]->name() );
00799 }
00800
00801
00802 void DatumLibraryImplementation::datumEllipsoidCode( const long index, char *code )
00803 {
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813 if( index < 0 || index >= datumList.size() )
00814 throw CoordinateConversionException( ErrorMessages::invalidIndex );
00815 else
00816 strcpy( code, datumList[index]->ellipsoidCode() );
00817 }
00818
00819
00820 void DatumLibraryImplementation::datumStandardErrors( const long index, double *sigmaX, double *sigmaY, double *sigmaZ )
00821 {
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832 if( index < 0 || index >= datumList.size() )
00833 throw CoordinateConversionException( ErrorMessages::invalidIndex );
00834 else
00835 {
00836 Datum* datum = datumList[index];
00837
00838 if( datum->datumType() == DatumType::threeParamDatum )
00839 {
00840 ThreeParameterDatum* threeParameterDatum = ( ThreeParameterDatum* )datum;
00841 *sigmaX = threeParameterDatum->sigmaX();
00842 *sigmaY = threeParameterDatum->sigmaY();
00843 *sigmaZ = threeParameterDatum->sigmaZ();
00844 }
00845 else
00846 throw CoordinateConversionException( ErrorMessages::invalidIndex );
00847 }
00848 }
00849
00850
00851 void DatumLibraryImplementation::datumSevenParameters( const long index, double *rotationX, double *rotationY, double *rotationZ, double *scaleFactor )
00852 {
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865 if( index < 0 || index >= datumList.size() )
00866 throw CoordinateConversionException( ErrorMessages::invalidIndex );
00867 else
00868 {
00869 Datum* datum = datumList[index];
00870
00871 if( datum->datumType() == DatumType::sevenParamDatum )
00872 {
00873 SevenParameterDatum* sevenParameterDatum = ( SevenParameterDatum* )datum;
00874
00875 *rotationX = sevenParameterDatum->rotationX();
00876 *rotationY = sevenParameterDatum->rotationY();
00877 *rotationZ = sevenParameterDatum->rotationZ();
00878 *scaleFactor = sevenParameterDatum->scaleFactor();
00879 }
00880 else
00881 throw CoordinateConversionException( ErrorMessages::invalidIndex );
00882 }
00883 }
00884
00885
00886 void DatumLibraryImplementation::datumTranslationValues( const long index, double *deltaX, double *deltaY, double *deltaZ )
00887 {
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898 if( index < 0 || index >= datumList.size() )
00899 throw CoordinateConversionException( ErrorMessages::invalidIndex );
00900 else
00901 {
00902 Datum* datum = datumList[index];
00903
00904 *deltaX = datum->deltaX();
00905 *deltaY = datum->deltaY();
00906 *deltaZ = datum->deltaZ();
00907 }
00908 }
00909
00910
00911 Accuracy* DatumLibraryImplementation::datumShiftError( const long sourceIndex, const long targetIndex,
00912 double longitude, double latitude, Accuracy* sourceAccuracy )
00913 {
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928 double sinlat = sin( latitude );
00929 double coslat = cos( latitude );
00930 double sinlon = sin( longitude );
00931 double coslon = cos( longitude );
00932 double sigma_delta_lat;
00933 double sigma_delta_lon;
00934 double sigma_delta_height;
00935 double sx, sy, sz;
00936 double ce90_in = -1.0;
00937 double le90_in = -1.0;
00938 double se90_in = -1.0;
00939 double ce90_out = -1.0;
00940 double le90_out = -1.0;
00941 double se90_out = -1.0;
00942 double circularError90 = sourceAccuracy->circularError90();
00943 double linearError90 = sourceAccuracy->linearError90();
00944 double sphericalError90 = sourceAccuracy->sphericalError90();
00945
00946 int numDatums = datumList.size();
00947
00948 if( ( sourceIndex < 0 ) || ( sourceIndex >= numDatums ) )
00949 throw CoordinateConversionException( ErrorMessages::invalidIndex );
00950 if( ( targetIndex < 0 ) || ( targetIndex >= numDatums ) )
00951 throw CoordinateConversionException( ErrorMessages::invalidIndex );
00952 if( ( latitude < ( -90 * PI_OVER_180 ) ) || ( latitude > ( 90 * PI_OVER_180 ) ) )
00953 throw CoordinateConversionException( ErrorMessages::latitude );
00954 if( ( longitude < ( -PI ) ) || ( longitude > TWO_PI ) )
00955 throw CoordinateConversionException( ErrorMessages::longitude );
00956
00957 if( sourceIndex == targetIndex )
00958 {
00959 circularError90 = circularError90;
00960 linearError90 = linearError90;
00961 sphericalError90 = sphericalError90;
00962 }
00963 else
00964 {
00965 Datum* sourceDatum = datumList[sourceIndex];
00966 Datum* targetDatum = datumList[targetIndex];
00967
00968
00969 switch( sourceDatum->datumType() )
00970 {
00971 case DatumType::wgs84Datum:
00972 case DatumType::wgs72Datum:
00973 case DatumType::sevenParamDatum:
00974 {
00975 ce90_in = 0.0;
00976 le90_in = 0.0;
00977 se90_in = 0.0;
00978 break;
00979 }
00980 case DatumType::threeParamDatum:
00981 {
00982 ThreeParameterDatum* sourceThreeParameterDatum = ( ThreeParameterDatum* )sourceDatum;
00983 if( ( sourceThreeParameterDatum->sigmaX() < 0 )
00984 || ( sourceThreeParameterDatum->sigmaY() < 0 )
00985 || ( sourceThreeParameterDatum->sigmaZ() < 0 ) )
00986 {
00987 ce90_in = -1.0;
00988 le90_in = -1.0;
00989 se90_in = -1.0;
00990 }
00991 else
00992 {
00993 sx = sourceThreeParameterDatum->sigmaX() * sinlat * coslon;
00994 sy = sourceThreeParameterDatum->sigmaY() * sinlat * sinlon;
00995 sz = sourceThreeParameterDatum->sigmaZ() * coslat;
00996 sigma_delta_lat = sqrt( ( sx * sx ) + ( sy * sy ) + ( sz * sz ) );
00997
00998 sx = sourceThreeParameterDatum->sigmaX() * sinlon;
00999 sy = sourceThreeParameterDatum->sigmaY() * coslon;
01000 sigma_delta_lon = sqrt( ( sx * sx ) + ( sy * sy ) );
01001
01002 sx = sourceThreeParameterDatum->sigmaX() * coslat * coslon;
01003 sy = sourceThreeParameterDatum->sigmaY() * coslat * sinlon;
01004 sz = sourceThreeParameterDatum->sigmaZ() * sinlat;
01005 sigma_delta_height = sqrt( ( sx * sx ) + ( sy * sy ) + ( sz * sz ) );
01006
01007 ce90_in = 2.146 * ( sigma_delta_lat + sigma_delta_lon ) / 2.0;
01008 le90_in = 1.6449 * sigma_delta_height;
01009 se90_in = 2.5003 * ( sourceThreeParameterDatum->sigmaX() + sourceThreeParameterDatum->sigmaY() + sourceThreeParameterDatum->sigmaZ() ) / 3.0;
01010 }
01011 break;
01012 }
01013 }
01014
01015
01016 switch( targetDatum->datumType() )
01017 {
01018 case DatumType::wgs84Datum:
01019 case DatumType::wgs72Datum:
01020 case DatumType::sevenParamDatum:
01021 {
01022 ce90_out = 0.0;
01023 le90_out = 0.0;
01024 se90_out = 0.0;
01025 break;
01026 }
01027 case DatumType::threeParamDatum:
01028 {
01029 ThreeParameterDatum* targetThreeParameterDatum = ( ThreeParameterDatum* )targetDatum;
01030 if ((targetThreeParameterDatum->sigmaX() < 0)
01031 ||(targetThreeParameterDatum->sigmaY() < 0)
01032 ||(targetThreeParameterDatum->sigmaZ() < 0))
01033 {
01034 ce90_out = -1.0;
01035 le90_out = -1.0;
01036 se90_out = -1.0;
01037 }
01038 else
01039 {
01040 sx = targetThreeParameterDatum->sigmaX() * sinlat * coslon;
01041 sy = targetThreeParameterDatum->sigmaY() * sinlat * sinlon;
01042 sz = targetThreeParameterDatum->sigmaZ() * coslat;
01043 sigma_delta_lat = sqrt((sx * sx) + (sy * sy) + (sz * sz));
01044
01045 sx = targetThreeParameterDatum->sigmaX() * sinlon;
01046 sy = targetThreeParameterDatum->sigmaY() * coslon;
01047 sigma_delta_lon = sqrt( ( sx * sx ) + ( sy * sy ) );
01048
01049 sx = targetThreeParameterDatum->sigmaX() * coslat * coslon;
01050 sy = targetThreeParameterDatum->sigmaY() * coslat * sinlon;
01051 sz = targetThreeParameterDatum->sigmaZ() * sinlat;
01052 sigma_delta_height = sqrt( ( sx * sx ) + ( sy * sy ) + ( sz * sz ) );
01053
01054 ce90_out = 2.146 * ( sigma_delta_lat + sigma_delta_lon ) / 2.0;
01055 le90_out = 1.6449 * sigma_delta_height;
01056 se90_out = 2.5003 * ( targetThreeParameterDatum->sigmaX() + targetThreeParameterDatum->sigmaY() + targetThreeParameterDatum->sigmaZ() ) / 3.0;
01057 }
01058 break;
01059 }
01060 }
01061
01062
01063 if( ( circularError90 < 0.0 ) || ( ce90_in < 0.0 ) || ( ce90_out < 0.0 ) )
01064 {
01065 circularError90 = -1.0;
01066 linearError90 = -1.0;
01067 sphericalError90 = -1.0;
01068 }
01069 else
01070 {
01071 circularError90 = sqrt( ( circularError90 * circularError90 ) + ( ce90_in * ce90_in ) + ( ce90_out * ce90_out ) );
01072 if( circularError90 < 1.0 )
01073 {
01074 circularError90 = 1.0;
01075 }
01076
01077 if( ( linearError90 < 0.0 ) || ( le90_in < 0.0 ) || ( le90_out < 0.0 ) )
01078 {
01079 linearError90 = -1.0;
01080 sphericalError90 = -1.0;
01081 }
01082 else
01083 {
01084 linearError90 = sqrt( ( linearError90 * linearError90 ) + ( le90_in * le90_in ) + ( le90_out * le90_out ) );
01085 if( linearError90 < 1.0 )
01086 {
01087 linearError90 = 1.0;
01088 }
01089
01090 if( ( sphericalError90 < 0.0 ) || ( se90_in < 0.0 ) || ( se90_out < 0.0 ) )
01091 sphericalError90 = -1.0;
01092 else
01093 {
01094 sphericalError90 = sqrt( ( sphericalError90 * sphericalError90 ) + ( se90_in * se90_in ) + ( se90_out * se90_out ) );
01095 if( sphericalError90 < 1.0 )
01096 {
01097 sphericalError90 = 1.0;
01098 }
01099 }
01100 }
01101 }
01102 }
01103
01104 return new Accuracy(circularError90, linearError90, sphericalError90);
01105 }
01106
01107
01108 void DatumLibraryImplementation::datumUserDefined( const long index, long *result )
01109 {
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121 *result = false;
01122
01123 if( index < 0 || index >= datumList.size() )
01124 throw CoordinateConversionException( ErrorMessages::invalidIndex );
01125 else
01126 {
01127 Datum* datum = datumList[index];
01128
01129 if( datum->datumType() == DatumType::threeParamDatum )
01130 {
01131 ThreeParameterDatum* threeParameterDatum = ( ThreeParameterDatum* )datum;
01132
01133 if( threeParameterDatum->userDefined() )
01134 *result = true;
01135 }
01136 else if( datum->datumType() == DatumType::sevenParamDatum )
01137 {
01138 SevenParameterDatum* sevenParamDatum = ( SevenParameterDatum* )datum;
01139
01140 if( sevenParamDatum->userDefined() )
01141 *result = true;
01142 }
01143 else
01144 throw CoordinateConversionException( ErrorMessages::invalidIndex );
01145 }
01146 }
01147
01148
01149 bool DatumLibraryImplementation::datumUsesEllipsoid( const char *ellipsoidCode )
01150 {
01151
01152
01153
01154
01155
01156
01157
01158 char temp_code[DATUM_CODE_LENGTH];
01159 long length;
01160 long pos = 0;
01161 long i = 0;
01162 bool ellipsoid_in_use = false;
01163
01164 length = strlen( ellipsoidCode );
01165 if( length <= ( ELLIPSOID_CODE_LENGTH-1 ) )
01166 {
01167 strcpy( temp_code, ellipsoidCode );
01168
01169
01170 for( i=0;i<length;i++ )
01171 temp_code[i] = ( char )toupper( temp_code[i] );
01172
01173
01174 while( pos < length )
01175 {
01176 if( isspace( temp_code[pos] ) )
01177 {
01178 for( i=pos; i<=length; i++ )
01179 temp_code[i] = temp_code[i+1];
01180 length -= 1;
01181 }
01182 else
01183 pos += 1;
01184 }
01185
01186
01187 i = 0;
01188 int numDatums = datumList.size();
01189 while( ( i < numDatums ) && ( !ellipsoid_in_use ) )
01190 {
01191 if( strcmp( temp_code, datumList[i]->ellipsoidCode() ) == 0 )
01192 ellipsoid_in_use = true;
01193 i++;
01194 }
01195 }
01196
01197 return ellipsoid_in_use;
01198 }
01199
01200
01201 void DatumLibraryImplementation::datumValidRectangle( const long index, double *westLongitude, double *eastLongitude, double *southLatitude, double *northLatitude )
01202 {
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215 if( index < 0 && index >= datumList.size() )
01216 throw CoordinateConversionException( ErrorMessages::invalidIndex );
01217 else
01218 {
01219 *westLongitude = datumList[index]->westLongitude();
01220 *eastLongitude = datumList[index]->eastLongitude();
01221 *southLatitude = datumList[index]->southLatitude();
01222 *northLatitude = datumList[index]->northLatitude();
01223 }
01224 }
01225
01226
01227 CartesianCoordinates* DatumLibraryImplementation::geocentricDatumShift( const long sourceIndex, const double sourceX, const double sourceY, const double sourceZ,
01228 const long targetIndex )
01229 {
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246 char errorStatus[50] = "";
01247
01248 int numDatums = datumList.size();
01249
01250 if( ( sourceIndex < 0 ) || ( sourceIndex >= numDatums ) )
01251 strcat( errorStatus, ErrorMessages::invalidIndex );
01252 if( ( targetIndex < 0 ) || ( targetIndex >= numDatums ) )
01253 strcat( errorStatus, ErrorMessages::invalidIndex );
01254
01255 if( strlen( errorStatus ) > 0)
01256 throw CoordinateConversionException( errorStatus );
01257
01258 if( sourceIndex == targetIndex )
01259 {
01260 return new CartesianCoordinates( CoordinateType::geocentric, sourceX, sourceY, sourceZ );
01261 }
01262 else
01263 {
01264 CartesianCoordinates* wgs84CartesianCoordinates = geocentricShiftToWGS84( sourceIndex, sourceX, sourceY, sourceZ );
01265
01266 CartesianCoordinates* targetCartesianCoordinates = geocentricShiftFromWGS84( wgs84CartesianCoordinates->x(), wgs84CartesianCoordinates->y(), wgs84CartesianCoordinates->z(), targetIndex );
01267
01268 delete wgs84CartesianCoordinates;
01269
01270 return targetCartesianCoordinates;
01271 }
01272 }
01273
01274
01275 CartesianCoordinates* DatumLibraryImplementation::geocentricShiftFromWGS84( const double WGS84X, const double WGS84Y, const double WGS84Z, const long targetIndex )
01276 {
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291 char errorStatus[50] = "";
01292
01293 int numDatums = datumList.size();
01294
01295 if( ( targetIndex < 0 ) || ( targetIndex >= numDatums ) )
01296 strcat( errorStatus, ErrorMessages::invalidIndex );
01297
01298 if( strlen( errorStatus ) > 0)
01299 throw CoordinateConversionException( errorStatus );
01300
01301 Datum* localDatum = datumList[targetIndex];
01302 switch( localDatum->datumType() )
01303 {
01304 case DatumType::wgs72Datum:
01305 {
01306 CartesianCoordinates* wgs72CartesianCoordinates = geocentricShiftWGS84ToWGS72( WGS84X, WGS84Y, WGS84Z );
01307
01308 return wgs72CartesianCoordinates;
01309 }
01310 case DatumType::wgs84Datum:
01311 {
01312 return new CartesianCoordinates( CoordinateType::geocentric, WGS84X, WGS84Y, WGS84Z );
01313 }
01314 case DatumType::sevenParamDatum:
01315 {
01316 SevenParameterDatum* sevenParameterDatum = ( SevenParameterDatum* )localDatum;
01317
01318 double targetX = WGS84X - sevenParameterDatum->deltaX() - sevenParameterDatum->rotationZ() * WGS84Y
01319 + sevenParameterDatum->rotationY() * WGS84Z - sevenParameterDatum->scaleFactor() * WGS84X;
01320
01321 double targetY = WGS84Y - sevenParameterDatum->deltaY() + sevenParameterDatum->rotationZ() * WGS84X
01322 - sevenParameterDatum->rotationX() * WGS84Z - sevenParameterDatum->scaleFactor() * WGS84Y;
01323
01324 double targetZ = WGS84Z - sevenParameterDatum->deltaZ() - sevenParameterDatum->rotationY() * WGS84X
01325 + sevenParameterDatum->rotationX() * WGS84Y - sevenParameterDatum->scaleFactor() * WGS84Z;
01326
01327 return new CartesianCoordinates( CoordinateType::geocentric, targetX, targetY, targetZ );
01328 }
01329 case DatumType::threeParamDatum:
01330 {
01331 ThreeParameterDatum* threeParameterDatum = ( ThreeParameterDatum* )localDatum;
01332
01333 double targetX = WGS84X - threeParameterDatum->deltaX();
01334 double targetY = WGS84Y - threeParameterDatum->deltaY();
01335 double targetZ = WGS84Z - threeParameterDatum->deltaZ();
01336
01337 return new CartesianCoordinates( CoordinateType::geocentric, targetX, targetY, targetZ );
01338 }
01339 default:
01340 throw CoordinateConversionException( ErrorMessages::datumType );
01341 }
01342 }
01343
01344
01345 CartesianCoordinates* DatumLibraryImplementation::geocentricShiftToWGS84( const long sourceIndex, const double sourceX, const double sourceY, const double sourceZ )
01346 {
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361 char errorStatus[50] = "";
01362
01363 int numDatums = datumList.size();
01364
01365 if( ( sourceIndex < 0 ) || (sourceIndex > numDatums ) )
01366 strcat( errorStatus, ErrorMessages::invalidIndex );
01367
01368 if( strlen( errorStatus ) > 0)
01369 throw CoordinateConversionException( errorStatus );
01370
01371 Datum* localDatum = datumList[sourceIndex];
01372 switch( localDatum->datumType() )
01373 {
01374 case DatumType::wgs72Datum:
01375 {
01376 CartesianCoordinates* wgs84CartesianCoordinates = geocentricShiftWGS72ToWGS84( sourceX, sourceY, sourceZ );
01377
01378 return wgs84CartesianCoordinates;
01379 }
01380 case DatumType::wgs84Datum:
01381 {
01382 return new CartesianCoordinates( CoordinateType::geocentric, sourceX, sourceY, sourceZ );
01383 }
01384 case DatumType::sevenParamDatum:
01385 {
01386 SevenParameterDatum* sevenParameterDatum = ( SevenParameterDatum* )localDatum;
01387
01388 double wgs84X = sourceX + sevenParameterDatum->deltaX() + sevenParameterDatum->rotationZ() * sourceY
01389 - sevenParameterDatum->rotationY() * sourceZ + sevenParameterDatum->scaleFactor() * sourceX;
01390
01391 double wgs84Y = sourceY + sevenParameterDatum->deltaY() - sevenParameterDatum->rotationZ() * sourceX
01392 + sevenParameterDatum->rotationX() * sourceZ + sevenParameterDatum->scaleFactor() * sourceY;
01393
01394 double wgs84Z = sourceZ + sevenParameterDatum->deltaZ() + sevenParameterDatum->rotationY() * sourceX
01395 - sevenParameterDatum->rotationX() * sourceY + sevenParameterDatum->scaleFactor() * sourceZ;
01396
01397 return new CartesianCoordinates( CoordinateType::geocentric, wgs84X, wgs84Y, wgs84Z );
01398 }
01399 case DatumType::threeParamDatum:
01400 {
01401 ThreeParameterDatum* threeParameterDatum = ( ThreeParameterDatum* )localDatum;
01402
01403 double wgs84X = sourceX + threeParameterDatum->deltaX();
01404 double wgs84Y = sourceY + threeParameterDatum->deltaY();
01405 double wgs84Z = sourceZ + threeParameterDatum->deltaZ();
01406
01407 return new CartesianCoordinates( CoordinateType::geocentric, wgs84X, wgs84Y, wgs84Z );
01408 }
01409 default:
01410 throw CoordinateConversionException( ErrorMessages::datumType );
01411 }
01412 }
01413
01414
01415 GeodeticCoordinates* DatumLibraryImplementation::geodeticDatumShift( const long sourceIndex, const GeodeticCoordinates* sourceCoordinates,
01416 const long targetIndex )
01417 {
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434 long E_Index;
01435 double a;
01436 double f;
01437 char errorStatus[50] = "";
01438
01439 double sourceLongitude = sourceCoordinates->longitude();
01440 double sourceLatitude = sourceCoordinates->latitude();
01441 double sourceHeight = sourceCoordinates->height();
01442
01443 int numDatums = datumList.size();
01444
01445 if( sourceIndex < 0 || sourceIndex >= numDatums )
01446 strcat( errorStatus, ErrorMessages::invalidIndex );
01447 if( targetIndex < 0 || targetIndex >= numDatums )
01448 strcat( errorStatus, ErrorMessages::invalidIndex );
01449
01450 if( ( sourceLatitude < ( -90 * PI_OVER_180 ) ) || ( sourceLatitude > ( 90 * PI_OVER_180 ) ) )
01451 strcat( errorStatus, ErrorMessages::latitude );
01452 if( ( sourceLongitude < ( -PI )) || ( sourceLongitude > TWO_PI ) )
01453 strcat( errorStatus, ErrorMessages::longitude );
01454
01455 if( strlen( errorStatus ) > 0)
01456 throw CoordinateConversionException( errorStatus );
01457
01458 Datum* sourceDatum = datumList[sourceIndex];
01459 Datum* targetDatum = datumList[targetIndex];
01460
01461 if ( sourceIndex == targetIndex )
01462 {
01463 return new GeodeticCoordinates(CoordinateType::geodetic, sourceLatitude, sourceLongitude, sourceHeight);
01464 }
01465 else
01466 {
01467 if( _ellipsoidLibraryImplementation )
01468 {
01469 if( sourceDatum->datumType() == DatumType::sevenParamDatum )
01470 {
01471 _ellipsoidLibraryImplementation->ellipsoidIndex( sourceDatum->ellipsoidCode(), &E_Index );
01472 _ellipsoidLibraryImplementation->ellipsoidParameters( E_Index, &a, &f );
01473
01474
01475
01476
01477
01478
01479 Geocentric geocentricFromGeodetic( a, f );
01480
01481 CartesianCoordinates* sourceCartesianCoordinates = geocentricFromGeodetic.convertFromGeodetic( sourceCoordinates );
01482
01483 if( targetDatum->datumType() == DatumType::sevenParamDatum )
01484 {
01485 CartesianCoordinates* targetCartesianCoordinates = geocentricDatumShift( sourceIndex, sourceCartesianCoordinates->x(), sourceCartesianCoordinates->y(), sourceCartesianCoordinates->z(), targetIndex );
01486
01487 _ellipsoidLibraryImplementation->ellipsoidIndex( targetDatum->ellipsoidCode(), &E_Index );
01488 _ellipsoidLibraryImplementation->ellipsoidParameters( E_Index, &a, &f );
01489
01490
01491
01492
01493
01494
01495 Geocentric geocentricToGeodetic( a, f );
01496
01497 GeodeticCoordinates* targetGeodeticCoordinates = geocentricToGeodetic.convertToGeodetic( targetCartesianCoordinates );
01498
01499 delete targetCartesianCoordinates;
01500 delete sourceCartesianCoordinates;
01501
01502 return targetGeodeticCoordinates;
01503 }
01504 else
01505 {
01506 CartesianCoordinates* wgs84CartesianCoordinates = geocentricShiftToWGS84( sourceIndex, sourceCartesianCoordinates->x(), sourceCartesianCoordinates->y(), sourceCartesianCoordinates->z() );
01507
01508 long wgs84EllipsoidIndex;
01509 _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
01510 _ellipsoidLibraryImplementation->ellipsoidParameters( wgs84EllipsoidIndex, &a, &f );
01511
01512 Geocentric geocentricToGeodetic( a, f );
01513
01514 GeodeticCoordinates* wgs84GeodeticCoordinates = geocentricToGeodetic.convertToGeodetic( wgs84CartesianCoordinates );
01515
01516 GeodeticCoordinates* targetGeodeticCoordinates = geodeticShiftFromWGS84( wgs84GeodeticCoordinates, targetIndex );
01517
01518 delete wgs84GeodeticCoordinates;
01519 delete wgs84CartesianCoordinates;
01520
01521 return targetGeodeticCoordinates;
01522 }
01523 }
01524 else if( targetDatum->datumType() == DatumType::sevenParamDatum )
01525 {
01526 GeodeticCoordinates* wgs84GeodeticCoordinates = geodeticShiftToWGS84( sourceIndex, sourceCoordinates );
01527
01528 long wgs84EllipsoidIndex;
01529 _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
01530 _ellipsoidLibraryImplementation->ellipsoidParameters( wgs84EllipsoidIndex, &a, &f );
01531
01532 Geocentric geocentricFromGeodetic( a, f );
01533
01534 CartesianCoordinates* wgs84CartesianCoordinates = geocentricFromGeodetic.convertFromGeodetic( wgs84GeodeticCoordinates );
01535
01536 CartesianCoordinates* targetCartesianCoordinates = geocentricShiftFromWGS84( wgs84CartesianCoordinates->x(), wgs84CartesianCoordinates->y(), wgs84CartesianCoordinates->z(), targetIndex );
01537
01538 _ellipsoidLibraryImplementation->ellipsoidIndex( targetDatum->ellipsoidCode(), &E_Index );
01539 _ellipsoidLibraryImplementation->ellipsoidParameters( E_Index, &a, &f );
01540
01541
01542
01543
01544
01545 Geocentric geocentricToGeodetic( a, f );
01546
01547 GeodeticCoordinates* targetGeodeticCoordinates = geocentricToGeodetic.convertToGeodetic( targetCartesianCoordinates );
01548
01549 delete targetCartesianCoordinates;
01550 delete wgs84CartesianCoordinates;
01551 delete wgs84GeodeticCoordinates;
01552
01553 return targetGeodeticCoordinates;
01554 }
01555 else
01556 {
01557 GeodeticCoordinates* wgs84GeodeticCoordinates = geodeticShiftToWGS84( sourceIndex, sourceCoordinates );
01558
01559 GeodeticCoordinates* targetGeodeticCoordinates = geodeticShiftFromWGS84( wgs84GeodeticCoordinates, targetIndex );
01560
01561 delete wgs84GeodeticCoordinates;
01562
01563 return targetGeodeticCoordinates;
01564 }
01565 }
01566 else
01567 throw CoordinateConversionException( ErrorMessages::ellipse );
01568 }
01569 }
01570
01571
01572 GeodeticCoordinates* DatumLibraryImplementation::geodeticShiftFromWGS84( const GeodeticCoordinates* sourceCoordinates, const long targetIndex )
01573 {
01574
01575
01576
01577
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587
01588 double WGS84_a;
01589 double WGS84_f;
01590 double a;
01591 double da;
01592 double f;
01593 double df;
01594 double dx;
01595 double dy;
01596 double dz;
01597 long E_Index;
01598 char errorStatus[50] = "";
01599
01600 double WGS84Longitude = sourceCoordinates->longitude();
01601 double WGS84Latitude = sourceCoordinates->latitude();
01602 double WGS84Height = sourceCoordinates->height();
01603
01604 if( ( targetIndex < 0) || (targetIndex >= datumList.size() ) )
01605 strcat( errorStatus, ErrorMessages::invalidIndex );
01606 if( ( WGS84Latitude < ( -90 * PI_OVER_180 ) ) || ( WGS84Latitude > ( 90 * PI_OVER_180 ) ) )
01607 strcat( errorStatus, ErrorMessages::latitude );
01608 if( ( WGS84Longitude < ( -PI ) ) || ( WGS84Longitude > TWO_PI ) )
01609 strcat( errorStatus, ErrorMessages::longitude );
01610
01611 if( strlen( errorStatus ) > 0)
01612 throw CoordinateConversionException( errorStatus );
01613
01614 Datum* localDatum = datumList[targetIndex];
01615 switch( localDatum->datumType() )
01616 {
01617 case DatumType::wgs72Datum:
01618 {
01619 GeodeticCoordinates* targetGeodeticCoordinates = geodeticShiftWGS84ToWGS72( WGS84Longitude, WGS84Latitude, WGS84Height );
01620 return targetGeodeticCoordinates;
01621 }
01622 case DatumType::wgs84Datum:
01623 {
01624 return new GeodeticCoordinates( CoordinateType::geodetic, WGS84Longitude, WGS84Latitude, WGS84Height );
01625 }
01626 case DatumType::sevenParamDatum:
01627 case DatumType::threeParamDatum:
01628 {
01629 if( _ellipsoidLibraryImplementation )
01630 {
01631 _ellipsoidLibraryImplementation->ellipsoidIndex( localDatum->ellipsoidCode(), &E_Index );
01632 _ellipsoidLibraryImplementation->ellipsoidParameters( E_Index, &a, &f );
01633
01634
01635
01636
01637
01638 long wgs84EllipsoidIndex;
01639 _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
01640 _ellipsoidLibraryImplementation->ellipsoidParameters( wgs84EllipsoidIndex, &WGS84_a, &WGS84_f );
01641
01642 if( ( localDatum->datumType() == DatumType::sevenParamDatum ) ||
01643 ( WGS84Latitude < ( -MOLODENSKY_MAX ) ) ||
01644 ( WGS84Latitude > MOLODENSKY_MAX ) )
01645 {
01646 Geocentric geocentricFromGeodetic( WGS84_a, WGS84_f );
01647
01648 CartesianCoordinates* wgs84CartesianCoordinates = geocentricFromGeodetic.convertFromGeodetic( sourceCoordinates );
01649
01650 CartesianCoordinates* localCartesianCoordinates = geocentricShiftFromWGS84( wgs84CartesianCoordinates->x(), wgs84CartesianCoordinates->y(),
01651 wgs84CartesianCoordinates->z(), targetIndex );
01652
01653 Geocentric geocentricToGeodetic( a, f );
01654
01655 GeodeticCoordinates* targetGeodeticCoordinates = geocentricToGeodetic.convertToGeodetic( localCartesianCoordinates );
01656
01657 delete localCartesianCoordinates;
01658 delete wgs84CartesianCoordinates;
01659
01660 return targetGeodeticCoordinates;
01661 }
01662 else
01663 {
01664 da = a - WGS84_a;
01665 df = f - WGS84_f;
01666 dx = -( localDatum->deltaX() );
01667 dy = -( localDatum->deltaY() );
01668 dz = -( localDatum->deltaZ() );
01669
01670 GeodeticCoordinates* targetGeodeticCoordinates = molodenskyShift( WGS84_a, da, WGS84_f, df, dx, dy, dz,
01671 WGS84Longitude, WGS84Latitude, WGS84Height );
01672
01673 return targetGeodeticCoordinates;
01674 }
01675 }
01676 }
01677 default:
01678 throw CoordinateConversionException( ErrorMessages::datumType );
01679 }
01680 }
01681
01682
01683 GeodeticCoordinates* DatumLibraryImplementation::geodeticShiftToWGS84( const long sourceIndex, const GeodeticCoordinates* sourceCoordinates )
01684 {
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699 double WGS84_a;
01700 double WGS84_f;
01701 double a;
01702 double da;
01703 double f;
01704 double df;
01705 double dx;
01706 double dy;
01707 double dz;
01708 long E_Index;
01709 char errorStatus[50] = "";
01710
01711 double sourceLongitude = sourceCoordinates->longitude();
01712 double sourceLatitude = sourceCoordinates->latitude();
01713 double sourceHeight = sourceCoordinates->height();
01714
01715 if( ( sourceIndex < 0 ) || ( sourceIndex >= datumList.size() ) )
01716 strcat( errorStatus, ErrorMessages::invalidIndex );
01717 if( ( sourceLatitude < ( -90 * PI_OVER_180 ) ) || ( sourceLatitude > ( 90 * PI_OVER_180 ) ) )
01718 strcat( errorStatus, ErrorMessages::latitude );
01719 if( ( sourceLongitude < ( -PI ) ) || ( sourceLongitude > TWO_PI ) )
01720 strcat( errorStatus, ErrorMessages::longitude );
01721
01722 if( strlen( errorStatus ) > 0)
01723 throw CoordinateConversionException( errorStatus );
01724
01725 Datum* localDatum = datumList[sourceIndex];
01726 switch( localDatum->datumType() )
01727 {
01728 case DatumType::wgs72Datum:
01729 {
01730 return geodeticShiftWGS72ToWGS84( sourceLongitude, sourceLatitude, sourceHeight );
01731 }
01732 case DatumType::wgs84Datum:
01733 {
01734 return new GeodeticCoordinates(CoordinateType::geodetic, sourceLongitude, sourceLatitude, sourceHeight);
01735 }
01736 case DatumType::sevenParamDatum:
01737 case DatumType::threeParamDatum:
01738 {
01739 if( _ellipsoidLibraryImplementation )
01740 {
01741 _ellipsoidLibraryImplementation->ellipsoidIndex( localDatum->ellipsoidCode(), &E_Index );
01742 _ellipsoidLibraryImplementation->ellipsoidParameters( E_Index, &a, &f );
01743
01744
01745
01746
01747
01748
01749
01750
01751 if( ( localDatum->datumType() == DatumType::sevenParamDatum ) ||
01752 ( sourceLatitude < (-MOLODENSKY_MAX ) ) ||
01753 ( sourceLatitude > MOLODENSKY_MAX ) )
01754 {
01755 Geocentric geocentricFromGeodetic( a, f );
01756
01757 CartesianCoordinates* localCartesianCoordinates = geocentricFromGeodetic.convertFromGeodetic( sourceCoordinates );
01758
01759 CartesianCoordinates* wgs84CartesianCoordinates = geocentricShiftToWGS84( sourceIndex, localCartesianCoordinates->x(), localCartesianCoordinates->y(), localCartesianCoordinates->z() );
01760
01761 long wgs84EllipsoidIndex;
01762 _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
01763 _ellipsoidLibraryImplementation->ellipsoidParameters( wgs84EllipsoidIndex, &WGS84_a, &WGS84_f );
01764
01765 Geocentric geocentricToGeodetic( WGS84_a, WGS84_f );
01766
01767 GeodeticCoordinates* wgs84GeodeticCoordinates = geocentricToGeodetic.convertToGeodetic( wgs84CartesianCoordinates );
01768
01769 delete wgs84CartesianCoordinates;
01770 delete localCartesianCoordinates;
01771
01772 return wgs84GeodeticCoordinates;
01773 }
01774 else
01775 {
01776 long wgs84EllipsoidIndex;
01777 _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
01778 _ellipsoidLibraryImplementation->ellipsoidParameters( wgs84EllipsoidIndex, &WGS84_a, &WGS84_f );
01779
01780 da = WGS84_a - a;
01781 df = WGS84_f - f;
01782 dx = localDatum->deltaX();
01783 dy = localDatum->deltaY();
01784 dz = localDatum->deltaZ();
01785
01786 GeodeticCoordinates* wgs84GeodeticCoordinates = molodenskyShift( a, da, f, df, dx, dy, dz, sourceLongitude, sourceLatitude, sourceHeight );
01787
01788 return wgs84GeodeticCoordinates;
01789 }
01790 }
01791
01792
01793
01794
01795 }
01796 default:
01797 throw CoordinateConversionException( ErrorMessages::datumType );
01798 }
01799 }
01800
01801
01802 void DatumLibraryImplementation::retrieveDatumType( const long index, DatumType::Enum *datumType )
01803 {
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813 if( index < 0 || index >= datumList.size() )
01814 throw CoordinateConversionException( ErrorMessages::invalidIndex );
01815 else
01816 *datumType = datumList[index]->datumType();
01817 }
01818
01819
01820 void DatumLibraryImplementation::validDatum( const long index, double longitude, double latitude, long *result )
01821 {
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834 *result = 0;
01835
01836 if( ( index < 0 ) || ( index >= datumList.size() ) )
01837 throw CoordinateConversionException( ErrorMessages::invalidIndex );
01838 if( ( latitude < MIN_LAT ) || ( latitude > MAX_LAT ) )
01839 throw CoordinateConversionException( ErrorMessages::latitude );
01840 if( ( longitude < MIN_LON ) || ( longitude > MAX_LON ) )
01841 throw CoordinateConversionException( ErrorMessages::longitude );
01842
01843 Datum* datum = datumList[index];
01844
01845 double west_longitude = datum->westLongitude();
01846 double east_longitude = datum->eastLongitude();
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856 if( ( west_longitude < 0 ) || ( east_longitude < 0 ) )
01857 {
01858 if( west_longitude > east_longitude )
01859 {
01860 if( west_longitude < 0 )
01861 west_longitude += TWO_PI;
01862 if( east_longitude < 0 )
01863 east_longitude += TWO_PI;
01864 if( longitude < 0 )
01865 longitude += TWO_PI;
01866 }
01867 }
01868 else if( ( west_longitude > PI ) || ( east_longitude > PI ) )
01869 {
01870 if( west_longitude > east_longitude )
01871 {
01872 if( west_longitude > PI )
01873 west_longitude -= TWO_PI;
01874 if( east_longitude > PI )
01875 east_longitude -= TWO_PI;
01876 }
01877 else
01878 {
01879 if( longitude < 0 )
01880 longitude += TWO_PI;
01881 }
01882 }
01883
01884 if( ( datum->southLatitude() <= latitude ) &&
01885 ( latitude <= datum->northLatitude() ) &&
01886 ( west_longitude <= longitude ) &&
01887 ( longitude <= east_longitude ) )
01888 {
01889 *result = 1;
01890 }
01891 else
01892 {
01893 *result = 0;
01894 }
01895 }
01896
01897
01898 void DatumLibraryImplementation::setEllipsoidLibraryImplementation( EllipsoidLibraryImplementation* __ellipsoidLibraryImplementation )
01899 {
01900
01901
01902
01903
01904
01905
01906
01907
01908 _ellipsoidLibraryImplementation = __ellipsoidLibraryImplementation;
01909 }
01910
01911
01912
01913
01914
01915
01916
01917 void DatumLibraryImplementation::loadDatums()
01918 {
01919
01920
01921
01922
01923
01924
01925
01926 long index = 0, i = 0;
01927 char *PathName = NULL;
01928 char* FileName7 = 0;
01929 FILE *fp_7param = NULL;
01930 FILE *fp_3param = NULL;
01931 char* FileName3 = 0;
01932 const int file_name_length = 23;
01933
01934 CCSThreadLock lock(&mutex);
01935
01936
01937
01938
01939 PathName = getenv( "MSPCCS_DATA" );
01940 if (PathName != NULL)
01941 {
01942 FileName7 = new char[ strlen( PathName ) + 13 ];
01943 strcpy( FileName7, PathName );
01944 strcat( FileName7, "/" );
01945 }
01946 else
01947 {
01948 FileName7 = new char[ file_name_length ];
01949 strcpy( FileName7, "../../data/" );
01950 }
01951 strcat( FileName7, "7_param.dat" );
01952
01953
01954
01955 if (( fp_7param = fopen( FileName7, "r" ) ) == NULL)
01956 {
01957 delete [] FileName7;
01958 FileName7 = 0;
01959
01960 char message[256] = "";
01961 strcpy( message, ErrorMessages::datumFileOpenError );
01962 strcat( message, ": 7_param.dat\n" );
01963 }
01964
01965 if (PathName != NULL)
01966 {
01967 FileName3 = new char[ strlen( PathName ) + 13 ];
01968 strcpy( FileName3, PathName );
01969 strcat( FileName3, "/" );
01970 }
01971 else
01972 {
01973 FileName3 = new char[ file_name_length ];
01974 strcpy( FileName3, "../../data/" );
01975 }
01976 strcat( FileName3, "3_param.dat" );
01977
01978
01979
01980 if (( fp_3param = fopen( FileName3, "r" ) ) == NULL)
01981 {
01982 delete [] FileName7;
01983 FileName7 = 0;
01984 delete [] FileName3;
01985 FileName3 = 0;
01986
01987 char message[256] = "";
01988 strcpy( message, ErrorMessages::datumFileOpenError );
01989 strcat( message, ": 3_param.dat\n" );
01990 throw CoordinateConversionException( message );
01991 }
01992
01993
01994 datumList.push_back( new Datum( index, "WGE", "WE", "World Geodetic System 1984", DatumType::wgs84Datum, 0.0, 0.0, 0.0, -PI, +PI, -PI / 2.0, +PI / 2.0, false ) );
01995 index ++;
01996
01997
01998
01999 datumList.push_back( new Datum( index, "WGC", "WD", "World Geodetic System 1972", DatumType::wgs72Datum, 0.0, 0.0, 0.0, -PI, +PI, -PI / 2.0, +PI / 2.0, false ) );
02000
02001 index ++;
02002
02003 datum7ParamCount = 0;
02004
02005 while ( !feof(fp_7param ) )
02006 {
02007 if ( datum7ParamCount < MAX_7PARAM )
02008 {
02009 char code[DATUM_CODE_LENGTH];
02010
02011 bool userDefined = false;
02012
02013 if (fscanf(fp_7param, "%s ", code) <= 0)
02014 throw CoordinateConversionException( ErrorMessages::datumFileParseError );
02015 else
02016 {
02017 if( code[0] == '*' )
02018 {
02019 userDefined = true;
02020 for( int i = 0; i < DATUM_CODE_LENGTH; i++ )
02021 code[i] = code[i+1];
02022 }
02023 }
02024
02025 char name[DATUM_NAME_LENGTH];
02026 if (fscanf(fp_7param, "\"%32[^\"]\"", name) <= 0)
02027 name[0] = '\0';
02028
02029 char ellipsoidCode[ELLIPSOID_CODE_LENGTH];
02030 double deltaX;
02031 double deltaY;
02032 double deltaZ;
02033 double rotationX;
02034 double rotationY;
02035 double rotationZ;
02036 double scaleFactor;
02037
02038 if( fscanf( fp_7param, " %s %lf %lf %lf %lf %lf %lf %lf ",
02039 ellipsoidCode, &deltaX, &deltaY, &deltaZ, &rotationX, &rotationY, &rotationZ, &scaleFactor ) <= 0 )
02040 {
02041 throw CoordinateConversionException( ErrorMessages::datumFileParseError );
02042 }
02043 else
02044 {
02045
02046 rotationX /= SECONDS_PER_RADIAN;
02047 rotationY /= SECONDS_PER_RADIAN;
02048 rotationZ /= SECONDS_PER_RADIAN;
02049
02050 datumList.push_back( new SevenParameterDatum( index, code, ellipsoidCode, name, DatumType::sevenParamDatum, deltaX, deltaY, deltaZ, -PI, +PI, -PI / 2.0, +PI / 2.0,
02051 rotationX, rotationY, rotationZ, scaleFactor, userDefined ) );
02052 }
02053 index++;
02054 datum7ParamCount++;
02055 }
02056 else
02057 {
02058 throw CoordinateConversionException( ErrorMessages::datumOverflow );
02059 }
02060 }
02061 fclose( fp_7param );
02062
02063 datum3ParamCount = 0;
02064
02065
02066 while( !feof( fp_3param ) )
02067 {
02068 if( datum3ParamCount < MAX_3PARAM )
02069 {
02070 char code[DATUM_CODE_LENGTH];
02071
02072 bool userDefined = false;
02073
02074 if( fscanf( fp_3param, "%s ", code ) <= 0 )
02075 throw CoordinateConversionException( ErrorMessages::datumFileParseError );
02076 else
02077 {
02078 if( code[0] == '*' )
02079 {
02080 userDefined = true;
02081 for( int i = 0; i < DATUM_CODE_LENGTH; i++ )
02082 code[i] = code[i+1];
02083 }
02084 }
02085
02086 char name[DATUM_NAME_LENGTH];
02087 if( fscanf( fp_3param, "\"%32[^\"]\"", name) <= 0 )
02088 name[0] = '\0';
02089
02090 char ellipsoidCode[ELLIPSOID_CODE_LENGTH];
02091 double deltaX;
02092 double deltaY;
02093 double deltaZ;
02094 double sigmaX;
02095 double sigmaY;
02096 double sigmaZ;
02097 double eastLongitude;
02098 double westLongitude;
02099 double northLatitude;
02100 double southLatitude;
02101
02102 if( fscanf( fp_3param, " %s %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf ",
02103 ellipsoidCode, &deltaX, &sigmaX, &deltaY, &sigmaY, &deltaZ, &sigmaZ,
02104 &southLatitude, &northLatitude, &westLongitude, &eastLongitude ) <= 0 )
02105 {
02106 throw CoordinateConversionException( ErrorMessages::datumFileParseError );
02107 }
02108 else
02109 {
02110 southLatitude *= PI_OVER_180;
02111 northLatitude *= PI_OVER_180;
02112 westLongitude *= PI_OVER_180;
02113 eastLongitude *= PI_OVER_180;
02114
02115 datumList.push_back( new ThreeParameterDatum( index, code, ellipsoidCode, name, DatumType::threeParamDatum, deltaX, deltaY, deltaZ,
02116 westLongitude, eastLongitude, southLatitude, northLatitude, sigmaX, sigmaY, sigmaZ, userDefined ) );
02117 }
02118
02119 index++;
02120 datum3ParamCount++;
02121 }
02122 else
02123 {
02124 throw CoordinateConversionException( ErrorMessages::datumOverflow );
02125 }
02126 }
02127 fclose( fp_3param );
02128
02129 delete [] FileName7;
02130 FileName7 = 0;
02131 delete [] FileName3;
02132 FileName3 = 0;
02133 }
02134
02135
02136 void DatumLibraryImplementation::write3ParamFile()
02137 {
02138
02139
02140
02141
02142
02143 char datum_name[DATUM_NAME_LENGTH+2];
02144 char *PathName = NULL;
02145 char FileName[FILENAME_LENGTH];
02146 FILE *fp_3param = NULL;
02147
02148 CCSThreadLock lock(&mutex);
02149
02150
02151 PathName = getenv( "MSPCCS_DATA" );
02152 if (PathName != NULL)
02153 {
02154 strcpy( FileName, PathName );
02155 strcat( FileName, "/" );
02156 }
02157 else
02158 {
02159 strcpy( FileName, "../../data/" );
02160 }
02161 strcat( FileName, "3_param.dat" );
02162
02163 if ((fp_3param = fopen(FileName, "w")) == NULL)
02164 {
02165 char message[256] = "";
02166 strcpy( message, ErrorMessages::datumFileOpenError );
02167 strcat( message, ": 3_param.dat\n" );
02168 throw CoordinateConversionException( message );
02169 }
02170
02171
02172 long index = MAX_WGS + datum7ParamCount;
02173 int size = datumList.size();
02174 while( index < size )
02175 {
02176 ThreeParameterDatum* _3parameterDatum = ( ThreeParameterDatum* )datumList[index];
02177 if( _3parameterDatum )
02178 {
02179 strcpy( datum_name, "\"" );
02180 strcat( datum_name, datumList[index]->name());
02181 strcat( datum_name, "\"" );
02182 if( _3parameterDatum->userDefined() )
02183 fprintf( fp_3param, "*");
02184
02185 fprintf( fp_3param, "%-6s %-33s%-2s %4.0f %4.0f %4.0f %4.0f %5.0f %4.0f %4.0f %4.0f %4.0f %4.0f \n",
02186 _3parameterDatum->code(),
02187 datum_name,
02188 _3parameterDatum->ellipsoidCode(),
02189 _3parameterDatum->deltaX(),
02190 _3parameterDatum->sigmaX(),
02191 _3parameterDatum->deltaY(),
02192 _3parameterDatum->sigmaY(),
02193 _3parameterDatum->deltaZ(),
02194 _3parameterDatum->sigmaZ(),
02195 ( _3parameterDatum->southLatitude() * _180_OVER_PI ),
02196 ( _3parameterDatum->northLatitude() * _180_OVER_PI ),
02197 ( _3parameterDatum->westLongitude() * _180_OVER_PI ),
02198 ( _3parameterDatum->eastLongitude() * _180_OVER_PI ) );
02199 }
02200
02201 index++;
02202 }
02203
02204 fclose( fp_3param );
02205
02206 }
02207
02208
02209 void DatumLibraryImplementation::write7ParamFile()
02210 {
02211
02212
02213
02214
02215
02216 char datum_name[DATUM_NAME_LENGTH+2];
02217 char *PathName = NULL;
02218 char FileName[FILENAME_LENGTH];
02219 FILE *fp_7param = NULL;
02220
02221 CCSThreadLock lock(&mutex);
02222
02223
02224 PathName = getenv( "MSPCCS_DATA" );
02225 if (PathName != NULL)
02226 {
02227 strcpy( FileName, PathName );
02228 strcat( FileName, "/" );
02229 }
02230 else
02231 {
02232 strcpy( FileName, "../../data/" );
02233 }
02234 strcat( FileName, "7_param.dat" );
02235
02236 if ((fp_7param = fopen(FileName, "w")) == NULL)
02237 {
02238 char message[256] = "";
02239 strcpy( message, ErrorMessages::datumFileOpenError );
02240 strcat( message, ": 7_param.dat\n" );
02241 throw CoordinateConversionException( message );
02242 }
02243
02244
02245 long index = MAX_WGS;
02246 int endIndex = datum7ParamCount + MAX_WGS;
02247 while( index < endIndex )
02248 {
02249 SevenParameterDatum* _7parameterDatum = ( SevenParameterDatum* )datumList[index];
02250 if( _7parameterDatum )
02251 {
02252 strcpy( datum_name, "\"" );
02253 strcat( datum_name, datumList[index]->name());
02254 strcat( datum_name, "\"" );
02255 if( _7parameterDatum->userDefined() )
02256 fprintf( fp_7param, "*");
02257
02258 fprintf( fp_7param, "%-6s %-33s%-2s %4.0f %4.0f %4.0f % 4.3f % 4.3f % 4.3f % 4.10f \n",
02259 _7parameterDatum->code(),
02260 datum_name,
02261 _7parameterDatum->ellipsoidCode(),
02262 _7parameterDatum->deltaX(),
02263 _7parameterDatum->deltaY(),
02264 _7parameterDatum->deltaZ(),
02265 _7parameterDatum->rotationX() * SECONDS_PER_RADIAN,
02266 _7parameterDatum->rotationY() * SECONDS_PER_RADIAN,
02267 _7parameterDatum->rotationZ() * SECONDS_PER_RADIAN,
02268 _7parameterDatum->scaleFactor() );
02269 }
02270
02271 index++;
02272 }
02273
02274 fclose( fp_7param );
02275 }
02276
02277
02278 GeodeticCoordinates* DatumLibraryImplementation::geodeticShiftWGS84ToWGS72( const double WGS84Longitude, const double WGS84Latitude, const double WGS84Height )
02279 {
02280
02281
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291
02292
02293 double Delta_Lat;
02294 double Delta_Lon;
02295 double Delta_Hgt;
02296 double WGS84_a;
02297 double WGS84_f;
02298 double WGS72_a;
02299 double WGS72_f;
02300 double da;
02301 double df;
02302 double Q;
02303 double sin_Lat;
02304 double sin2_Lat;
02305
02306 long wgs84EllipsoidIndex;
02307 _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
02308 _ellipsoidLibraryImplementation->ellipsoidParameters( wgs84EllipsoidIndex, &WGS84_a, &WGS84_f );
02309
02310 long wgs72EllipsoidIndex;
02311 _ellipsoidLibraryImplementation->ellipsoidIndex( "WD", &wgs72EllipsoidIndex );
02312 _ellipsoidLibraryImplementation->ellipsoidParameters( wgs72EllipsoidIndex, &WGS72_a, &WGS72_f );
02313
02314 da = WGS72_a - WGS84_a;
02315 df = WGS72_f - WGS84_f;
02316 Q = PI / 648000;
02317 sin_Lat = sin(WGS84Latitude);
02318 sin2_Lat = sin_Lat * sin_Lat;
02319
02320 Delta_Lat = (-4.5 * cos(WGS84Latitude)) / (WGS84_a*Q)
02321 + (df * sin(2*WGS84Latitude)) / Q;
02322 Delta_Lat /= SECONDS_PER_RADIAN;
02323 Delta_Lon = -0.554 / SECONDS_PER_RADIAN;
02324 Delta_Hgt = -4.5 * sin_Lat + WGS84_a * df * sin2_Lat - da - 1.4;
02325
02326 double wgs72Latitude = WGS84Latitude + Delta_Lat;
02327 double wgs72Longitude = WGS84Longitude + Delta_Lon;
02328 double wgs72Height = WGS84Height + Delta_Hgt;
02329
02330 if (wgs72Latitude > PI_OVER_2)
02331 wgs72Latitude = PI_OVER_2 - (wgs72Latitude - PI_OVER_2);
02332 else if (wgs72Latitude < -PI_OVER_2)
02333 wgs72Latitude = -PI_OVER_2 - (wgs72Latitude + PI_OVER_2);
02334
02335 if (wgs72Longitude > PI)
02336 wgs72Longitude -= TWO_PI;
02337 if (wgs72Longitude < -PI)
02338 wgs72Longitude += TWO_PI;
02339
02340 return new GeodeticCoordinates(CoordinateType::geodetic, wgs72Longitude, wgs72Latitude, wgs72Height);
02341 }
02342
02343
02344 GeodeticCoordinates* DatumLibraryImplementation::geodeticShiftWGS72ToWGS84( const double WGS72Longitude, const double WGS72Latitude, const double WGS72Height )
02345 {
02346
02347
02348
02349
02350
02351
02352
02353
02354
02355
02356
02357
02358
02359 double Delta_Lat;
02360 double Delta_Lon;
02361 double Delta_Hgt;
02362 double WGS84_a;
02363 double WGS84_f;
02364 double WGS72_a;
02365 double WGS72_f;
02366 double da;
02367 double df;
02368 double Q;
02369 double sin_Lat;
02370 double sin2_Lat;
02371
02372 long wgs84EllipsoidIndex;
02373 _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
02374 _ellipsoidLibraryImplementation->ellipsoidParameters( wgs84EllipsoidIndex, &WGS84_a, &WGS84_f );
02375
02376 long wgs72EllipsoidIndex;
02377 _ellipsoidLibraryImplementation->ellipsoidIndex( "WD", &wgs72EllipsoidIndex );
02378 _ellipsoidLibraryImplementation->ellipsoidParameters( wgs72EllipsoidIndex, &WGS72_a, &WGS72_f );
02379
02380 da = WGS84_a - WGS72_a;
02381 df = WGS84_f - WGS72_f;
02382 Q = PI / 648000;
02383 sin_Lat = sin(WGS72Latitude);
02384 sin2_Lat = sin_Lat * sin_Lat;
02385
02386 Delta_Lat = (4.5 * cos(WGS72Latitude)) / (WGS72_a*Q) + (df * sin(2*WGS72Latitude)) / Q;
02387 Delta_Lat /= SECONDS_PER_RADIAN;
02388 Delta_Lon = 0.554 / SECONDS_PER_RADIAN;
02389 Delta_Hgt = 4.5 * sin_Lat + WGS72_a * df * sin2_Lat - da + 1.4;
02390
02391 double wgs84Latitude = WGS72Latitude + Delta_Lat;
02392 double wgs84Longitude = WGS72Longitude + Delta_Lon;
02393 double wgs84Height = WGS72Height + Delta_Hgt;
02394
02395 if (wgs84Latitude > PI_OVER_2)
02396 wgs84Latitude = PI_OVER_2 - (wgs84Latitude - PI_OVER_2);
02397 else if (wgs84Latitude < -PI_OVER_2)
02398 wgs84Latitude = -PI_OVER_2 - (wgs84Latitude + PI_OVER_2);
02399
02400 if (wgs84Longitude > PI)
02401 wgs84Longitude -= TWO_PI;
02402 if (wgs84Longitude < -PI)
02403 wgs84Longitude += TWO_PI;
02404
02405 return new GeodeticCoordinates(CoordinateType::geodetic, wgs84Longitude, wgs84Latitude, wgs84Height);
02406 }
02407
02408
02409 CartesianCoordinates* DatumLibraryImplementation::geocentricShiftWGS84ToWGS72( const double X_WGS84, const double Y_WGS84, const double Z_WGS84 )
02410 {
02411
02412
02413
02414
02415
02416
02417
02418
02419
02420
02421
02422
02423 double a_72;
02424 double f_72;
02425 double a_84;
02426 double f_84;
02427
02428
02429 long wgs84EllipsoidIndex;
02430 _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
02431 _ellipsoidLibraryImplementation->ellipsoidParameters( wgs84EllipsoidIndex, &a_84, &f_84 );
02432
02433 Geocentric geocentric84( a_84, f_84 );
02434
02435
02436 GeodeticCoordinates* wgs84GeodeticCoordinates = geocentric84.convertToGeodetic( new CartesianCoordinates( CoordinateType::geocentric, X_WGS84, Y_WGS84, Z_WGS84 ) );
02437
02438 GeodeticCoordinates* wgs72GeodeticCoordinates = geodeticShiftWGS84ToWGS72( wgs84GeodeticCoordinates->longitude(), wgs84GeodeticCoordinates->latitude(), wgs84GeodeticCoordinates->height() );
02439
02440
02441 long wgs72EllipsoidIndex;
02442 _ellipsoidLibraryImplementation->ellipsoidIndex( "WD", &wgs72EllipsoidIndex );
02443 _ellipsoidLibraryImplementation->ellipsoidParameters( wgs72EllipsoidIndex, &a_72, &f_72 );
02444
02445 Geocentric geocentric72( a_72, f_72 );
02446
02447
02448 CartesianCoordinates* wgs72GeocentricCoordinates = geocentric72.convertFromGeodetic( wgs72GeodeticCoordinates );
02449
02450 delete wgs72GeodeticCoordinates;
02451 delete wgs84GeodeticCoordinates;
02452
02453 return wgs72GeocentricCoordinates;
02454 }
02455
02456
02457 CartesianCoordinates* DatumLibraryImplementation::geocentricShiftWGS72ToWGS84( const double X, const double Y, const double Z )
02458 {
02459
02460
02461
02462
02463
02464
02465
02466
02467
02468
02469
02470
02471 double a_72;
02472 double f_72;
02473 double a_84;
02474 double f_84;
02475
02476
02477 long wgs72EllipsoidIndex;
02478 _ellipsoidLibraryImplementation->ellipsoidIndex( "WD", &wgs72EllipsoidIndex );
02479 _ellipsoidLibraryImplementation->ellipsoidParameters( wgs72EllipsoidIndex, &a_72, &f_72 );
02480
02481 Geocentric geocentric72( a_72, f_72 );
02482
02483 GeodeticCoordinates* wgs72GeodeticCoordinates = geocentric72.convertToGeodetic( new CartesianCoordinates( CoordinateType::geocentric, X, Y, Z ) );
02484
02485 GeodeticCoordinates* wgs84GeodeticCoordinates = geodeticShiftWGS72ToWGS84( wgs72GeodeticCoordinates->longitude(), wgs72GeodeticCoordinates->latitude(), wgs72GeodeticCoordinates->height() );
02486
02487
02488 long wgs84EllipsoidIndex;
02489 _ellipsoidLibraryImplementation->ellipsoidIndex( "WE", &wgs84EllipsoidIndex );
02490 _ellipsoidLibraryImplementation->ellipsoidParameters( wgs84EllipsoidIndex, &a_84, &f_84 );
02491
02492 Geocentric geocentric84( a_84, f_84 );
02493
02494 CartesianCoordinates* wgs84GeocentricCoordinates = geocentric84.convertFromGeodetic( wgs84GeodeticCoordinates );
02495
02496 delete wgs84GeodeticCoordinates;
02497 delete wgs72GeodeticCoordinates;
02498
02499 return wgs84GeocentricCoordinates;
02500 }
02501
02502
02503