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