00001
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
00033
00034
00035
00036
00037 #include <fstream>
00038
00039 #ifdef IRIXN32
00040 #include <math.h>
00041 #else
00042 #include <cmath>
00043 #endif
00044
00045 #include "CCSThreadLock.h"
00046 #include "CoordinateConversionException.h"
00047
00048 #include "egm2008_aoi_grid_package.h"
00049
00050 namespace
00051 {
00052 const int BYTES_PER_DOUBLE = sizeof( double );
00053 const int BYTES_PER_FLOAT = sizeof( float );
00054 const int BYTES_PER_INT = sizeof( int );
00055
00056 const int BYTES_IN_HEADER =
00057 ( 3 * BYTES_PER_INT + 2 * BYTES_PER_DOUBLE );
00058
00059 const int NOM_AOI_COLS1 = 100;
00060 const int NOM_AOI_ROWS1 = 100;
00061 const int NOM_AOI_COLS2_5 = 50;
00062 const int NOM_AOI_ROWS2_5 = 50;
00063
00064 const double PI = 3.14159265358979323;
00065
00066 const double PIDIV2 = PI / 2.0;
00067 const double TWOPI = 2.0 * PI;
00068
00069 const double DEG_PER_RAD = 180.0 / PI;
00070 const double RAD_PER_DEG = PI / 180.0;
00071
00072 const double DLAT1 = RAD_PER_DEG * (1.0 / 60.0);
00073 const double DLON1 = RAD_PER_DEG * (1.0 / 60.0);
00074 const double DLAT2_5 = RAD_PER_DEG * (2.5 / 60.0);
00075 const double DLON2_5 = RAD_PER_DEG * (2.5 / 60.0);
00076
00077 const double EPSILON = 1.0e-15;
00078
00079 const double SEMI_MAJOR_AXIS = 6378137.0;
00080 const double FLATTENING = 1.0 / 298.257223563;
00081 }
00082
00083 using namespace MSP;
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093 Egm2008AoiGrid::Egm2008AoiGrid (void)
00094 {
00095
00096
00097
00098 int status;
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108 _minAoiRowIndex = 10000000;
00109 _maxAoiRowIndex = -10000000;
00110 _minAoiColIndex = 10000000;
00111 _maxAoiColIndex = -10000000;
00112
00113 _nAoiCols = 0;
00114 _nAoiRows = 0;
00115
00116 _nomAoiCols = 0;
00117 _nomAoiRows = 0;
00118
00119 _heightGrid = NULL;
00120
00121
00122
00123 status = this->loadGridMetadata();
00124
00125 if ( status != 0 )
00126 {
00127 throw MSP::CCS::CoordinateConversionException(
00128 "Error: Egm2008AoiGrid: constructor failed." );
00129 }
00130
00131 }
00132
00133
00134
00135
00136
00137
00138 Egm2008AoiGrid::Egm2008AoiGrid
00139 ( const Egm2008AoiGrid& oldGrid )
00140
00141 : Egm2008GeoidGrid( oldGrid )
00142
00143 {
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153 int i;
00154 int kount;
00155
00156
00157
00158 _maxAoiColIndex = oldGrid._maxAoiColIndex;
00159 _minAoiColIndex = oldGrid._minAoiColIndex;
00160
00161 _maxAoiRowIndex = oldGrid._maxAoiRowIndex;
00162 _minAoiRowIndex = oldGrid._minAoiRowIndex;
00163
00164 _nAoiCols = oldGrid._nAoiCols;
00165 _nAoiRows = oldGrid._nAoiRows;
00166
00167 _nomAoiCols = oldGrid._nomAoiCols;
00168 _nomAoiRows = oldGrid._nomAoiRows;
00169
00170 _heightGrid = NULL;
00171
00172
00173
00174
00175
00176
00177
00178 try
00179 {
00180 if ( NULL != oldGrid._heightGrid )
00181 {
00182 kount = _nAoiRows * _nAoiCols;
00183
00184 _heightGrid = new float[ kount ];
00185
00186 for ( i = 0; i < kount; i++ )
00187 {
00188 _heightGrid[ i ] = oldGrid._heightGrid[ i ];
00189 }
00190 }
00191 }
00192 catch (...)
00193 {
00194 oldGrid._mutex.unlock();
00195
00196 throw MSP::CCS::CoordinateConversionException(
00197 "Error: Egm2008AoiGrid: copy contructor failed");
00198 }
00199
00200 oldGrid._mutex.unlock();
00201
00202 }
00203
00204
00205
00206
00207
00208
00209 Egm2008AoiGrid::~Egm2008AoiGrid (void)
00210 {
00211
00212
00213
00214
00215
00216 delete[] _heightGrid;
00217
00218 }
00219
00220
00221
00222
00223
00224
00225 Egm2008AoiGrid&
00226 Egm2008AoiGrid::operator= ( const Egm2008AoiGrid& oldGrid )
00227 {
00228
00229
00230
00231
00232
00233
00234 int i;
00235 int kount;
00236
00237
00238
00239
00240 if ( this == & oldGrid ) return ( *this );
00241
00242
00243
00244 Egm2008GeoidGrid::operator= ( oldGrid );
00245
00246
00247
00248 _maxAoiColIndex = oldGrid._maxAoiColIndex;
00249 _minAoiColIndex = oldGrid._minAoiColIndex;
00250
00251 _maxAoiRowIndex = oldGrid._maxAoiRowIndex;
00252 _minAoiRowIndex = oldGrid._minAoiRowIndex;
00253
00254 _nAoiCols = oldGrid._nAoiCols;
00255 _nAoiRows = oldGrid._nAoiRows;
00256
00257 _nomAoiCols = oldGrid._nomAoiCols;
00258 _nomAoiRows = oldGrid._nomAoiRows;
00259
00260
00261
00262
00263
00264
00265
00266 try
00267 {
00268 delete[] _heightGrid; _heightGrid = NULL;
00269
00270 if ( NULL != oldGrid._heightGrid )
00271 {
00272 kount = _nAoiRows * _nAoiCols;
00273
00274 _heightGrid = new float[ kount ];
00275
00276 for ( i = 0; i < kount; i++ )
00277 {
00278 _heightGrid[ i ] = oldGrid._heightGrid[ i ];
00279 }
00280 }
00281 }
00282 catch (...)
00283 {
00284 _mutex.unlock();
00285 oldGrid._mutex.unlock();
00286
00287 throw MSP::CCS::CoordinateConversionException(
00288 "Error: Egm2008AoiGrid: assignment operator failed");
00289 }
00290
00291 _mutex.unlock();
00292 oldGrid._mutex.unlock();
00293
00294 return( *this );
00295
00296 }
00297
00298
00299
00300
00301
00302
00303 int
00304 Egm2008AoiGrid::geoidHeight(
00305 int wSize,
00306 double latitude,
00307 double longitude,
00308 double& gHeight )
00309 {
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359 try {
00360
00361 MSP::CCSThreadLock aLock( &_mutex );
00362
00363 const int TWENTY = 20;
00364
00365 bool oddSize;
00366
00367 int i;
00368 int i0;
00369 int iIndex;
00370 int iMax;
00371 int iMin;
00372 int j;
00373 int j0;
00374 int jIndex;
00375 int jMax;
00376 int jMin;
00377 int kIndex;
00378 int offset;
00379 int status;
00380
00381 double latIndex;
00382 double lonIndex;
00383 double temp;
00384
00385 double latSupport[ TWENTY ];
00386 double lonSupport[ TWENTY ];
00387 double moments [ TWENTY ];
00388
00389
00390
00391 gHeight = 0.0;
00392
00393 if ( wSize > MAX_WSIZE ) wSize = MAX_WSIZE;
00394
00395 if (( latitude < -PIDIV2 ) || ( latitude > PIDIV2 )) return( 1 );
00396
00397 while ( longitude < 0.0 ) longitude += TWOPI;
00398 while ( longitude > TWOPI ) longitude -= TWOPI;
00399
00400 if ( TWENTY != MAX_WSIZE ) return( 1 );
00401
00402 for (i = 0; i < TWENTY; i++)
00403 {
00404 latSupport[ i ] = lonSupport[ i ] = moments[ i ] = 0.0;
00405 }
00406
00407
00408
00409
00410 latIndex =
00411 double( _nGridPad ) +
00412 ( latitude + PIDIV2 ) / _dLat;
00413 lonIndex =
00414 double( _nGridPad ) +
00415 ( longitude - 0.0 ) / _dLon;
00416
00417
00418
00419
00420 oddSize = ( wSize != (( wSize / 2 ) * 2 ));
00421
00422 if ( oddSize ) {
00423 i0 = int( latIndex + 0.5 );
00424 j0 = int( lonIndex + 0.5 );
00425
00426 iMin = i0 - ( wSize / 2 );
00427 jMin = j0 - ( wSize / 2 );
00428 }
00429 else {
00430 i0 = int( latIndex );
00431 j0 = int( lonIndex );
00432
00433 iMin = i0 - ( wSize / 2 ) + 1;
00434 jMin = j0 - ( wSize / 2 ) + 1;
00435 }
00436
00437 iMax = iMin + wSize - 1;
00438 jMax = jMin + wSize - 1;
00439
00440
00441
00442
00443
00444 if (( iMin < _minAoiRowIndex ) ||
00445 ( iMax > _maxAoiRowIndex ) ||
00446 ( jMin < _minAoiColIndex ) ||
00447 ( jMax > _maxAoiColIndex ))
00448 {
00449 status = this->loadAoiParms( i0, j0 );
00450
00451 if ( status != 0 ) return( 1 );
00452
00453 status = this->loadGrid();
00454
00455 if ( status != 0 ) return( 1 );
00456 }
00457
00458
00459
00460
00461
00462
00463 if ( wSize < 3 )
00464 {
00465 status =
00466 this->geoidHeight( latitude, longitude, gHeight );
00467
00468 if ( status != 0 ) return( 1 );
00469
00470 ; return( 0 );
00471 }
00472
00473
00474
00475
00476
00477
00478
00479 temp =
00480 lonIndex - double( jMin );
00481
00482
00483
00484
00485 for ( i = 0; i < wSize; i++ )
00486 {
00487 iIndex = iMin + i - _minAoiRowIndex;
00488
00489 offset = ( _nAoiRows - iIndex - 1 ) * _nAoiCols;
00490
00491
00492
00493
00494 for ( j = 0; j < wSize; j++ )
00495 {
00496 jIndex = jMin + j - _minAoiColIndex;
00497
00498 kIndex = jIndex + offset;
00499
00500 lonSupport[ j ] = _heightGrid[ kIndex ];
00501 }
00502
00503
00504
00505
00506
00507
00508 status =
00509 this->initSpline( wSize, lonSupport, moments );
00510
00511 if ( status != 0 ) return( 1 );
00512
00513 latSupport[ i ] =
00514 this->spline( wSize, temp, lonSupport, moments );
00515 }
00516
00517
00518
00519
00520
00521 temp =
00522 latIndex - double( iMin );
00523
00524
00525
00526
00527 status =
00528 this->initSpline( wSize, latSupport, moments );
00529
00530 if ( status != 0 ) return( 1 );
00531
00532 gHeight =
00533 this->spline( wSize, temp, latSupport, moments );
00534
00535
00536
00537
00538 }
00539
00540 catch ( ... ) { return( 1 ); }
00541
00542 return ( 0 );
00543
00544 }
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555 int
00556 Egm2008AoiGrid::geoidHeight(
00557 double latitude,
00558 double longitude,
00559 double& gHeight )
00560 {
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582 try {
00583
00584 int i1;
00585 int i1ww;
00586 int i2;
00587 int i3;
00588 int i4;
00589 int index;
00590 int j1;
00591 int j1ww;
00592 int j2;
00593 int j3;
00594 int j4;
00595 int status;
00596
00597 double a0;
00598 double a1;
00599 double a2;
00600 double a3;
00601 double lat1;
00602 double lon1;
00603 double n1;
00604 double n2;
00605 double n3;
00606 double n4;
00607 double x;
00608 double y;
00609
00610
00611
00612 gHeight = 0.0;
00613
00614 if (( latitude < -PIDIV2 ) || ( latitude > PIDIV2 )) return( 1 );
00615
00616 while ( longitude < 0.0 ) longitude += TWOPI;
00617 while ( longitude > TWOPI ) longitude -= TWOPI;
00618
00619
00620
00621
00622 status =
00623 this->swGridIndices
00624 ( latitude, longitude, i1, j1 );
00625
00626 if ( status != 0 ) return( 1 );
00627
00628 i1ww = i1;
00629 j1ww = j1;
00630
00631 i2 = i1;
00632 j2 = j1 + 1;
00633
00634 i3 = i1 + 1;
00635 j3 = j2;
00636
00637 i4 = i3;
00638 j4 = j1;
00639
00640
00641
00642
00643 i1 -= _minAoiRowIndex;
00644 j1 -= _minAoiColIndex;
00645
00646 i2 -= _minAoiRowIndex;
00647 j2 -= _minAoiColIndex;
00648
00649 i3 -= _minAoiRowIndex;
00650 j3 -= _minAoiColIndex;
00651
00652 i4 -= _minAoiRowIndex;
00653 j4 -= _minAoiColIndex;
00654
00655
00656
00657 index =
00658 j1 + ( _nAoiRows - i1 - 1 ) * _nAoiCols;
00659 n1 = _heightGrid[ index ];
00660
00661 index =
00662 j2 + ( _nAoiRows - i2 - 1 ) * _nAoiCols;
00663 n2 = _heightGrid[ index ];
00664
00665 index =
00666 j3 + ( _nAoiRows - i3 - 1 ) * _nAoiCols;
00667 n3 = _heightGrid[ index ];
00668
00669 index =
00670 j4 + ( _nAoiRows - i4 - 1 ) * _nAoiCols;
00671 n4 = _heightGrid[ index ];
00672
00673
00674
00675
00676
00677 a0 = n1;
00678 a1 = n2 - n1;
00679 a2 = n4 - n1;
00680 a3 = n1 + n3 - n2 - n4;
00681
00682
00683
00684 lat1 =
00685 _baseLatitude + _dLat * double( i1ww );
00686 lon1 =
00687 _baseLongitude + _dLon * double( j1ww );
00688
00689
00690
00691
00692
00693
00694 x = ( longitude - lon1 ) / _dLon;
00695 y = ( latitude - lat1 ) / _dLat;
00696
00697
00698
00699 gHeight = a0 + a1 * x + a2 * y + a3 * x * y;
00700
00701 }
00702
00703 catch ( ... ) { return( 1 ); }
00704
00705 return( 0 );
00706
00707 }
00708
00709
00710
00711
00712
00713
00714 int
00715 Egm2008AoiGrid::loadAoiParms(
00716 int i0,
00717 int j0 )
00718 {
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747 try {
00748
00749 const double THRESHOLD = 0.05;
00750
00751 int status;
00752
00753 double cLat;
00754 double eSquared;
00755 double ewDelta;
00756 double ewDelta0;
00757 double latitude;
00758 double longitude;
00759 double nRadius;
00760 double ratio;
00761 double sLat;
00762 double temp;
00763
00764
00765
00766
00767
00768
00769 status =
00770 this->loadGridCoords(
00771 i0, j0, latitude, longitude );
00772
00773 if ( status != 0 ) return( 1 );
00774
00775 cLat = cos( latitude );
00776 sLat = sin( latitude );
00777
00778
00779
00780
00781 eSquared = FLATTENING * ( 2.0 - FLATTENING );
00782
00783 temp = 1.0 - eSquared * sLat * sLat;
00784
00785 nRadius = SEMI_MAJOR_AXIS / sqrt( temp );
00786
00787
00788
00789 ewDelta0 = SEMI_MAJOR_AXIS * _dLon;
00790
00791 ewDelta = nRadius * _dLon * cLat;
00792
00793 ratio = ewDelta / ewDelta0;
00794
00795 if ( ratio < THRESHOLD ) ratio = THRESHOLD;
00796
00797 _nAoiCols = int( double( _nomAoiCols ) / ratio );
00798
00799 _nAoiCols = 2 * ( _nAoiCols / 2 );
00800 _nAoiRows = 2 * ( _nomAoiRows / 2 );
00801
00802
00803
00804
00805
00806 _minAoiRowIndex = i0 - (( _nAoiRows - 2 ) / 2 ) + 1;
00807 _maxAoiRowIndex = _minAoiRowIndex + _nAoiRows - 1;
00808
00809 _minAoiColIndex = j0 - (( _nAoiCols - 2 ) / 2 ) + 1;
00810 _maxAoiColIndex = _minAoiColIndex + _nAoiCols - 1;
00811
00812
00813
00814
00815 if ( _minAoiRowIndex < 0 )
00816 {
00817 _minAoiRowIndex = 0;
00818 _maxAoiRowIndex = _minAoiRowIndex + _nAoiRows - 1;
00819 }
00820
00821 if ( _maxAoiRowIndex >= _nGridRows )
00822 {
00823 _maxAoiRowIndex = _nGridRows - 1;
00824 _minAoiRowIndex = _maxAoiRowIndex - _nAoiRows + 1;
00825 }
00826
00827 if ( _minAoiColIndex < 0 )
00828 {
00829 _minAoiColIndex = 0;
00830 _maxAoiColIndex = _minAoiColIndex + _nAoiCols - 1;
00831 }
00832
00833 if ( _maxAoiColIndex >= _nGridCols )
00834 {
00835 _maxAoiColIndex = _nGridCols - 1;
00836 _minAoiColIndex = _maxAoiColIndex - _nAoiCols + 1;
00837 }
00838
00839 }
00840
00841 catch ( ... ) { return( 1 ); }
00842
00843 return( 0 );
00844
00845 }
00846
00847
00848
00849
00850
00851
00852
00853 int
00854 Egm2008AoiGrid::loadGrid( void )
00855 {
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872 try {
00873
00874 int endRow;
00875 int i;
00876 int index0;
00877 int index1;
00878 int index2;
00879 int kount;
00880 int startCol;
00881 int startRow;
00882
00883 std::ifstream fin;
00884
00885
00886
00887 fin.open(
00888 _gridFname.c_str(),
00889 std::ifstream::binary | std::ifstream::in );
00890
00891 if ( fin.fail() ) return( 1 );
00892
00893 delete[] _heightGrid; _heightGrid = NULL;
00894
00895 _heightGrid = new float[ _nAoiRows * _nAoiCols ];
00896
00897 if ( NULL == _heightGrid ) return ( 1 );
00898
00899
00900
00901 index0 = 0;
00902
00903 startRow = _maxAoiRowIndex;
00904 endRow = _minAoiRowIndex;
00905
00906 index1 =
00907 BYTES_IN_HEADER +
00908 BYTES_PER_FLOAT *
00909 ( _minAoiColIndex +
00910 ( _nGridRows - startRow - 1 ) * _nGridCols );
00911
00912 for ( i = startRow; i >= endRow; i-- )
00913 {
00914
00915
00916 fin.seekg( index1 );
00917
00918 fin.read(
00919 (char*) ( &_heightGrid[ index0 ] ),
00920 ( _nAoiCols * BYTES_PER_FLOAT ));
00921
00922 if ( fin.fail() ) { fin.close(); return( 1 ); }
00923
00924
00925
00926 #if LITTLE_ENDIAN
00927
00928 this->swapBytes(
00929 &_heightGrid[ index0 ],
00930 BYTES_PER_FLOAT, _nAoiCols );
00931
00932 #endif
00933
00934 index0 += _nAoiCols;
00935 index1 += ( BYTES_PER_FLOAT * _nGridCols );
00936 }
00937
00938 fin.close();
00939
00940 }
00941
00942 catch ( ... ) { return( 1 ); }
00943
00944 return ( 0 );
00945
00946 }
00947
00948
00949
00950
00951
00952
00953
00954
00955 int
00956 Egm2008AoiGrid::loadGridMetadata( void )
00957 {
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973 try {
00974
00975 std::ifstream fin;
00976
00977
00978
00979 fin.open(
00980 _gridFname.c_str(),
00981 std::ifstream::binary | std::ifstream::in );
00982
00983 if ( fin.fail() ) { return( 1 ); }
00984
00985
00986
00987 fin.read((char*) &_nGridPad, BYTES_PER_INT );
00988
00989 if ( fin.fail() ) { fin.close(); return( 1 ); }
00990
00991 fin.read((char*) &_nOrigRows, BYTES_PER_INT );
00992
00993 if ( fin.fail() ) { fin.close(); return( 1 ); }
00994
00995 fin.read((char*) &_nOrigCols, BYTES_PER_INT );
00996
00997 if ( fin.fail() ) { fin.close(); return( 1 ); }
00998
00999 fin.read((char*) &_dLat, BYTES_PER_DOUBLE );
01000
01001 if ( fin.fail() ) { fin.close(); return( 1 ); }
01002
01003 fin.read((char*) &_dLon, BYTES_PER_DOUBLE );
01004
01005 if ( fin.fail() ) { fin.close(); return( 1 ); }
01006
01007 fin.close();
01008
01009 #if LITTLE_ENDIAN
01010
01011 this->swapBytes( &_nGridPad, BYTES_PER_INT, 1 );
01012 this->swapBytes( &_nOrigRows, BYTES_PER_INT, 1 );
01013 this->swapBytes( &_nOrigCols, BYTES_PER_INT, 1 );
01014 this->swapBytes( &_dLat, BYTES_PER_DOUBLE, 1 );
01015 this->swapBytes( &_dLon, BYTES_PER_DOUBLE, 1 );
01016
01017 #endif
01018
01019 _dLat *= RAD_PER_DEG;
01020 _dLon *= RAD_PER_DEG;
01021
01022
01023
01024 _nGridRows = _nOrigRows + ( 2 * _nGridPad );
01025 _nGridCols = _nOrigCols + ( 2 * _nGridPad ) + 1;
01026
01027 if ((( _dLat - EPSILON ) < DLAT1 ) &&
01028 (( _dLat + EPSILON ) > DLAT1 ) &&
01029 (( _dLon - EPSILON ) < DLON1 ) &&
01030 (( _dLon + EPSILON ) > DLON1 ))
01031 {
01032 _nomAoiCols = NOM_AOI_COLS1;
01033 _nomAoiRows = NOM_AOI_ROWS1;
01034 }
01035 else if ((( _dLat - EPSILON ) < DLAT2_5 ) &&
01036 (( _dLat + EPSILON ) > DLAT2_5 ) &&
01037 (( _dLon - EPSILON ) < DLON2_5 ) &&
01038 (( _dLon + EPSILON ) > DLON2_5 ))
01039 {
01040 _nomAoiCols = NOM_AOI_COLS2_5;
01041 _nomAoiRows = NOM_AOI_ROWS2_5;
01042 }
01043 else
01044 {
01045 ; return( 1 );
01046 }
01047
01048 _baseLatitude =
01049 -PIDIV2 - _dLat * double( _nGridPad );
01050 _baseLongitude =
01051 0.0 - _dLon * double( _nGridPad );
01052
01053 }
01054
01055 catch ( ... ) { return( 1 ); }
01056
01057 return( 0 );
01058
01059 }
01060
01062
01064