00001
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00026
00027
00028
00029
00030 #include <fstream>
00031
00032 #ifdef IRIXN32
00033 #include <math.h>
00034 #else
00035 #include <cmath>
00036 #endif
00037
00038 #include "CCSThreadLock.h"
00039 #include "CoordinateConversionException.h"
00040
00041 #include "egm2008_full_grid_package.h"
00042
00043 using namespace MSP;
00044
00045 namespace
00046 {
00047 const int BYTES_PER_DOUBLE = sizeof( double );
00048 const int BYTES_PER_FLOAT = sizeof( float );
00049 const int BYTES_PER_INT = sizeof( int );
00050
00051 const double EPSILON = 1.0e-15;
00052
00053 const double PI = 3.14159265358979323;
00054
00055 const double PIDIV2 = PI / 2.0;
00056 const double TWOPI = 2.0 * PI;
00057
00058 const double DEG_PER_RAD = 180.0 / PI;
00059 const double RAD_PER_DEG = PI / 180.0;
00060 }
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070 Egm2008FullGrid::Egm2008FullGrid (void)
00071 {
00072
00073
00074
00075
00076
00077
00078 int status;
00079
00080
00081
00082
00083
00084
00085
00086 _heightGrid = NULL;
00087
00088
00089
00090 status = this->loadGrid();
00091
00092 if ( status != 0 )
00093 {
00094 throw MSP::CCS::CoordinateConversionException(
00095 "Error: Egm2008GeoidGrid: constructor failed.");
00096 }
00097
00098 }
00099
00100
00101
00102
00103
00104
00105 Egm2008FullGrid::Egm2008FullGrid
00106 ( const Egm2008FullGrid& oldGrid )
00107
00108 : Egm2008GeoidGrid( oldGrid )
00109
00110 {
00111
00112
00113 int i;
00114 int kount;
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128 try
00129 {
00130 _heightGrid = NULL;
00131
00132 kount = _nGridRows * _nGridCols;
00133
00134 _heightGrid = new float[ kount ];
00135
00136 for ( i = 0; i < kount; i++ )
00137 {
00138 _heightGrid[ i ] = oldGrid._heightGrid[ i ];
00139 }
00140 }
00141 catch (...)
00142 {
00143 oldGrid._mutex.unlock();
00144
00145 throw MSP::CCS::CoordinateConversionException(
00146 "Error: Egm2008GeoidGrid: copy contructor failed");
00147 }
00148
00149 oldGrid._mutex.unlock();
00150
00151 }
00152
00153
00154
00155
00156
00157
00158 Egm2008FullGrid::~Egm2008FullGrid (void)
00159 {
00160
00161
00162
00163
00164
00165 delete[] _heightGrid;
00166
00167 }
00168
00169
00170
00171
00172
00173
00174 Egm2008FullGrid&
00175 Egm2008FullGrid::operator= ( const Egm2008FullGrid& oldGrid )
00176 {
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186 int i;
00187 int kount;
00188
00189
00190
00191
00192 if ( this == & oldGrid ) return ( *this );
00193
00194
00195
00196 Egm2008GeoidGrid::operator= ( oldGrid );
00197
00198
00199
00200 try
00201 {
00202 delete[] _heightGrid; _heightGrid = NULL;
00203
00204 kount = _nGridRows * _nGridCols;
00205
00206 _heightGrid = new float[ kount ];
00207
00208 for ( i = 0; i < kount; i++ )
00209 {
00210 _heightGrid[ i ] = oldGrid._heightGrid[ i ];
00211 }
00212 }
00213 catch (...)
00214 {
00215 _mutex.unlock();
00216 oldGrid._mutex.unlock();
00217
00218 throw MSP::CCS::CoordinateConversionException(
00219 "Error: Egm2008GeoidGrid: assignment operator failed");
00220 }
00221
00222 _mutex.unlock();
00223 oldGrid._mutex.unlock();
00224
00225 return( *this );
00226
00227 }
00228
00229
00230
00231
00232
00233
00234 int
00235 Egm2008FullGrid::geoidHeight(
00236 int wSize,
00237 double latitude,
00238 double longitude,
00239 double& gHeight )
00240 {
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289 try {
00290
00291
00292
00293 const int TWENTY = 20;
00294
00295 bool oddSize;
00296
00297 int i;
00298 int i0;
00299 int iIndex;
00300 int iMin;
00301 int j;
00302 int j0;
00303 int jIndex;
00304 int jMin;
00305 int kIndex;
00306 int offset;
00307 int status;
00308
00309 double latIndex;
00310 double lonIndex;
00311 double temp;
00312
00313 double latSupport[ TWENTY ];
00314 double lonSupport[ TWENTY ];
00315 double moments [ TWENTY ];
00316
00317
00318
00319 gHeight = 0.0;
00320
00321 if ( TWENTY != MAX_WSIZE ) return( 1 );
00322
00323 if ( wSize > MAX_WSIZE ) wSize = MAX_WSIZE;
00324
00325 if (( latitude < -PIDIV2 ) || ( latitude > PIDIV2 )) return( 1 );
00326
00327
00328
00329 while ( longitude < 0.0 ) longitude += TWOPI;
00330 while ( longitude > TWOPI ) longitude -= TWOPI;
00331
00332 for (i = 0; i < TWENTY; i++)
00333 {
00334 latSupport[ i ] = lonSupport[ i ] = moments[ i ] = 0.0;
00335 }
00336
00337
00338
00339
00340
00341
00342
00343 if ( wSize < 3 )
00344 {
00345 status =
00346 this->geoidHeight( latitude, longitude, gHeight );
00347
00348 if ( status != 0 ) return( 1 );
00349
00350 ; return( 0 );
00351 }
00352
00353
00354
00355
00356
00357
00358 latIndex =
00359 double( _nGridPad ) + ( latitude + PIDIV2 ) / _dLat;
00360 lonIndex =
00361 double( _nGridPad ) + ( longitude - 0.0 ) / _dLon;
00362
00363 oddSize = ( wSize != (( wSize / 2 ) * 2 ));
00364
00365 if ( oddSize ) {
00366 i0 = int( latIndex + 0.5 );
00367 j0 = int( lonIndex + 0.5 );
00368
00369 iMin = i0 - ( wSize / 2 );
00370 jMin = j0 - ( wSize / 2 );
00371 }
00372 else {
00373 i0 = int( latIndex );
00374 j0 = int( lonIndex );
00375
00376 iMin = i0 - ( wSize / 2 ) + 1;
00377 jMin = j0 - ( wSize / 2 ) + 1;
00378 }
00379
00380
00381
00382
00383
00384
00385
00386 temp =
00387 lonIndex - double( jMin );
00388
00389
00390
00391
00392 for ( i = 0; i < wSize; i++ )
00393 {
00394 iIndex = iMin + i;
00395
00396 offset =
00397 ( _nGridRows - iIndex - 1 ) * _nGridCols;
00398
00399
00400
00401
00402 for ( j = 0; j < wSize; j++ )
00403 {
00404 jIndex = jMin + j;
00405
00406 kIndex = jIndex + offset;
00407
00408 lonSupport[j] = _heightGrid[ kIndex ];
00409 }
00410
00411
00412
00413
00414
00415
00416 status =
00417 this->initSpline( wSize, lonSupport, moments );
00418
00419 if ( status != 0 ) return( 1 );
00420
00421 latSupport[i] =
00422 this->spline( wSize, temp, lonSupport, moments );
00423 }
00424
00425
00426
00427
00428
00429 temp =
00430 latIndex - double( iMin );
00431
00432
00433
00434
00435 status =
00436 this->initSpline( wSize, latSupport, moments );
00437
00438 if ( status != 0 ) return( 1 );
00439
00440 gHeight =
00441 this->spline( wSize, temp, latSupport, moments );
00442
00443
00444
00445
00446 }
00447
00448 catch ( ... ) { gHeight = 0.0; return( 1 ); }
00449
00450 return( 0 );
00451
00452 }
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463 int
00464 Egm2008FullGrid::geoidHeight(
00465 double latitude,
00466 double longitude,
00467 double& gHeight )
00468 {
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490 try {
00491
00492 int i1;
00493 int i2;
00494 int i3;
00495 int i4;
00496 int index;
00497 int j1;
00498 int j2;
00499 int j3;
00500 int j4;
00501 int status;
00502
00503 double a0;
00504 double a1;
00505 double a2;
00506 double a3;
00507 double lat1;
00508 double lon1;
00509 double n1;
00510 double n2;
00511 double n3;
00512 double n4;
00513 double x;
00514 double y;
00515
00516
00517
00518 gHeight = 0.0;
00519
00520 if (( latitude < -PIDIV2 ) || ( latitude > PIDIV2 )) return( 1 );
00521
00522 while ( longitude < 0.0 ) longitude += TWOPI;
00523 while ( longitude > TWOPI ) longitude -= TWOPI;
00524
00525
00526
00527 status =
00528 this->swGridIndices
00529 ( latitude, longitude, i1, j1 );
00530
00531 if ( status != 0 ) return( 1 );
00532
00533 i2 = i1;
00534 j2 = j1 + 1;
00535
00536 i3 = i1 + 1;
00537 j3 = j2;
00538
00539 i4 = i3;
00540 j4 = j1;
00541
00542
00543
00544 index =
00545 j1 + ( _nGridRows - i1 - 1 ) * _nGridCols;
00546 n1 = _heightGrid[ index ];
00547
00548 index =
00549 j2 + ( _nGridRows - i2 - 1 ) * _nGridCols;
00550 n2 = _heightGrid[ index ];
00551
00552 index =
00553 j3 + ( _nGridRows - i3 - 1 ) * _nGridCols;
00554 n3 = _heightGrid[ index ];
00555
00556 index =
00557 j4 + ( _nGridRows - i4 - 1 ) * _nGridCols;
00558 n4 = _heightGrid[ index ];
00559
00560
00561
00562
00563
00564 a0 = n1;
00565 a1 = n2 - n1;
00566 a2 = n4 - n1;
00567 a3 = n1 + n3 - n2 - n4;
00568
00569
00570
00571 lat1 =
00572 _baseLatitude + _dLat * double( i1 );
00573 lon1 =
00574 _baseLongitude + _dLon * double( j1 );
00575
00576
00577
00578
00579
00580
00581 x = ( longitude - lon1 ) / _dLon;
00582 y = ( latitude - lat1 ) / _dLat;
00583
00584
00585
00586 gHeight = a0 + a1 * x + a2 * y + a3 * x * y;
00587
00588 }
00589
00590 catch ( ... ) { gHeight = 0.0; return( 1 ); }
00591
00592 return( 0 );
00593
00594 }
00595
00596
00597
00598
00599
00600
00601
00602 int
00603 Egm2008FullGrid::loadGrid( void )
00604 {
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620 try {
00621
00622 int endRow;
00623 int i;
00624 int index0;
00625 int index1;
00626 int index2;
00627 int j;
00628 int k;
00629 int kount;
00630 int startRow;
00631
00632 std::ifstream fin;
00633
00634
00635
00636 fin.open(
00637 _gridFname.c_str(),
00638 std::ifstream::binary | std::ifstream::in );
00639
00640 if ( fin.fail() ) return( 1 );
00641
00642
00643
00644 fin.read((char*) &_nGridPad, BYTES_PER_INT );
00645
00646 if (fin.fail()) { fin.close(); return( 1 ); }
00647
00648 fin.read((char*) &_nOrigRows, BYTES_PER_INT );
00649
00650 if (fin.fail()) { fin.close(); return( 1 ); }
00651
00652 fin.read((char*) &_nOrigCols, BYTES_PER_INT );
00653
00654 if (fin.fail()) { fin.close(); return( 1 ); }
00655
00656 fin.read((char*) &_dLat, BYTES_PER_DOUBLE );
00657
00658 if (fin.fail()) { fin.close(); return( 1 ); }
00659
00660 fin.read((char*) &_dLon, BYTES_PER_DOUBLE );
00661
00662 if (fin.fail()) { fin.close(); return( 1 ); }
00663
00664
00665
00666 #if LITTLE_ENDIAN
00667
00668 this->swapBytes( &_nGridPad, BYTES_PER_INT, 1 );
00669 this->swapBytes( &_nOrigRows, BYTES_PER_INT, 1 );
00670 this->swapBytes( &_nOrigCols, BYTES_PER_INT, 1 );
00671 this->swapBytes( &_dLat, BYTES_PER_DOUBLE, 1 );
00672 this->swapBytes( &_dLon, BYTES_PER_DOUBLE, 1 );
00673
00674 #endif
00675
00676 _dLat *= RAD_PER_DEG;
00677 _dLon *= RAD_PER_DEG;
00678
00679
00680
00681 _nGridRows = _nOrigRows + ( 2 * _nGridPad );
00682 _nGridCols = _nOrigCols + ( 2 * _nGridPad ) + 1;
00683
00684 _baseLatitude =
00685 -PIDIV2 - _dLat * double( _nGridPad );
00686 _baseLongitude =
00687 0.0 - _dLon * double( _nGridPad );
00688
00689
00690
00691
00692
00693
00694
00695 _heightGrid = NULL;
00696
00697 _heightGrid =
00698 new float[ _nGridRows * _nGridCols ];
00699
00700 if ( NULL == _heightGrid ) return( 1 );
00701
00702 index0 = 0;
00703
00704 startRow = _nGridRows - 1;
00705 endRow = 0;
00706
00707 for ( i = startRow; i >= endRow; i-- )
00708 {
00709
00710
00711 fin.read(
00712 (char*) (&_heightGrid[ index0 ]),
00713 ( _nGridCols * BYTES_PER_FLOAT ));
00714
00715 if ( fin.fail() ) { fin.close(); return( 1 ); }
00716
00717
00718
00719 #if LITTLE_ENDIAN
00720
00721 this->swapBytes(
00722 &_heightGrid[ index0 ], BYTES_PER_FLOAT, _nGridCols );
00723
00724 #endif
00725
00726 index0 += _nGridCols;
00727 }
00728
00729 fin.close();
00730
00731 }
00732
00733 catch ( ... ) { return( 1 ); }
00734
00735 return ( 0 );
00736
00737 }
00738
00740
00742