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 #include <string.h>
00081 #include <stdlib.h>
00082 #include <stdio.h>
00083 #include "GeoidLibrary.h"
00084 #include "CoordinateConversionException.h"
00085 #include "ErrorMessages.h"
00086 #include "CCSThreadMutex.h"
00087 #include "CCSThreadLock.h"
00088
00089 #include <vector>
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102 using namespace MSP::CCS;
00103 using MSP::CCSThreadMutex;
00104 using MSP::CCSThreadLock;
00105
00106
00107
00108
00109
00110
00111
00112 const double PI = 3.14159265358979323e0;
00113 const double PI_OVER_2 = PI / 2.0;
00114 const double TWO_PI = 2.0 * PI;
00115 const double _180_OVER_PI = (180.0 / PI);
00116 const int EGM96_COLS = 1441;
00117 const int EGM96_ROWS = 721;
00118 const int EGM84_COLS = 37;
00119 const int EGM84_ROWS = 19;
00120 const int EGM84_30_MIN_COLS = 721;
00121 const int EGM84_30_MIN_ROWS = 361;
00122 const int EGM96_HEADER_ITEMS = 6;
00123 const double SCALE_FACTOR_15_MINUTES = .25;
00124 const double SCALE_FACTOR_10_DEGREES = 10;
00125 const double SCALE_FACTOR_30_MINUTES = .5;
00126 const double SCALE_FACTOR_1_DEGREE = 1;
00127 const double SCALE_FACTOR_2_DEGREES = 2;
00128 const int EGM96_ELEVATIONS = EGM96_COLS * EGM96_ROWS;
00129 const int EGM84_ELEVATIONS = EGM84_COLS * EGM84_ROWS;
00130 const int EGM84_30_MIN_ELEVATIONS = EGM84_30_MIN_COLS * EGM84_30_MIN_ROWS;
00131 const int EGM96_INSET_AREAS = 53;
00132
00133
00134
00135 struct EGM96_Variable_Grid
00136 {
00137 double min_lat;
00138 double max_lat;
00139 double min_lon;
00140 double max_lon;
00141 };
00142
00143
00144 const EGM96_Variable_Grid EGM96_Variable_Grid_Table[EGM96_INSET_AREAS] =
00145 {{74.5, 75.5, 273.5, 280.0},
00146 {66.5, 67.5, 293.5, 295.0},
00147 {62.5, 64.0, 133.0, 136.5},
00148 {60.5, 61.5, 208.5, 210.0},
00149 {60.5, 61.0, 219.0, 220.5},
00150 {51.0, 53.0, 172.0, 174.5},
00151 {52.0, 53.0, 192.5, 194.0},
00152 {51.0, 52.0, 188.5, 191.0},
00153 {50.0, 52.0, 178.0, 182.5},
00154 {43.0, 46.0, 148.0, 153.5},
00155 {43.0, 45.0, 84.0, 89.5},
00156 {40.0, 41.0, 70.5, 72.0},
00157 {36.5, 37.0, 78.5, 79.0},
00158 {36.0, 37.0, 348.0, 349.5},
00159 {35.0, 36.0, 171.0, 172.5},
00160 {34.0, 35.0, 140.5, 142.0},
00161 {29.5, 31.0, 78.5, 81.0},
00162 {28.5, 30.0, 81.5, 83.0},
00163 {26.5, 30.0, 142.0, 143.5},
00164 {26.0, 29.0, 91.5, 96.0},
00165 {27.5, 29.0, 84.0, 86.5},
00166 {28.0, 29.0, 342.5, 344.0},
00167 {26.5, 28.0, 88.5, 90.0},
00168 {25.0, 26.0, 189.0, 190.5},
00169 {23.0, 24.0, 195.0, 196.5},
00170 {21.0, 21.5, 204.0, 204.5},
00171 {20.0, 21.0, 283.5, 288.0},
00172 {18.5, 20.5, 204.0, 205.5},
00173 {18.0, 20.0, 291.0, 296.5},
00174 {17.0, 18.0, 298.0, 299.5},
00175 {15.0, 16.0, 122.0, 123.5},
00176 {12.0, 14.0, 144.5, 147.0},
00177 {11.0, 12.0, 141.5, 144.0},
00178 {9.5, 11.5, 125.0, 127.5},
00179 {10.0, 11.0, 286.0, 287.5},
00180 {6.0, 9.5, 287.0, 289.5},
00181 {5.0, 7.0, 124.0, 128.5},
00182 {-1.0, 1.0, 125.0, 128.5},
00183 {-3.0, -1.5, 281.0, 282.5},
00184 {-7.0, -5.0, 150.5, 155.0},
00185 {-8.0, -7.0, 107.0, 108.5},
00186 {-9.0, -7.0, 147.0, 149.5},
00187 {-11.0, -10.0, 161.5, 163.0},
00188 {-14.5, -13.5, 166.0, 167.5},
00189 {-18.5, -17.0, 186.5, 188.0},
00190 {-20.5, -20.0, 168.0, 169.5},
00191 {-23.0, -20.0, 184.5, 187.0},
00192 {-27.0, -24.0, 288.0, 290.5},
00193 {-53.0, -52.0, 312.0, 313.5},
00194 {-56.0, -55.0, 333.0, 334.5},
00195 {-61.5, -60.0, 312.5, 317.0},
00196 {-61.5, -60.5, 300.5, 303.0},
00197 {-73.0, -72.0, 24.5, 26.0}};
00198
00199
00200
00201
00202
00203
00204
00205 void swapBytes(
00206 void *buffer,
00207 size_t size,
00208 size_t count)
00209 {
00210 char *b = (char *) buffer;
00211 char temp;
00212 for (size_t c = 0; c < count; c ++)
00213 {
00214 for (size_t s = 0; s < size / 2; s ++)
00215 {
00216 temp = b[c*size + s];
00217 b[c*size + s] = b[c*size + size - s - 1];
00218 b[c*size + size - s - 1] = temp;
00219 }
00220 }
00221 }
00222
00223
00224 size_t readBinary(
00225 void *buffer,
00226 size_t size,
00227 size_t count,
00228 FILE *stream )
00229 {
00230 size_t actual_count = fread(buffer, size, count, stream);
00231
00232 if(ferror(stream) || actual_count < count )
00233 {
00234 char message[256] = "";
00235 strcpy( message, "Error reading binary file." );
00236 throw CoordinateConversionException( message );
00237 }
00238
00239 #ifdef LITTLE_ENDIAN
00240 swapBytes(buffer, size, count);
00241 #endif
00242
00243 return actual_count;
00244 }
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254 class MSP::CCS::GeoidLibraryCleaner
00255 {
00256 public:
00257
00258 ~GeoidLibraryCleaner()
00259 {
00260 CCSThreadLock lock(&GeoidLibrary::mutex);
00261 GeoidLibrary::deleteInstance();
00262 }
00263
00264 } geoidLibraryCleanerInstance;
00265
00266
00267
00268 CCSThreadMutex GeoidLibrary::mutex;
00269 GeoidLibrary* GeoidLibrary::instance = 0;
00270 int GeoidLibrary::instanceCount = 0;
00271
00272
00273 GeoidLibrary* GeoidLibrary::getInstance()
00274 {
00275 CCSThreadLock lock(&mutex);
00276 if( instance == 0 )
00277 instance = new GeoidLibrary;
00278
00279 instanceCount++;
00280
00281 return instance;
00282 }
00283
00284
00285 void GeoidLibrary::removeInstance()
00286 {
00287
00288
00289
00290
00291 CCSThreadLock lock(&mutex);
00292 if ( --instanceCount < 1 )
00293 {
00294 deleteInstance();
00295 }
00296 }
00297
00298
00299 void GeoidLibrary::deleteInstance()
00300 {
00301
00302
00303
00304
00305 if( instance != 0 )
00306 {
00307 delete instance;
00308 instance = 0;
00309 }
00310 }
00311
00312
00313 GeoidLibrary::GeoidLibrary()
00314 {
00315 loadGeoids();
00316 }
00317
00318
00319 GeoidLibrary::GeoidLibrary( const GeoidLibrary &gl )
00320 {
00321 *this = gl;
00322 }
00323
00324
00325 GeoidLibrary::~GeoidLibrary()
00326 {
00327 delete [] egm96GeoidList;
00328 delete [] egm84GeoidList;
00329 delete [] egm84ThirtyMinGeoidList;
00330 }
00331
00332
00333 GeoidLibrary& GeoidLibrary::operator=( const GeoidLibrary &gl )
00334 {
00335 if ( &gl == this )
00336 return *this;
00337
00338 egm96GeoidList = new float[EGM96_ELEVATIONS];
00339 egm84GeoidList = new float[EGM84_ELEVATIONS];
00340 egm84ThirtyMinGeoidList = new double[EGM84_30_MIN_ELEVATIONS];
00341
00342 for( int i = 0; i < EGM96_ELEVATIONS; i++ )
00343 {
00344 egm96GeoidList[i] = gl.egm96GeoidList[i];
00345 }
00346
00347 for( int j = 0; j < EGM84_ELEVATIONS; j++ )
00348 {
00349 egm84GeoidList[j] = gl.egm84GeoidList[j];
00350 }
00351
00352 for( int k = 0; k < EGM84_30_MIN_ELEVATIONS; k++ )
00353 {
00354 egm84ThirtyMinGeoidList[k] = gl.egm84ThirtyMinGeoidList[k];
00355 }
00356
00357 return *this;
00358 }
00359
00360
00361 void GeoidLibrary::loadGeoids()
00362 {
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373 initializeEGM96Geoid();
00374 initializeEGM84Geoid();
00375 initializeEGM84ThirtyMinGeoid();
00376 }
00377
00378
00379 void GeoidLibrary::convertEllipsoidToEGM96FifteenMinBilinearGeoidHeight(
00380 double longitude,
00381 double latitude,
00382 double ellipsoidHeight,
00383 double *geoidHeight )
00384 {
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397 double delta_height;
00398
00399 bilinearInterpolate(
00400 longitude, latitude,
00401 SCALE_FACTOR_15_MINUTES, EGM96_COLS, EGM96_ROWS, egm96GeoidList,
00402 &delta_height );
00403
00404 *geoidHeight = ellipsoidHeight - delta_height;
00405 }
00406
00407
00408 void GeoidLibrary::convertEllipsoidToEGM96VariableNaturalSplineHeight(
00409 double longitude,
00410 double latitude,
00411 double ellipsoidHeight,
00412 double *geoidHeight )
00413 {
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426 int i = 0;
00427 int num_cols = EGM96_COLS;
00428 int num_rows = EGM96_ROWS;
00429 double latitude_degrees = latitude * _180_OVER_PI;
00430 double longitude_degrees = longitude * _180_OVER_PI;
00431 double scale_factor = SCALE_FACTOR_15_MINUTES;
00432 double delta_height;
00433 long found = 0;
00434
00435 if( longitude_degrees < 0.0 )
00436 longitude_degrees += 360.0;
00437
00438 while( !found && i < EGM96_INSET_AREAS )
00439 {
00440 if(( latitude_degrees >= EGM96_Variable_Grid_Table[i].min_lat ) &&
00441 ( longitude_degrees >= EGM96_Variable_Grid_Table[i].min_lon ) &&
00442 ( latitude_degrees < EGM96_Variable_Grid_Table[i].max_lat ) &&
00443 ( longitude_degrees < EGM96_Variable_Grid_Table[i].max_lon ) )
00444 {
00445 scale_factor = SCALE_FACTOR_30_MINUTES;
00446 num_cols = 721;
00447 num_rows = 361;
00448 found = 1;
00449 }
00450
00451 i++;
00452 }
00453
00454 if( !found )
00455 {
00456 if( latitude_degrees >= -60.0 && latitude_degrees < 60.0 )
00457 {
00458 scale_factor = SCALE_FACTOR_1_DEGREE;
00459 num_cols = 361;
00460 num_rows = 181;
00461 }
00462 else
00463 {
00464 scale_factor = SCALE_FACTOR_2_DEGREES;
00465 num_cols = 181;
00466 num_rows = 91;
00467 }
00468 }
00469
00470 naturalSplineInterpolate(
00471 longitude, latitude,
00472 scale_factor, num_cols, num_rows, EGM96_ELEVATIONS-1, egm96GeoidList,
00473 &delta_height );
00474
00475 *geoidHeight = ellipsoidHeight - delta_height;
00476 }
00477
00478
00479 void GeoidLibrary::convertEllipsoidToEGM84TenDegBilinearHeight(
00480 double longitude,
00481 double latitude,
00482 double ellipsoidHeight,
00483 double *geoidHeight )
00484 {
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497 double delta_height;
00498
00499 bilinearInterpolate(
00500 longitude, latitude,
00501 SCALE_FACTOR_10_DEGREES, EGM84_COLS, EGM84_ROWS, egm84GeoidList,
00502 &delta_height );
00503
00504 *geoidHeight = ellipsoidHeight - delta_height;
00505 }
00506
00507
00508 void GeoidLibrary::convertEllipsoidToEGM84TenDegNaturalSplineHeight(
00509 double longitude,
00510 double latitude,
00511 double ellipsoidHeight,
00512 double *geoidHeight )
00513 {
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526 double delta_height;
00527
00528 naturalSplineInterpolate(
00529 longitude, latitude,
00530 SCALE_FACTOR_10_DEGREES, EGM84_COLS, EGM84_ROWS, EGM84_ELEVATIONS-1,
00531 egm84GeoidList, &delta_height );
00532
00533 *geoidHeight = ellipsoidHeight - delta_height;
00534 }
00535
00536
00537 void GeoidLibrary::convertEllipsoidToEGM84ThirtyMinBiLinearHeight(
00538 double longitude, double latitude, double ellipsoidHeight,
00539 double *geoidHeight )
00540 {
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554 double delta_height;
00555
00556 bilinearInterpolateDoubleHeights(
00557 longitude, latitude,
00558 SCALE_FACTOR_30_MINUTES, EGM84_30_MIN_COLS, EGM84_30_MIN_ROWS,
00559 egm84ThirtyMinGeoidList, &delta_height );
00560
00561 *geoidHeight = ellipsoidHeight - delta_height;
00562 }
00563
00564
00565 void GeoidLibrary::convertEGM96FifteenMinBilinearGeoidToEllipsoidHeight( double longitude, double latitude, double geoidHeight, double *ellipsoidHeight )
00566 {
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579 double delta_height;
00580
00581 bilinearInterpolate(
00582 longitude, latitude,
00583 SCALE_FACTOR_15_MINUTES, EGM96_COLS, EGM96_ROWS, egm96GeoidList,
00584 &delta_height );
00585
00586 *ellipsoidHeight = geoidHeight + delta_height;
00587 }
00588
00589
00590 void GeoidLibrary::convertEGM96VariableNaturalSplineToEllipsoidHeight(
00591 double longitude,
00592 double latitude,
00593 double geoidHeight,
00594 double *ellipsoidHeight )
00595 {
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608 int i = 0;
00609 int num_cols = EGM96_COLS;
00610 int num_rows = EGM96_ROWS;
00611 double latitude_degrees = latitude * _180_OVER_PI;
00612 double longitude_degrees = longitude * _180_OVER_PI;
00613 double scale_factor = SCALE_FACTOR_15_MINUTES;
00614 double delta_height;
00615 long found = 0;
00616
00617 if( longitude_degrees < 0.0 )
00618 longitude_degrees += 360.0;
00619
00620 while( !found && i < EGM96_INSET_AREAS )
00621 {
00622 if(( latitude_degrees >= EGM96_Variable_Grid_Table[i].min_lat ) &&
00623 ( longitude_degrees >= EGM96_Variable_Grid_Table[i].min_lon ) &&
00624 ( latitude_degrees < EGM96_Variable_Grid_Table[i].max_lat ) &&
00625 ( longitude_degrees < EGM96_Variable_Grid_Table[i].max_lon ) )
00626 {
00627 scale_factor = SCALE_FACTOR_30_MINUTES;
00628 num_cols = 721;
00629 num_rows = 361;
00630 found = 1;
00631 }
00632
00633 i++;
00634 }
00635
00636 if( !found )
00637 {
00638 if( latitude_degrees >= -60.0 && latitude_degrees < 60.0 )
00639 {
00640 scale_factor = SCALE_FACTOR_1_DEGREE;
00641 num_cols = 361;
00642 num_rows = 181;
00643 }
00644 else
00645 {
00646 scale_factor = SCALE_FACTOR_2_DEGREES;
00647 num_cols = 181;
00648 num_rows = 91;
00649 }
00650 }
00651
00652 naturalSplineInterpolate(
00653 longitude, latitude,
00654 scale_factor, num_cols, num_rows, EGM96_ELEVATIONS-1, egm96GeoidList,
00655 &delta_height );
00656
00657 *ellipsoidHeight = geoidHeight + delta_height;
00658 }
00659
00660
00661 void GeoidLibrary::convertEGM84TenDegBilinearToEllipsoidHeight( double longitude, double latitude, double geoidHeight, double *ellipsoidHeight )
00662 {
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675 double delta_height;
00676
00677 bilinearInterpolate(
00678 longitude, latitude,
00679 SCALE_FACTOR_10_DEGREES, EGM84_COLS, EGM84_ROWS, egm84GeoidList,
00680 &delta_height );
00681
00682 *ellipsoidHeight = geoidHeight + delta_height;
00683 }
00684
00685
00686 void GeoidLibrary::convertEGM84TenDegNaturalSplineToEllipsoidHeight(
00687 double longitude,
00688 double latitude,
00689 double geoidHeight,
00690 double *ellipsoidHeight )
00691 {
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704 double delta_height;
00705
00706 naturalSplineInterpolate(
00707 longitude, latitude, SCALE_FACTOR_10_DEGREES,
00708 EGM84_COLS, EGM84_ROWS, EGM84_ELEVATIONS-1,
00709 egm84GeoidList, &delta_height );
00710
00711 *ellipsoidHeight = geoidHeight + delta_height;
00712 }
00713
00714
00715 void GeoidLibrary::convertEGM84ThirtyMinBiLinearToEllipsoidHeight(
00716 double longitude, double latitude, double geoidHeight,
00717 double *ellipsoidHeight )
00718 {
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732 double delta_height;
00733
00734 bilinearInterpolateDoubleHeights( longitude, latitude, SCALE_FACTOR_30_MINUTES,
00735 EGM84_30_MIN_COLS, EGM84_30_MIN_ROWS, egm84ThirtyMinGeoidList,
00736 &delta_height );
00737 *ellipsoidHeight = geoidHeight + delta_height;
00738 }
00739
00740
00741
00742
00743
00744
00745
00746 void GeoidLibrary::initializeEGM96Geoid()
00747 {
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757 int items_read = 0;
00758 char* file_name = 0;
00759 char* path_name = getenv( "MSPCCS_DATA" );
00760 long elevations_read = 0;
00761 long items_discarded = 0;
00762 long num = 0;
00763 FILE* geoid_height_file;
00764
00765 CCSThreadLock lock(&mutex);
00766
00767
00768
00769
00770 if (path_name != NULL)
00771 {
00772 file_name = new char[ strlen( path_name ) + 11 ];
00773 strcpy( file_name, path_name );
00774 strcat( file_name, "/" );
00775 }
00776 else
00777 {
00778 file_name = new char[ 21 ];
00779 strcpy( file_name, "../../data/" );
00780 }
00781 strcat( file_name, "egm96.grd" );
00782
00783
00784
00785 if ( ( geoid_height_file = fopen( file_name, "rb" ) ) == NULL )
00786 {
00787 delete [] file_name;
00788 file_name = 0;
00789
00790 char message[256] = "";
00791 strcpy( message, ErrorMessages::geoidFileOpenError );
00792 strcat( message, ": egm96.grd\n" );
00793 throw CoordinateConversionException( message );
00794 }
00795
00796
00797
00798 float egm96GeoidHeaderList[EGM96_HEADER_ITEMS];
00799 items_discarded = readBinary(
00800 egm96GeoidHeaderList, 4, EGM96_HEADER_ITEMS, geoid_height_file );
00801
00802
00803
00804 if( egm96GeoidHeaderList[0] != -90.0 ||
00805 egm96GeoidHeaderList[1] != 90.0 ||
00806 egm96GeoidHeaderList[2] != 0.0 ||
00807 egm96GeoidHeaderList[3] != 360.0 ||
00808 egm96GeoidHeaderList[4] != SCALE_FACTOR_15_MINUTES ||
00809 egm96GeoidHeaderList[5] != SCALE_FACTOR_15_MINUTES ||
00810 items_discarded != EGM96_HEADER_ITEMS )
00811 {
00812 fclose( geoid_height_file );
00813 delete [] file_name;
00814 file_name = 0;
00815
00816 char message[256] = "";
00817 strcpy( message, ErrorMessages::geoidFileParseError );
00818 strcat( message, ": egm96.grd\n" );
00819 throw CoordinateConversionException( message );
00820 }
00821
00822
00823 egm96GeoidList = new float[EGM96_ELEVATIONS];
00824 elevations_read = readBinary(
00825 egm96GeoidList, 4, EGM96_ELEVATIONS, geoid_height_file );
00826
00827 fclose( geoid_height_file );
00828 geoid_height_file = 0;
00829
00830 delete [] file_name;
00831 file_name = 0;
00832 }
00833
00834
00835 void GeoidLibrary::initializeEGM84Geoid()
00836 {
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846 int items_read = 0;
00847 char* file_name = 0;
00848 char* path_name = getenv( "MSPCCS_DATA" );
00849 long elevations_read = 0;
00850 long num = 0;
00851 FILE* geoid_height_file;
00852
00853 CCSThreadLock lock(&mutex);
00854
00855
00856
00857
00858 if (path_name != NULL)
00859 {
00860 file_name = new char[ strlen( path_name ) + 11 ];
00861 strcpy( file_name, path_name );
00862 strcat( file_name, "/" );
00863 }
00864 else
00865 {
00866 file_name =new char[ 21 ];
00867 strcpy( file_name, "../../data/" );
00868 }
00869 strcat( file_name, "egm84.grd" );
00870
00871
00872
00873 if( ( geoid_height_file = fopen( file_name, "rb" ) ) == NULL )
00874 {
00875 delete [] file_name;
00876 file_name = 0;
00877
00878 char message[256] = "";
00879 strcpy( message, ErrorMessages::geoidFileOpenError );
00880 strcat( message, ": egm84.grd\n" );
00881 throw CoordinateConversionException( message );
00882 }
00883
00884
00885
00886 egm84GeoidList = new float[EGM84_ELEVATIONS];
00887 elevations_read = readBinary(
00888 egm84GeoidList, 4, EGM84_ELEVATIONS, geoid_height_file );
00889
00890 fclose( geoid_height_file );
00891
00892 delete [] file_name;
00893 file_name = 0;
00894 }
00895
00896 void GeoidLibrary::initializeEGM84ThirtyMinGeoid()
00897 {
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907 int items_read = 0;
00908 char* file_name = 0;
00909 char* path_name = getenv( "MSPCCS_DATA" );
00910 long elevations_read = 0;
00911 long num = 0;
00912 FILE* geoid_height_file;
00913
00914 CCSThreadLock lock(&mutex);
00915
00916
00917
00918
00919 if (path_name != NULL)
00920 {
00921 file_name = new char[ strlen( path_name ) + 12 ];
00922 strcpy( file_name, path_name );
00923 strcat( file_name, "/" );
00924 }
00925 else
00926 {
00927 file_name =new char[ 22 ];
00928 strcpy( file_name, "../../data/" );
00929 }
00930 strcat( file_name, "wwgrid.bin" );
00931
00932
00933
00934 if( ( geoid_height_file = fopen( file_name, "rb" ) ) == NULL )
00935 {
00936 delete [] file_name;
00937 file_name = 0;
00938
00939 char message[256] = "";
00940 strcpy( message, ErrorMessages::geoidFileOpenError );
00941 strcat( message, ": wwgrid.bin\n" );
00942 throw CoordinateConversionException( message );
00943 }
00944
00945
00946
00947
00948 egm84ThirtyMinGeoidList = new double[EGM84_30_MIN_ELEVATIONS];
00949 elevations_read = readBinary(
00950 egm84ThirtyMinGeoidList, 8, EGM84_30_MIN_ELEVATIONS, geoid_height_file );
00951
00952 fclose( geoid_height_file );
00953
00954 delete [] file_name;
00955 file_name = 0;
00956 }
00957
00958 void GeoidLibrary::bilinearInterpolateDoubleHeights(
00959 double longitude,
00960 double latitude,
00961 double scale_factor,
00962 int num_cols,
00963 int num_rows,
00964 double *height_buffer,
00965 double *delta_height )
00966 {
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983 int index;
00984 int post_x, post_y;
00985 double offset_x, offset_y;
00986 double delta_x, delta_y;
00987 double _1_minus_delta_x, _1_minus_delta_y;
00988 double latitude_dd, longitude_dd;
00989 double height_se, height_sw, height_ne, height_nw;
00990 double w_sw, w_se, w_ne, w_nw;
00991 double south_lat, west_lon;
00992 int end_index = 0;
00993 int max_index = num_rows * num_cols - 1;
00994 char errorStatus[50] = "";
00995
00996 if( ( latitude < -PI_OVER_2 ) || ( latitude > PI_OVER_2 ) )
00997 {
00998 strcat( errorStatus, ErrorMessages::latitude );
00999 }
01000 if( ( longitude < -PI ) || ( longitude > TWO_PI ) )
01001 {
01002 strcat( errorStatus, ErrorMessages::longitude );
01003 }
01004
01005 if( strlen( errorStatus ) > 0)
01006 throw CoordinateConversionException( errorStatus );
01007
01008 latitude_dd = latitude * _180_OVER_PI;
01009 longitude_dd = longitude * _180_OVER_PI;
01010
01011
01012
01013 if( longitude_dd < 0.0 )
01014 {
01015 offset_x = ( longitude_dd + 360.0 ) / scale_factor;
01016 }
01017 else
01018 {
01019 offset_x = longitude_dd / scale_factor;
01020 }
01021 offset_y = ( 90 - latitude_dd ) / scale_factor;
01022
01023
01024
01025
01026 post_x = ( int )( offset_x );
01027 if( ( post_x + 1 ) == num_cols )
01028 post_x--;
01029 post_y = ( int )( offset_y + 1.0e-11 );
01030 if( ( post_y + 1 ) == num_rows )
01031 post_y--;
01032
01033
01034 index = post_y * num_cols + post_x;
01035 if( index < 0 )
01036 height_nw = height_buffer[ 0 ];
01037 else if( index > max_index )
01038 height_nw = height_buffer[ max_index ];
01039 else
01040 height_nw = height_buffer[ index ];
01041
01042 end_index = index + 1;
01043 if( end_index > max_index )
01044 height_ne = height_buffer[ max_index ];
01045 else
01046 height_ne = height_buffer[ end_index ];
01047
01048
01049 index = ( post_y + 1 ) * num_cols + post_x;
01050 if( index < 0 )
01051 height_sw = height_buffer[ 0 ];
01052 else if( index > max_index )
01053 height_sw = height_buffer[ max_index ];
01054 else
01055 height_sw = height_buffer[ index ];
01056
01057
01058 end_index = index + 1;
01059 if( end_index > max_index )
01060 height_se = height_buffer[ max_index ];
01061 else
01062 height_se = height_buffer[ end_index ];
01063
01064 west_lon = post_x * scale_factor;
01065
01066
01067 south_lat = ( 90 - ( post_y * scale_factor ) ) - scale_factor;
01068
01069
01070
01071 if( longitude_dd < 0.0 )
01072 delta_x = ( longitude_dd + 360.0 - west_lon ) / scale_factor;
01073 else
01074 delta_x = ( longitude_dd - west_lon ) / scale_factor;
01075 delta_y = ( latitude_dd - south_lat ) / scale_factor;
01076
01077 _1_minus_delta_x = 1 - delta_x;
01078 _1_minus_delta_y = 1 - delta_y;
01079
01080 w_sw = _1_minus_delta_x * _1_minus_delta_y;
01081 w_se = delta_x * _1_minus_delta_y;
01082 w_ne = delta_x * delta_y;
01083 w_nw = _1_minus_delta_x * delta_y;
01084
01085 *delta_height =
01086 height_sw * w_sw + height_se * w_se + height_ne * w_ne + height_nw * w_nw;
01087 }
01088
01089
01090 void GeoidLibrary::bilinearInterpolate(
01091 double longitude,
01092 double latitude,
01093 double scale_factor,
01094 int num_cols,
01095 int num_rows,
01096 float *height_buffer,
01097 double *delta_height )
01098 {
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115 int index;
01116 int post_x, post_y;
01117 double offset_x, offset_y;
01118 double delta_x, delta_y;
01119 double _1_minus_delta_x, _1_minus_delta_y;
01120 double latitude_dd, longitude_dd;
01121 double height_se, height_sw, height_ne, height_nw;
01122 double w_sw, w_se, w_ne, w_nw;
01123 double south_lat, west_lon;
01124 int end_index = 0;
01125 int max_index = num_rows * num_cols - 1;
01126 char errorStatus[50] = "";
01127
01128 if( ( latitude < -PI_OVER_2 ) || ( latitude > PI_OVER_2 ) )
01129 {
01130 strcat( errorStatus, ErrorMessages::latitude );
01131 }
01132 if( ( longitude < -PI ) || ( longitude > TWO_PI ) )
01133 {
01134 strcat( errorStatus, ErrorMessages::longitude );
01135 }
01136
01137 if( strlen( errorStatus ) > 0)
01138 throw CoordinateConversionException( errorStatus );
01139
01140 latitude_dd = latitude * _180_OVER_PI;
01141 longitude_dd = longitude * _180_OVER_PI;
01142
01143
01144
01145 if( longitude_dd < 0.0 )
01146 {
01147 offset_x = ( longitude_dd + 360.0 ) / scale_factor;
01148 }
01149 else
01150 {
01151 offset_x = longitude_dd / scale_factor;
01152 }
01153 offset_y = ( 90 - latitude_dd ) / scale_factor;
01154
01155
01156
01157
01158 post_x = ( int )( offset_x );
01159 if( ( post_x + 1 ) == num_cols )
01160 post_x--;
01161 post_y = ( int )( offset_y + 1.0e-11 );
01162 if( ( post_y + 1 ) == num_rows )
01163 post_y--;
01164
01165
01166 index = post_y * num_cols + post_x;
01167 if( index < 0 )
01168 height_nw = height_buffer[ 0 ];
01169 else if( index > max_index )
01170 height_nw = height_buffer[ max_index ];
01171 else
01172 height_nw = height_buffer[ index ];
01173
01174 end_index = index + 1;
01175 if( end_index > max_index )
01176 height_ne = height_buffer[ max_index ];
01177 else
01178 height_ne = height_buffer[ end_index ];
01179
01180
01181 index = ( post_y + 1 ) * num_cols + post_x;
01182 if( index < 0 )
01183 height_sw = height_buffer[ 0 ];
01184 else if( index > max_index )
01185 height_sw = height_buffer[ max_index ];
01186 else
01187 height_sw = height_buffer[ index ];
01188
01189
01190 end_index = index + 1;
01191 if( end_index > max_index )
01192 height_se = height_buffer[ max_index ];
01193 else
01194 height_se = height_buffer[ end_index ];
01195
01196 west_lon = post_x * scale_factor;
01197
01198
01199 south_lat = ( 90 - ( post_y * scale_factor ) ) - scale_factor;
01200
01201
01202
01203 if( longitude_dd < 0.0 )
01204 delta_x = ( longitude_dd + 360.0 - west_lon ) / scale_factor;
01205 else
01206 delta_x = ( longitude_dd - west_lon ) / scale_factor;
01207 delta_y = ( latitude_dd - south_lat ) / scale_factor;
01208
01209 _1_minus_delta_x = 1 - delta_x;
01210 _1_minus_delta_y = 1 - delta_y;
01211
01212 w_sw = _1_minus_delta_x * _1_minus_delta_y;
01213 w_se = delta_x * _1_minus_delta_y;
01214 w_ne = delta_x * delta_y;
01215 w_nw = _1_minus_delta_x * delta_y;
01216
01217 *delta_height =
01218 height_sw * w_sw + height_se * w_se + height_ne * w_ne + height_nw * w_nw;
01219 }
01220
01221
01222 void GeoidLibrary::naturalSplineInterpolate(
01223 double longitude,
01224 double latitude,
01225 double scale_factor,
01226 int num_cols,
01227 int num_rows,
01228 int max_index,
01229 float *height_buffer,
01230 double *delta_height )
01231 {
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248 int index;
01249 int post_x, post_y;
01250 int temp_offset_x, temp_offset_y;
01251 double offset_x, offset_y;
01252 double delta_x, delta_y;
01253 double delta_x2, delta_y2;
01254 double _1_minus_delta_x, _1_minus_delta_y;
01255 double _1_minus_delta_x2, _1_minus_delta_y2;
01256 double _3_minus_2_times_1_minus_delta_x, _3_minus_2_times_1_minus_delta_y;
01257 double _3_minus_2_times_delta_x, _3_minus_2_times_delta_y;
01258 double latitude_dd, longitude_dd;
01259 double height_se, height_sw, height_ne, height_nw;
01260 double w_sw, w_se, w_ne, w_nw;
01261 double south_lat, west_lon;
01262 int end_index = 0;
01263 double skip_factor = 1.0;
01264 char errorStatus[50] = "";
01265
01266 if( ( latitude < -PI_OVER_2 ) || ( latitude > PI_OVER_2 ) )
01267 {
01268 strcat( errorStatus, ErrorMessages::latitude );
01269 }
01270 if( ( longitude < -PI ) || ( longitude > TWO_PI ) )
01271 {
01272 strcat( errorStatus, ErrorMessages::longitude );
01273 }
01274
01275 if( strlen( errorStatus ) > 0)
01276 throw CoordinateConversionException( errorStatus );
01277
01278 latitude_dd = latitude * _180_OVER_PI;
01279 longitude_dd = longitude * _180_OVER_PI;
01280
01281
01282
01283 if( longitude_dd < 0.0 )
01284 {
01285 offset_x = ( longitude_dd + 360.0 ) / scale_factor;
01286 }
01287 else
01288 {
01289 offset_x = longitude_dd / scale_factor;
01290 }
01291 offset_y = ( 90.0 - latitude_dd ) / scale_factor;
01292
01293
01294
01295
01296 post_x = ( int ) offset_x;
01297 if ( ( post_x + 1 ) == num_cols)
01298 post_x--;
01299 post_y = ( int )( offset_y + 1.0e-11 );
01300 if ( ( post_y + 1 ) == num_rows)
01301 post_y--;
01302
01303 if( scale_factor == SCALE_FACTOR_30_MINUTES )
01304 {
01305 skip_factor = 2.0;
01306 num_rows = EGM96_ROWS;
01307 num_cols = EGM96_COLS;
01308 }
01309 else if( scale_factor == SCALE_FACTOR_1_DEGREE )
01310 {
01311 skip_factor = 4.0;
01312 num_rows = EGM96_ROWS;
01313 num_cols = EGM96_COLS;
01314 }
01315 else if( scale_factor == SCALE_FACTOR_2_DEGREES )
01316 {
01317 skip_factor = 8.0;
01318 num_rows = EGM96_ROWS;
01319 num_cols = EGM96_COLS;
01320 }
01321
01322 temp_offset_x = ( int )( post_x * skip_factor );
01323 temp_offset_y = ( int )( post_y * skip_factor + 1.0e-11 );
01324 if ( ( temp_offset_x + 1 ) == num_cols )
01325 temp_offset_x--;
01326 if ( ( temp_offset_y + 1 ) == num_rows )
01327 temp_offset_y--;
01328
01329
01330 index = ( int )( temp_offset_y * num_cols + temp_offset_x );
01331 if( index < 0 )
01332 height_nw = height_buffer[ 0 ];
01333 else if( index > max_index )
01334 height_nw = height_buffer[ max_index ];
01335 else
01336 height_nw = height_buffer[ index ];
01337
01338 end_index = index + (int)skip_factor;
01339 if( end_index < 0 )
01340 height_ne = height_buffer[ 0 ];
01341 else if( end_index > max_index )
01342 height_ne = height_buffer[ max_index ];
01343 else
01344 height_ne = height_buffer[ end_index ];
01345
01346
01347 index = ( int )( ( temp_offset_y + skip_factor ) * num_cols + temp_offset_x );
01348 if( index < 0 )
01349 height_sw = height_buffer[ 0 ];
01350 else if( index > max_index )
01351 height_sw = height_buffer[ max_index ];
01352 else
01353 height_sw = height_buffer[ index ];
01354
01355 end_index = index + (int)skip_factor;
01356 if( end_index < 0 )
01357 height_se = height_buffer[ 0 ];
01358 else if( end_index > max_index )
01359 height_se = height_buffer[ max_index ];
01360 else
01361 height_se = height_buffer[ end_index ];
01362
01363 west_lon = post_x * scale_factor;
01364
01365
01366 south_lat = ( 90 - ( post_y * scale_factor ) ) - scale_factor;
01367
01368
01369
01370 if( longitude_dd < 0.0 )
01371 delta_x = ( longitude_dd + 360.0 - west_lon ) / scale_factor;
01372 else
01373 delta_x = ( longitude_dd - west_lon ) / scale_factor;
01374 delta_y = ( latitude_dd - south_lat ) / scale_factor;
01375
01376 delta_x2 = delta_x * delta_x;
01377 delta_y2 = delta_y * delta_y;
01378
01379 _1_minus_delta_x = 1 - delta_x;
01380 _1_minus_delta_y = 1 - delta_y;
01381
01382 _1_minus_delta_x2 = _1_minus_delta_x * _1_minus_delta_x;
01383 _1_minus_delta_y2 = _1_minus_delta_y * _1_minus_delta_y;
01384
01385 _3_minus_2_times_1_minus_delta_x = 3 - 2 * _1_minus_delta_x;
01386 _3_minus_2_times_1_minus_delta_y = 3 - 2 * _1_minus_delta_y;
01387
01388 _3_minus_2_times_delta_x = 3 - 2 * delta_x;
01389 _3_minus_2_times_delta_y = 3 - 2 * delta_y;
01390
01391 w_sw = _1_minus_delta_x2 * _1_minus_delta_y2 *
01392 ( _3_minus_2_times_1_minus_delta_x * _3_minus_2_times_1_minus_delta_y );
01393
01394 w_se = delta_x2 * _1_minus_delta_y2 *
01395 ( _3_minus_2_times_delta_x * _3_minus_2_times_1_minus_delta_y );
01396
01397 w_ne = delta_x2 * delta_y2 *
01398 ( _3_minus_2_times_delta_x * _3_minus_2_times_delta_y );
01399
01400 w_nw = _1_minus_delta_x2 * delta_y2 *
01401 ( _3_minus_2_times_1_minus_delta_x * _3_minus_2_times_delta_y );
01402
01403 *delta_height =
01404 height_sw * w_sw + height_se * w_se + height_ne * w_ne + height_nw * w_nw;
01405 }
01406
01407