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 <ctype.h>
00081 #include <math.h>
00082 #include <string.h>
00083 #include "UPS.h"
00084 #include "UTM.h"
00085 #include "USNG.h"
00086 #include "EllipsoidParameters.h"
00087 #include "MGRSorUSNGCoordinates.h"
00088 #include "GeodeticCoordinates.h"
00089 #include "UPSCoordinates.h"
00090 #include "UTMCoordinates.h"
00091 #include "CoordinateConversionException.h"
00092 #include "ErrorMessages.h"
00093 #include "WarningMessages.h"
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112 using namespace MSP::CCS;
00113
00114
00115
00116
00117
00118
00119 #define EPSILON 1.75e-7
00120
00121 const int LETTER_A = 0;
00122 const int LETTER_B = 1;
00123 const int LETTER_C = 2;
00124 const int LETTER_D = 3;
00125 const int LETTER_E = 4;
00126 const int LETTER_F = 5;
00127 const int LETTER_G = 6;
00128 const int LETTER_H = 7;
00129 const int LETTER_I = 8;
00130 const int LETTER_J = 9;
00131 const int LETTER_K = 10;
00132 const int LETTER_L = 11;
00133 const int LETTER_M = 12;
00134 const int LETTER_N = 13;
00135 const int LETTER_O = 14;
00136 const int LETTER_P = 15;
00137 const int LETTER_Q = 16;
00138 const int LETTER_R = 17;
00139 const int LETTER_S = 18;
00140 const int LETTER_T = 19;
00141 const int LETTER_U = 20;
00142 const int LETTER_V = 21;
00143 const int LETTER_W = 22;
00144 const int LETTER_X = 23;
00145 const int LETTER_Y = 24;
00146 const int LETTER_Z = 25;
00147 const double ONEHT = 100000.e0;
00148 const double TWOMIL = 2000000.e0;
00149 const double PI = 3.14159265358979323e0;
00150 const double PI_OVER_2 = (PI / 2.0e0);
00151 const double PI_OVER_180 = (PI / 180.0e0);
00152
00153 const double MIN_EASTING = 100000.0;
00154 const double MAX_EASTING = 900000.0;
00155 const double MIN_NORTHING = 0.0;
00156 const double MAX_NORTHING = 10000000.0;
00157 const int MAX_PRECISION = 5;
00158 const double MIN_USNG_NON_POLAR_LAT = -80.0 * ( PI / 180.0 );
00159 const double MAX_USNG_NON_POLAR_LAT = 84.0 * ( PI / 180.0 );
00160
00161 const double MIN_EAST_NORTH = 0.0;
00162 const double MAX_EAST_NORTH = 3999999.0;
00163
00164 const double _6 = (6.0 * (PI / 180.0));
00165 const double _8 = (8.0 * (PI / 180.0));
00166 const double _72 = (72.0 * (PI / 180.0));
00167 const double _80 = (80.0 * (PI / 180.0));
00168 const double _80_5 = (80.5 * (PI / 180.0));
00169 const double _84_5 = (84.5 * (PI / 180.0));
00170
00171 #define _500000 500000.0
00172
00173 struct Latitude_Band
00174 {
00175 long letter;
00176 double min_northing;
00177 double north;
00178 double south;
00179 double northing_offset;
00180 };
00181
00182 const Latitude_Band Latitude_Band_Table[20] =
00183 {{LETTER_C, 1100000.0, -72.0, -80.5, 0.0},
00184 {LETTER_D, 2000000.0, -64.0, -72.0, 2000000.0},
00185 {LETTER_E, 2800000.0, -56.0, -64.0, 2000000.0},
00186 {LETTER_F, 3700000.0, -48.0, -56.0, 2000000.0},
00187 {LETTER_G, 4600000.0, -40.0, -48.0, 4000000.0},
00188 {LETTER_H, 5500000.0, -32.0, -40.0, 4000000.0},
00189 {LETTER_J, 6400000.0, -24.0, -32.0, 6000000.0},
00190 {LETTER_K, 7300000.0, -16.0, -24.0, 6000000.0},
00191 {LETTER_L, 8200000.0, -8.0, -16.0, 8000000.0},
00192 {LETTER_M, 9100000.0, 0.0, -8.0, 8000000.0},
00193 {LETTER_N, 0.0, 8.0, 0.0, 0.0},
00194 {LETTER_P, 800000.0, 16.0, 8.0, 0.0},
00195 {LETTER_Q, 1700000.0, 24.0, 16.0, 0.0},
00196 {LETTER_R, 2600000.0, 32.0, 24.0, 2000000.0},
00197 {LETTER_S, 3500000.0, 40.0, 32.0, 2000000.0},
00198 {LETTER_T, 4400000.0, 48.0, 40.0, 4000000.0},
00199 {LETTER_U, 5300000.0, 56.0, 48.0, 4000000.0},
00200 {LETTER_V, 6200000.0, 64.0, 56.0, 6000000.0},
00201 {LETTER_W, 7000000.0, 72.0, 64.0, 6000000.0},
00202 {LETTER_X, 7900000.0, 84.5, 72.0, 6000000.0}};
00203
00204
00205 struct UPS_Constant
00206 {
00207 long letter;
00208 long ltr2_low_value;
00209 long ltr2_high_value;
00210 long ltr3_high_value;
00211 double false_easting;
00212 double false_northing;
00213 };
00214
00215 const UPS_Constant UPS_Constant_Table[4] =
00216 {{LETTER_A, LETTER_J, LETTER_Z, LETTER_Z, 800000.0, 800000.0},
00217 {LETTER_B, LETTER_A, LETTER_R, LETTER_Z, 2000000.0, 800000.0},
00218 {LETTER_Y, LETTER_J, LETTER_Z, LETTER_P, 800000.0, 1300000.0},
00219 {LETTER_Z, LETTER_A, LETTER_J, LETTER_P, 2000000.0, 1300000.0}};
00220
00221
00222
00223
00224
00225
00226
00227 void makeUSNGString(
00228 char* USNGString,
00229 long zone,
00230 int letters[USNG_LETTERS],
00231 double easting,
00232 double northing,
00233 long precision )
00234 {
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247 long i;
00248 long j;
00249 double divisor;
00250 long east;
00251 long north;
00252 char alphabet[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
00253
00254 i = 0;
00255 if (zone)
00256 i = sprintf (USNGString+i,"%2.2ld",zone);
00257 else
00258 strncpy(USNGString, " ", 2);
00259
00260 for (j=0;j<3;j++)
00261 USNGString[i++] = alphabet[letters[j]];
00262 divisor = pow (10.0, (5.0 - precision));
00263 easting = fmod (easting, 100000.0);
00264 if (easting >= 99999.5)
00265 easting = 99999.0;
00266 east = (long)(easting/divisor);
00267 i += sprintf (USNGString+i, "%*.*ld", precision, precision, east);
00268 northing = fmod (northing, 100000.0);
00269 if (northing >= 99999.5)
00270 northing = 99999.0;
00271 north = (long)(northing/divisor);
00272 i += sprintf (USNGString+i, "%*.*ld", precision, precision, north);
00273 }
00274
00275
00276 void breakUSNGString(
00277 char* USNGString,
00278 long* zone,
00279 long letters[USNG_LETTERS],
00280 double* easting,
00281 double* northing,
00282 long* precision )
00283 {
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296 long num_digits;
00297 long num_letters;
00298 long i = 0;
00299 long j = 0;
00300
00301 while (USNGString[i] == ' ')
00302 i++;
00303 j = i;
00304 while (isdigit(USNGString[i]))
00305 i++;
00306 num_digits = i - j;
00307 if (num_digits <= 2)
00308 if (num_digits > 0)
00309 {
00310 char zone_string[3];
00311
00312 strncpy (zone_string, USNGString+j, 2);
00313 zone_string[2] = 0;
00314 sscanf (zone_string, "%ld", zone);
00315 if ((*zone < 1) || (*zone > 60))
00316 throw CoordinateConversionException( ErrorMessages::usngString );
00317 }
00318 else
00319 *zone = 0;
00320 else
00321 throw CoordinateConversionException( ErrorMessages::usngString );
00322 j = i;
00323
00324 while (isalpha(USNGString[i]))
00325 i++;
00326 num_letters = i - j;
00327 if (num_letters == 3)
00328 {
00329
00330 letters[0] = (toupper(USNGString[j]) - (long)'A');
00331 if ((letters[0] == LETTER_I) || (letters[0] == LETTER_O))
00332 throw CoordinateConversionException( ErrorMessages::usngString );
00333 letters[1] = (toupper(USNGString[j+1]) - (long)'A');
00334 if ((letters[1] == LETTER_I) || (letters[1] == LETTER_O))
00335 throw CoordinateConversionException( ErrorMessages::usngString );
00336 letters[2] = (toupper(USNGString[j+2]) - (long)'A');
00337 if ((letters[2] == LETTER_I) || (letters[2] == LETTER_O))
00338 throw CoordinateConversionException( ErrorMessages::usngString );
00339 }
00340 else
00341 throw CoordinateConversionException( ErrorMessages::usngString );
00342 j = i;
00343 while (isdigit(USNGString[i]))
00344 i++;
00345 num_digits = i - j;
00346 if ((num_digits <= 10) && (num_digits%2 == 0))
00347 {
00348 long n;
00349 char east_string[6];
00350 char north_string[6];
00351 long east;
00352 long north;
00353 double multiplier;
00354
00355 n = num_digits/2;
00356 *precision = n;
00357 if (n > 0)
00358 {
00359 strncpy (east_string, USNGString+j, n);
00360 east_string[n] = 0;
00361 sscanf (east_string, "%ld", &east);
00362 strncpy (north_string, USNGString+j+n, n);
00363 north_string[n] = 0;
00364 sscanf (north_string, "%ld", &north);
00365 multiplier = pow (10.0, 5.0 - n);
00366 *easting = east * multiplier;
00367 *northing = north * multiplier;
00368 }
00369 else
00370 {
00371 *easting = 0.0;
00372 *northing = 0.0;
00373 }
00374 }
00375 else
00376 throw CoordinateConversionException( ErrorMessages::usngString );
00377 }
00378
00379
00380
00381
00382
00383
00384
00385 USNG::USNG( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, char* ellipsoidCode ) :
00386 CoordinateSystem(),
00387 ups( 0 ),
00388 utm( 0 )
00389 {
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400 double inv_f = 1 / ellipsoidFlattening;
00401 char errorStatus[500] = "";
00402
00403 if (ellipsoidSemiMajorAxis <= 0.0)
00404 {
00405 strcat( errorStatus, ErrorMessages::semiMajorAxis );
00406 }
00407 if ((inv_f < 250) || (inv_f > 350))
00408 {
00409 strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
00410 }
00411
00412 if( strlen( errorStatus ) > 0)
00413 throw CoordinateConversionException( errorStatus );
00414
00415 semiMajorAxis = ellipsoidSemiMajorAxis;
00416 flattening = ellipsoidFlattening;
00417
00418 strncpy (USNGEllipsoidCode, ellipsoidCode, 2);
00419 USNGEllipsoidCode[2] = '\0';
00420
00421 ups = new UPS( semiMajorAxis, flattening );
00422
00423 utm = new UTM( semiMajorAxis, flattening, 0 );
00424 }
00425
00426
00427 USNG::USNG( const USNG &u )
00428 {
00429 ups = new UPS( *( u.ups ) );
00430 utm = new UTM( *( u.utm ) );
00431
00432 semiMajorAxis = u.semiMajorAxis;
00433 flattening = u.flattening;
00434 strcpy( USNGEllipsoidCode, u.USNGEllipsoidCode );
00435 }
00436
00437
00438 USNG::~USNG()
00439 {
00440 delete ups;
00441 ups = 0;
00442
00443 delete utm;
00444 utm = 0;
00445 }
00446
00447
00448 USNG& USNG::operator=( const USNG &u )
00449 {
00450 if( this != &u )
00451 {
00452 ups->operator=( *u.ups );
00453 utm->operator=( *u.utm );
00454
00455 semiMajorAxis = u.semiMajorAxis;
00456 flattening = u.flattening;
00457 strcpy( USNGEllipsoidCode, u.USNGEllipsoidCode );
00458 }
00459
00460 return *this;
00461 }
00462
00463
00464 EllipsoidParameters* USNG::getParameters() const
00465 {
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475 return new EllipsoidParameters(
00476 semiMajorAxis, flattening, (char*)USNGEllipsoidCode );
00477 }
00478
00479
00480 MSP::CCS::MGRSorUSNGCoordinates* USNG::convertFromGeodetic(
00481 MSP::CCS::GeodeticCoordinates* geodeticCoordinates,
00482 long precision )
00483 {
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497 MGRSorUSNGCoordinates* mgrsorUSNGCoordinates = 0;
00498 UTMCoordinates* utmCoordinates = 0;
00499 UPSCoordinates* upsCoordinates = 0;
00500
00501 char errorStatus[50] = "";
00502
00503 double latitude = geodeticCoordinates->latitude();
00504 double longitude = geodeticCoordinates->longitude();
00505
00506 if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
00507 {
00508 strcat( errorStatus, ErrorMessages::latitude );
00509 }
00510 if ((longitude < -PI) || (longitude > (2*PI)))
00511 {
00512 strcat( errorStatus, ErrorMessages::longitude );
00513 }
00514 if ((precision < 0) || (precision > MAX_PRECISION))
00515 strcat( errorStatus, ErrorMessages::precision );
00516
00517 if( strlen( errorStatus ) > 0)
00518 throw CoordinateConversionException( errorStatus );
00519
00520
00521
00522
00523 try
00524 {
00525 if((latitude >= MIN_USNG_NON_POLAR_LAT) &&
00526 (latitude < MAX_USNG_NON_POLAR_LAT))
00527 {
00528 utmCoordinates = utm->convertFromGeodetic( geodeticCoordinates );
00529 mgrsorUSNGCoordinates = fromUTM(
00530 utmCoordinates, longitude, latitude, precision );
00531 delete utmCoordinates;
00532 }
00533 else
00534 {
00535 upsCoordinates = ups->convertFromGeodetic( geodeticCoordinates );
00536 mgrsorUSNGCoordinates = fromUPS( upsCoordinates, precision );
00537 delete upsCoordinates;
00538 }
00539 }
00540 catch ( CoordinateConversionException e ) {
00541 delete utmCoordinates;
00542 delete upsCoordinates;
00543 throw e;
00544 }
00545
00546 return mgrsorUSNGCoordinates;
00547 }
00548
00549
00550 MSP::CCS::GeodeticCoordinates* USNG::convertToGeodetic(
00551 MSP::CCS::MGRSorUSNGCoordinates* usngCoordinates )
00552 {
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565 long zone;
00566 long letters[USNG_LETTERS];
00567 double usng_easting;
00568 double usng_northing;
00569 long precision;
00570 GeodeticCoordinates* geodeticCoordinates = 0;
00571 UTMCoordinates* utmCoordinates = 0;
00572 UPSCoordinates* upsCoordinates = 0;
00573
00574 breakUSNGString(
00575 usngCoordinates->MGRSString(),
00576 &zone, letters, &usng_easting, &usng_northing, &precision );
00577
00578 try
00579 {
00580 if( zone )
00581 {
00582 utmCoordinates = toUTM(
00583 zone, letters, usng_easting, usng_northing, precision );
00584 geodeticCoordinates = utm->convertToGeodetic( utmCoordinates );
00585 delete utmCoordinates;
00586 }
00587 else
00588 {
00589 upsCoordinates = toUPS( letters, usng_easting, usng_northing );
00590 geodeticCoordinates = ups->convertToGeodetic( upsCoordinates );
00591 delete upsCoordinates;
00592 }
00593 }
00594 catch ( CoordinateConversionException e )
00595 {
00596 delete utmCoordinates;
00597 delete upsCoordinates;
00598 throw e;
00599 }
00600
00601 return geodeticCoordinates;
00602 }
00603
00604
00605 MSP::CCS::MGRSorUSNGCoordinates* USNG::convertFromUTM (
00606 UTMCoordinates* utmCoordinates, long precision )
00607 {
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622 char errorStatus[50] = "";
00623
00624 long zone = utmCoordinates->zone();
00625 char hemisphere = utmCoordinates->hemisphere();
00626 double easting = utmCoordinates->easting();
00627 double northing = utmCoordinates->northing();
00628
00629 if ((zone < 1) || (zone > 60))
00630 strcat( errorStatus, ErrorMessages::zone );
00631 if ((hemisphere != 'S') && (hemisphere != 'N'))
00632 strcat( errorStatus, ErrorMessages::hemisphere );
00633 if ((easting < MIN_EASTING) || (easting > MAX_EASTING))
00634 strcat( errorStatus, ErrorMessages::easting );
00635 if ((northing < MIN_NORTHING) || (northing > MAX_NORTHING))
00636 strcat( errorStatus, ErrorMessages::northing );
00637 if ((precision < 0) || (precision > MAX_PRECISION))
00638 strcat( errorStatus, ErrorMessages::precision );
00639
00640 if( strlen( errorStatus ) > 0)
00641 throw CoordinateConversionException( errorStatus );
00642
00643 GeodeticCoordinates* geodeticCoordinates = utm->convertToGeodetic(
00644 utmCoordinates );
00645
00646
00647
00648
00649 MGRSorUSNGCoordinates* mgrsorUSNGCoordinates = 0;
00650 UPSCoordinates* upsCoordinates = 0;
00651 double latitude = geodeticCoordinates->latitude();
00652
00653 try
00654 {
00655 if((latitude >= (MIN_USNG_NON_POLAR_LAT - EPSILON)) &&
00656 (latitude < (MAX_USNG_NON_POLAR_LAT + EPSILON)))
00657 mgrsorUSNGCoordinates = fromUTM(
00658 utmCoordinates, geodeticCoordinates->longitude(),
00659 latitude, precision );
00660 else
00661 {
00662 upsCoordinates = ups->convertFromGeodetic( geodeticCoordinates );
00663 mgrsorUSNGCoordinates = fromUPS( upsCoordinates, precision );
00664 }
00665 }
00666 catch ( CoordinateConversionException e )
00667 {
00668 delete geodeticCoordinates;
00669 delete upsCoordinates;
00670 throw e;
00671 }
00672
00673 delete geodeticCoordinates;
00674 delete upsCoordinates;
00675
00676 return mgrsorUSNGCoordinates;
00677 }
00678
00679
00680 MSP::CCS::UTMCoordinates* USNG::convertToUTM(
00681 MSP::CCS::MGRSorUSNGCoordinates* mgrsorUSNGCoordinates )
00682 {
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696 long zone;
00697 long letters[USNG_LETTERS];
00698 double usng_easting, usng_northing;
00699 long precision;
00700 UTMCoordinates* utmCoordinates = 0;
00701 GeodeticCoordinates* geodeticCoordinates = 0;
00702 UPSCoordinates* upsCoordinates = 0;
00703 char errorStatus[50] = "";
00704
00705 breakUSNGString( mgrsorUSNGCoordinates->MGRSString(),
00706 &zone, letters, &usng_easting, &usng_northing, &precision );
00707 try
00708 {
00709 if (zone)
00710 {
00711 utmCoordinates = toUTM(
00712 zone, letters, usng_easting, usng_northing, precision );
00713
00714 geodeticCoordinates = utm->convertToGeodetic( utmCoordinates );
00715 }
00716 else
00717 {
00718 upsCoordinates = toUPS( letters, usng_easting, usng_northing );
00719 geodeticCoordinates = ups->convertToGeodetic( upsCoordinates );
00720 utmCoordinates = utm->convertFromGeodetic( geodeticCoordinates );
00721 }
00722 }
00723 catch ( CoordinateConversionException e )
00724 {
00725 delete utmCoordinates;
00726 delete geodeticCoordinates;
00727 delete upsCoordinates;
00728 throw e;
00729 }
00730
00731 delete geodeticCoordinates;
00732 delete upsCoordinates;
00733
00734 return utmCoordinates;
00735 }
00736
00737
00738 MSP::CCS::MGRSorUSNGCoordinates* USNG::convertFromUPS(
00739 MSP::CCS::UPSCoordinates* upsCoordinates,
00740 long precision )
00741 {
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755 int index = 0;
00756 char errorStatus[50] = "";
00757
00758 char hemisphere = upsCoordinates->hemisphere();
00759 double easting = upsCoordinates->easting();
00760 double northing = upsCoordinates->northing();
00761
00762 if ((hemisphere != 'N') && (hemisphere != 'S'))
00763 strcat( errorStatus, ErrorMessages::hemisphere );
00764 if ((easting < MIN_EAST_NORTH) || (easting > MAX_EAST_NORTH))
00765 strcat( errorStatus, ErrorMessages::easting );
00766 if ((northing < MIN_EAST_NORTH) || (northing > MAX_EAST_NORTH))
00767 strcat( errorStatus, ErrorMessages::northing );
00768 if ((precision < 0) || (precision > MAX_PRECISION))
00769 strcat( errorStatus, ErrorMessages::precision );
00770
00771 if( strlen( errorStatus ) > 0)
00772 throw CoordinateConversionException( errorStatus );
00773
00774 GeodeticCoordinates* geodeticCoordinates = ups->convertToGeodetic(
00775 upsCoordinates );
00776
00777
00778
00779
00780 MGRSorUSNGCoordinates* mgrsorUSNGCoordinates = 0;
00781 UTMCoordinates* utmCoordinates = 0;
00782 double latitude = geodeticCoordinates->latitude();
00783
00784 try
00785 {
00786 if((latitude < (MIN_USNG_NON_POLAR_LAT + EPSILON)) ||
00787 (latitude >= (MAX_USNG_NON_POLAR_LAT - EPSILON)))
00788 mgrsorUSNGCoordinates = fromUPS( upsCoordinates, precision );
00789 else
00790 {
00791 utmCoordinates = utm->convertFromGeodetic( geodeticCoordinates );
00792 double longitude = geodeticCoordinates->longitude();
00793 mgrsorUSNGCoordinates = fromUTM(
00794 utmCoordinates, longitude, latitude, precision );
00795 }
00796 }
00797 catch ( CoordinateConversionException e )
00798 {
00799 delete geodeticCoordinates;
00800 delete utmCoordinates;
00801 throw e;
00802 }
00803
00804 delete geodeticCoordinates;
00805 delete utmCoordinates;
00806
00807 return mgrsorUSNGCoordinates;
00808 }
00809
00810
00811 MSP::CCS::UPSCoordinates* USNG::convertToUPS(
00812 MSP::CCS::MGRSorUSNGCoordinates* mgrsorUSNGCoordinates )
00813 {
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826 long zone;
00827 long letters[USNG_LETTERS];
00828 long precision;
00829 double usng_easting;
00830 double usng_northing;
00831 int index = 0;
00832 UPSCoordinates* upsCoordinates = 0;
00833 GeodeticCoordinates* geodeticCoordinates = 0;
00834 UTMCoordinates* utmCoordinates = 0;
00835
00836 breakUSNGString( mgrsorUSNGCoordinates->MGRSString(),
00837 &zone, letters, &usng_easting, &usng_northing, &precision );
00838 try
00839 {
00840 if( !zone )
00841 {
00842 upsCoordinates = toUPS( letters, usng_easting, usng_northing );
00843
00844 geodeticCoordinates = ups->convertToGeodetic( upsCoordinates );
00845 }
00846 else
00847 {
00848 utmCoordinates = toUTM(
00849 zone, letters, usng_easting, usng_northing, precision );
00850 geodeticCoordinates = utm->convertToGeodetic( utmCoordinates );
00851 upsCoordinates = ups->convertFromGeodetic( geodeticCoordinates );
00852 }
00853 }
00854 catch ( CoordinateConversionException e )
00855 {
00856 delete geodeticCoordinates;
00857 delete utmCoordinates;
00858 delete upsCoordinates;
00859 throw e;
00860 }
00861
00862 delete geodeticCoordinates;
00863 delete utmCoordinates;
00864
00865 return upsCoordinates;
00866 }
00867
00868
00869 MSP::CCS::MGRSorUSNGCoordinates* USNG::fromUTM(
00870 MSP::CCS::UTMCoordinates* utmCoordinates,
00871 double longitude,
00872 double latitude,
00873 long precision )
00874 {
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887 double pattern_offset;
00888 double grid_northing;
00889 long ltr2_low_value;
00890 long ltr2_high_value;
00891 int letters[USNG_LETTERS];
00892 char USNGString[21];
00893 long override = 0;
00894 long natural_zone;
00895
00896 long zone = utmCoordinates->zone();
00897 char hemisphere = utmCoordinates->hemisphere();
00898 double easting = utmCoordinates->easting();
00899 double northing = utmCoordinates->northing();
00900
00901 getLatitudeLetter( latitude, &letters[0] );
00902
00903
00904
00905 if (longitude < PI)
00906 natural_zone = (long)(31 + ((longitude) / _6));
00907 else
00908 natural_zone = (long)(((longitude) / _6) - 29);
00909
00910 if (natural_zone > 60)
00911 natural_zone = 1;
00912
00913 if (zone != natural_zone)
00914 {
00915 UTM utmOverride( semiMajorAxis, flattening, natural_zone );
00916 GeodeticCoordinates geodeticCoordinates(
00917 CoordinateType::geodetic, longitude, latitude );
00918 UTMCoordinates* utmCoordinatesOverride = utmOverride.convertFromGeodetic(
00919 &geodeticCoordinates );
00920
00921 zone = utmCoordinatesOverride->zone();
00922 hemisphere = utmCoordinatesOverride->hemisphere();
00923 easting = utmCoordinatesOverride->easting();
00924 northing = utmCoordinatesOverride->northing();
00925
00926 delete utmCoordinatesOverride;
00927 utmCoordinatesOverride = 0;
00928 }
00929
00930
00931 if (letters[0] == LETTER_V)
00932 {
00933 if ((zone == 31) && (easting >= _500000))
00934 override = 32;
00935 }
00936 else if (letters[0] == LETTER_X)
00937 {
00938 if ((zone == 32) && (easting < _500000))
00939 override = 31;
00940 else if (((zone == 32) && (easting >= _500000)) ||
00941 ((zone == 34) && (easting < _500000)))
00942 override = 33;
00943 else if (((zone == 34) && (easting >= _500000)) ||
00944 ((zone == 36) && (easting < _500000)))
00945 override = 35;
00946 else if ((zone == 36) && (easting >= _500000))
00947 override = 37;
00948 }
00949
00950 if (override)
00951 {
00952 UTM utmOverride( semiMajorAxis, flattening, override );
00953 GeodeticCoordinates geodeticCoordinates(
00954 CoordinateType::geodetic, longitude, latitude );
00955 UTMCoordinates* utmCoordinatesOverride =
00956 utmOverride.convertFromGeodetic( &geodeticCoordinates );
00957
00958 zone = utmCoordinatesOverride->zone();
00959 hemisphere = utmCoordinatesOverride->hemisphere();
00960 easting = utmCoordinatesOverride->easting();
00961 northing = utmCoordinatesOverride->northing();
00962
00963 delete utmCoordinatesOverride;
00964 utmCoordinatesOverride = 0;
00965 }
00966
00967
00968 double divisor = pow (10.0, (5.0 - precision));
00969 easting = (long)(easting/divisor) * divisor;
00970 northing = (long)(northing/divisor) * divisor;
00971
00972 if( latitude <= 0.0 && northing == 1.0e7)
00973 {
00974 latitude = 0.0;
00975 northing = 0.0;
00976 }
00977
00978 getGridValues( zone, <r2_low_value, <r2_high_value, &pattern_offset );
00979
00980 grid_northing = northing;
00981
00982 while (grid_northing >= TWOMIL)
00983 {
00984 grid_northing = grid_northing - TWOMIL;
00985 }
00986 grid_northing = grid_northing + pattern_offset;
00987 if(grid_northing >= TWOMIL)
00988 grid_northing = grid_northing - TWOMIL;
00989
00990 letters[2] = (long)(grid_northing / ONEHT);
00991 if (letters[2] > LETTER_H)
00992 letters[2] = letters[2] + 1;
00993
00994 if (letters[2] > LETTER_N)
00995 letters[2] = letters[2] + 1;
00996
00997 letters[1] = ltr2_low_value + ((long)(easting / ONEHT) -1);
00998 if ((ltr2_low_value == LETTER_J) && (letters[1] > LETTER_N))
00999 letters[1] = letters[1] + 1;
01000
01001 makeUSNGString( USNGString, zone, letters, easting, northing, precision );
01002
01003 return new MGRSorUSNGCoordinates( CoordinateType::usNationalGrid, USNGString);
01004 }
01005
01006
01007 MSP::CCS::UTMCoordinates* USNG::toUTM(
01008 long zone,
01009 long letters[USNG_LETTERS],
01010 double easting,
01011 double northing,
01012 long precision )
01013 {
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027 char hemisphere;
01028 double min_northing;
01029 double northing_offset;
01030 long ltr2_low_value;
01031 long ltr2_high_value;
01032 double pattern_offset;
01033 double upper_lat_limit;
01034 double lower_lat_limit;
01035 double grid_easting;
01036 double grid_northing;
01037 double temp_grid_northing = 0.0;
01038 double fabs_grid_northing = 0.0;
01039 double latitude = 0.0;
01040 double longitude = 0.0;
01041 double divisor = 1.0;
01042 UTMCoordinates* utmCoordinates = 0;
01043 char errorStatus[50] = "";
01044
01045 if((letters[0] == LETTER_X) && ((zone == 32) || (zone == 34) || (zone == 36)))
01046 throw CoordinateConversionException( ErrorMessages::usngString );
01047 else if ((letters[0] == LETTER_V) && (zone == 31) && (letters[1] > LETTER_D))
01048 throw CoordinateConversionException( ErrorMessages::usngString );
01049 else
01050 {
01051 if (letters[0] < LETTER_N)
01052 hemisphere = 'S';
01053 else
01054 hemisphere = 'N';
01055
01056 getGridValues(zone, <r2_low_value, <r2_high_value, &pattern_offset);
01057
01058
01059
01060
01061 if((letters[1] < ltr2_low_value) ||
01062 (letters[1] > ltr2_high_value) ||
01063 (letters[2] > LETTER_V))
01064 throw CoordinateConversionException( ErrorMessages::usngString );
01065
01066 grid_easting = (double)((letters[1]) - ltr2_low_value + 1) * ONEHT;
01067 if ((ltr2_low_value == LETTER_J) && (letters[1] > LETTER_O))
01068 grid_easting = grid_easting - ONEHT;
01069
01070 double row_letter_northing = (double)(letters[2]) * ONEHT;
01071 if (letters[2] > LETTER_O)
01072 row_letter_northing = row_letter_northing - ONEHT;
01073
01074 if (letters[2] > LETTER_I)
01075 row_letter_northing = row_letter_northing - ONEHT;
01076
01077 if (row_letter_northing >= TWOMIL)
01078 row_letter_northing = row_letter_northing - TWOMIL;
01079
01080 getLatitudeBandMinNorthing(letters[0], &min_northing, &northing_offset);
01081
01082 grid_northing = row_letter_northing - pattern_offset;
01083 if(grid_northing < 0)
01084 grid_northing += TWOMIL;
01085
01086 grid_northing += northing_offset;
01087
01088 if(grid_northing < min_northing)
01089 grid_northing += TWOMIL;
01090
01091 easting = grid_easting + easting;
01092 northing = grid_northing + northing;
01093
01094 utmCoordinates = new UTMCoordinates(
01095 CoordinateType::universalTransverseMercator,
01096 zone, hemisphere, easting, northing );
01097
01098
01099 GeodeticCoordinates* geodeticCoordinates = 0;
01100 try
01101 {
01102 geodeticCoordinates = utm->convertToGeodetic( utmCoordinates );
01103 }
01104 catch ( CoordinateConversionException e )
01105 {
01106 delete utmCoordinates;
01107 throw e;
01108 }
01109
01110 divisor = pow (10.0, (double)precision);
01111 getLatitudeRange(letters[0], &upper_lat_limit, &lower_lat_limit);
01112
01113 double latitude = geodeticCoordinates->latitude();
01114
01115 delete geodeticCoordinates;
01116
01117 if(!(((lower_lat_limit - PI_OVER_180/divisor) <= latitude) &&
01118 (latitude <= (upper_lat_limit + PI_OVER_180/divisor))))
01119 utmCoordinates->setWarningMessage( MSP::CCS::WarningMessages::latitude);
01120 }
01121
01122 return utmCoordinates;
01123 }
01124
01125
01126 MSP::CCS::MGRSorUSNGCoordinates* USNG::fromUPS(
01127 MSP::CCS::UPSCoordinates* upsCoordinates,
01128 long precision )
01129 {
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142 double false_easting;
01143 double false_northing;
01144 double grid_easting;
01145 double grid_northing;
01146 long ltr2_low_value;
01147 int letters[USNG_LETTERS];
01148 double divisor;
01149 int index = 0;
01150 char USNGString[21];
01151
01152 char hemisphere = upsCoordinates->hemisphere();
01153 double easting = upsCoordinates->easting();
01154 double northing = upsCoordinates->northing();
01155
01156 divisor = pow (10.0, (5.0 - precision));
01157 easting = (long)(easting/divisor) * divisor;
01158 northing = (long)(northing/divisor) * divisor;
01159
01160 if (hemisphere == 'N')
01161 {
01162 if (easting >= TWOMIL)
01163 letters[0] = LETTER_Z;
01164 else
01165 letters[0] = LETTER_Y;
01166
01167 index = letters[0] - 22;
01168 ltr2_low_value = UPS_Constant_Table[index].ltr2_low_value;
01169 false_easting = UPS_Constant_Table[index].false_easting;
01170 false_northing = UPS_Constant_Table[index].false_northing;
01171 }
01172 else
01173 {
01174 if (easting >= TWOMIL)
01175 letters[0] = LETTER_B;
01176 else
01177 letters[0] = LETTER_A;
01178
01179 ltr2_low_value = UPS_Constant_Table[letters[0]].ltr2_low_value;
01180 false_easting = UPS_Constant_Table[letters[0]].false_easting;
01181 false_northing = UPS_Constant_Table[letters[0]].false_northing;
01182 }
01183
01184 grid_northing = northing;
01185 grid_northing = grid_northing - false_northing;
01186 letters[2] = (long)(grid_northing / ONEHT);
01187
01188 if (letters[2] > LETTER_H)
01189 letters[2] = letters[2] + 1;
01190
01191 if (letters[2] > LETTER_N)
01192 letters[2] = letters[2] + 1;
01193
01194 grid_easting = easting;
01195 grid_easting = grid_easting - false_easting;
01196 letters[1] = ltr2_low_value + ((long)(grid_easting / ONEHT));
01197
01198 if (easting < TWOMIL)
01199 {
01200 if (letters[1] > LETTER_L)
01201 letters[1] = letters[1] + 3;
01202
01203 if (letters[1] > LETTER_U)
01204 letters[1] = letters[1] + 2;
01205 }
01206 else
01207 {
01208 if (letters[1] > LETTER_C)
01209 letters[1] = letters[1] + 2;
01210
01211 if (letters[1] > LETTER_H)
01212 letters[1] = letters[1] + 1;
01213
01214 if (letters[1] > LETTER_L)
01215 letters[1] = letters[1] + 3;
01216 }
01217
01218 makeUSNGString( USNGString, 0, letters, easting, northing, precision );
01219
01220 return new MGRSorUSNGCoordinates( CoordinateType::usNationalGrid, USNGString);
01221 }
01222
01223
01224 MSP::CCS::UPSCoordinates* USNG::toUPS(
01225 long letters[USNG_LETTERS],
01226 double easting,
01227 double northing )
01228 {
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241 long ltr2_high_value;
01242 long ltr3_high_value;
01243 long ltr2_low_value;
01244 double false_easting;
01245 double false_northing;
01246 double grid_easting;
01247 double grid_northing;
01248 char hemisphere;
01249 int index = 0;
01250
01251 if ((letters[0] == LETTER_Y) || (letters[0] == LETTER_Z))
01252 {
01253 hemisphere = 'N';
01254
01255 index = letters[0] - 22;
01256 ltr2_low_value = UPS_Constant_Table[index].ltr2_low_value;
01257 ltr2_high_value = UPS_Constant_Table[index].ltr2_high_value;
01258 ltr3_high_value = UPS_Constant_Table[index].ltr3_high_value;
01259 false_easting = UPS_Constant_Table[index].false_easting;
01260 false_northing = UPS_Constant_Table[index].false_northing;
01261 }
01262 else if ((letters[0] == LETTER_A) || (letters[0] == LETTER_B))
01263 {
01264 hemisphere = 'S';
01265
01266 ltr2_low_value = UPS_Constant_Table[letters[0]].ltr2_low_value;
01267 ltr2_high_value = UPS_Constant_Table[letters[0]].ltr2_high_value;
01268 ltr3_high_value = UPS_Constant_Table[letters[0]].ltr3_high_value;
01269 false_easting = UPS_Constant_Table[letters[0]].false_easting;
01270 false_northing = UPS_Constant_Table[letters[0]].false_northing;
01271 }
01272 else
01273 throw CoordinateConversionException( ErrorMessages::usngString );
01274
01275
01276
01277
01278 if ((letters[1] < ltr2_low_value) || (letters[1] > ltr2_high_value) ||
01279 ((letters[1] == LETTER_D) || (letters[1] == LETTER_E) ||
01280 (letters[1] == LETTER_M) || (letters[1] == LETTER_N) ||
01281 (letters[1] == LETTER_V) || (letters[1] == LETTER_W)) ||
01282 (letters[2] > ltr3_high_value))
01283 throw CoordinateConversionException( ErrorMessages::usngString );
01284
01285 grid_northing = (double)letters[2] * ONEHT + false_northing;
01286 if (letters[2] > LETTER_I)
01287 grid_northing = grid_northing - ONEHT;
01288
01289 if (letters[2] > LETTER_O)
01290 grid_northing = grid_northing - ONEHT;
01291
01292 grid_easting = (double)((letters[1]) - ltr2_low_value) * ONEHT +false_easting;
01293 if (ltr2_low_value != LETTER_A)
01294 {
01295 if (letters[1] > LETTER_L)
01296 grid_easting = grid_easting - 300000.0;
01297
01298 if (letters[1] > LETTER_U)
01299 grid_easting = grid_easting - 200000.0;
01300 }
01301 else
01302 {
01303 if (letters[1] > LETTER_C)
01304 grid_easting = grid_easting - 200000.0;
01305
01306 if (letters[1] > LETTER_I)
01307 grid_easting = grid_easting - ONEHT;
01308
01309 if (letters[1] > LETTER_L)
01310 grid_easting = grid_easting - 300000.0;
01311 }
01312
01313 easting = grid_easting + easting;
01314 northing = grid_northing + northing;
01315
01316 return new UPSCoordinates(
01317 CoordinateType::universalPolarStereographic,
01318 hemisphere, easting, northing);
01319 }
01320
01321
01322 void USNG::getGridValues(
01323 long zone,
01324 long* ltr2_low_value,
01325 long* ltr2_high_value,
01326 double *pattern_offset )
01327 {
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341 long set_number;
01342
01343 set_number = zone % 6;
01344
01345 if (!set_number)
01346 set_number = 6;
01347
01348 if ((set_number == 1) || (set_number == 4))
01349 {
01350 *ltr2_low_value = LETTER_A;
01351 *ltr2_high_value = LETTER_H;
01352 }
01353 else if ((set_number == 2) || (set_number == 5))
01354 {
01355 *ltr2_low_value = LETTER_J;
01356 *ltr2_high_value = LETTER_R;
01357 }
01358 else if ((set_number == 3) || (set_number == 6))
01359 {
01360 *ltr2_low_value = LETTER_S;
01361 *ltr2_high_value = LETTER_Z;
01362 }
01363
01364
01365 if ((set_number % 2) == 0)
01366 *pattern_offset = 500000.0;
01367 else
01368 *pattern_offset = 0.0;
01369 }
01370
01371
01372 void USNG::getLatitudeBandMinNorthing(
01373 long letter,
01374 double* min_northing,
01375 double* northing_offset )
01376 {
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387 if ((letter >= LETTER_C) && (letter <= LETTER_H))
01388 {
01389 *min_northing = Latitude_Band_Table[letter-2].min_northing;
01390 *northing_offset = Latitude_Band_Table[letter-2].northing_offset;
01391 }
01392 else if ((letter >= LETTER_J) && (letter <= LETTER_N))
01393 {
01394 *min_northing = Latitude_Band_Table[letter-3].min_northing;
01395 *northing_offset = Latitude_Band_Table[letter-3].northing_offset;
01396 }
01397 else if ((letter >= LETTER_P) && (letter <= LETTER_X))
01398 {
01399 *min_northing = Latitude_Band_Table[letter-4].min_northing;
01400 *northing_offset = Latitude_Band_Table[letter-4].northing_offset;
01401 }
01402 else
01403 throw CoordinateConversionException( ErrorMessages::usngString );
01404 }
01405
01406
01407 void USNG::getLatitudeRange( long letter, double* north, double* south )
01408 {
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419 if ((letter >= LETTER_C) && (letter <= LETTER_H))
01420 {
01421 *north = Latitude_Band_Table[letter-2].north * PI_OVER_180;
01422 *south = Latitude_Band_Table[letter-2].south * PI_OVER_180;
01423 }
01424 else if ((letter >= LETTER_J) && (letter <= LETTER_N))
01425 {
01426 *north = Latitude_Band_Table[letter-3].north * PI_OVER_180;
01427 *south = Latitude_Band_Table[letter-3].south * PI_OVER_180;
01428 }
01429 else if ((letter >= LETTER_P) && (letter <= LETTER_X))
01430 {
01431 *north = Latitude_Band_Table[letter-4].north * PI_OVER_180;
01432 *south = Latitude_Band_Table[letter-4].south * PI_OVER_180;
01433 }
01434 else
01435 throw CoordinateConversionException( ErrorMessages::usngString );
01436 }
01437
01438
01439 void USNG::getLatitudeLetter( double latitude, int* letter )
01440 {
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450 long band = 0;
01451
01452 if (latitude >= _72 && latitude < _84_5)
01453 *letter = LETTER_X;
01454 else if (latitude > -_80_5 && latitude < _72)
01455 {
01456 band = (long)(((latitude + _80) / _8) + 1.0e-12);
01457 if(band < 0)
01458 band = 0;
01459 *letter = Latitude_Band_Table[band].letter;
01460 }
01461 else
01462 throw CoordinateConversionException( ErrorMessages::latitude );
01463 }
01464
01465