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