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