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