123 using namespace MSP::CCS;
131 #define EPSILON 1.75e-7
160 #define ONEHT 100000.e0
161 #define TWOMIL 2000000.e0
164 #define PI 3.14159265358979323e0
165 #define PI_OVER_2 (PI / 2.0e0)
166 #define PI_OVER_180 (PI / 180.0e0)
168 #define MIN_EASTING 100000.0
169 #define MAX_EASTING 900000.0
170 #define MIN_NORTHING 0.0
171 #define MAX_NORTHING 10000000.0
172 #define MAX_PRECISION 5
173 #define MIN_MGRS_NON_POLAR_LAT (-80.0 * ( PI / 180.0 ))
174 #define MAX_MGRS_NON_POLAR_LAT ( 84.0 * ( PI / 180.0 ))
176 #define MIN_EAST_NORTH 0.0
177 #define MAX_EAST_NORTH 3999999.0
179 #define _6 (6.0 * (PI / 180.0))
180 #define _8 (8.0 * (PI / 180.0))
181 #define _72 (72.0 * (PI / 180.0))
182 #define _80 (80.0 * (PI / 180.0))
183 #define _80_5 (80.5 * (PI / 180.0))
184 #define _84_5 (84.5 * (PI / 180.0))
186 #define _500000 500000.0
194 #define CLARKE_1866 "CC"
195 #define CLARKE_1880 "CD"
196 #define BESSEL_1841 "BR"
197 #define BESSEL_1841_NAMIBIA "BN"
199 #define EPSILON2 4.99e-4
211 {
LETTER_C, 1100000.0, -72.0, -80.5, 0.0},
212 {
LETTER_D, 2000000.0, -64.0, -72.0, 2000000.0},
213 {
LETTER_E, 2800000.0, -56.0, -64.0, 2000000.0},
214 {
LETTER_F, 3700000.0, -48.0, -56.0, 2000000.0},
215 {
LETTER_G, 4600000.0, -40.0, -48.0, 4000000.0},
216 {
LETTER_H, 5500000.0, -32.0, -40.0, 4000000.0},
217 {
LETTER_J, 6400000.0, -24.0, -32.0, 6000000.0},
218 {
LETTER_K, 7300000.0, -16.0, -24.0, 6000000.0},
219 {
LETTER_L, 8200000.0, -8.0, -16.0, 8000000.0},
220 {
LETTER_M, 9100000.0, 0.0, -8.0, 8000000.0},
222 {
LETTER_P, 800000.0, 16.0, 8.0, 0.0},
223 {
LETTER_Q, 1700000.0, 24.0, 16.0, 0.0},
224 {
LETTER_R, 2600000.0, 32.0, 24.0, 2000000.0},
225 {
LETTER_S, 3500000.0, 40.0, 32.0, 2000000.0},
226 {
LETTER_T, 4400000.0, 48.0, 40.0, 4000000.0},
227 {
LETTER_U, 5300000.0, 56.0, 48.0, 4000000.0},
228 {
LETTER_V, 6200000.0, 64.0, 56.0, 6000000.0},
229 {
LETTER_W, 7000000.0, 72.0, 64.0, 6000000.0},
230 {
LETTER_X, 7900000.0, 84.5, 72.0, 6000000.0}};
258 double scale = 1.0e5;
316 char alphabet[] =
"ABCDEFGHIJKLMNOPQRSTUVWXYZ";
320 i = sprintf (MGRSString+i,
"%2.2ld",zone);
322 strncpy(MGRSString,
" ", 2);
325 MGRSString[i++] = alphabet[letters[j]];
329 easting = fmod (easting, 100000.0);
330 if (easting >= 99999.5)
332 east = (long)((easting+
EPSILON2) /divisor);
333 i += sprintf (MGRSString+i,
"%*.*ld", precision, precision, east);
334 northing = fmod (northing, 100000.0);
335 if (northing >= 99999.5)
337 north = (long)((northing+
EPSILON2) /divisor);
338 i += sprintf (MGRSString+i,
"%*.*ld", precision, precision, north);
368 char tempMGRSString[100+1];
369 while( MGRSString[i] !=
'\0' || j == 100 )
371 if( MGRSString[i] !=
' ' )
374 if (!isdigit(MGRSString[i]) && !isalpha(MGRSString[i]) )
376 tempMGRSString[j] = MGRSString[i];
381 tempMGRSString[j] =
'\0';
386 while (tempMGRSString[i] ==
' ')
389 while (isdigit(tempMGRSString[i]))
397 strncpy (zone_string, tempMGRSString+j, 2);
399 sscanf (zone_string,
"%ld", zone);
400 if ((*zone < 1) || (*zone > 60))
409 while (isalpha(tempMGRSString[i]))
412 if (num_letters == 3)
415 letters[0] = (toupper(tempMGRSString[j]) - (long)
'A');
418 letters[1] = (toupper(tempMGRSString[j+1]) - (long)
'A');
421 letters[2] = (toupper(tempMGRSString[j+2]) - (long)
'A');
428 while (isdigit(tempMGRSString[i]))
431 if ((num_digits <= 10) && (num_digits%2 == 0))
435 char north_string[6];
444 strncpy (east_string, tempMGRSString+j, n);
446 sscanf (east_string,
"%ld", &east);
447 strncpy (north_string, tempMGRSString+j+n, n);
449 sscanf (north_string,
"%ld", &north);
452 *easting = east * multiplier;
453 *northing = north * multiplier;
473 double ellipsoidSemiMajorAxis,
474 double ellipsoidFlattening,
475 char* ellipsoidCode ) :
490 double inv_f = 1 / ellipsoidFlattening;
492 if (ellipsoidSemiMajorAxis <= 0.0)
496 if ((inv_f < 250) || (inv_f > 350))
504 strncpy (MGRSEllipsoidCode, ellipsoidCode, 2);
505 MGRSEllipsoidCode[2] =
'\0';
515 ups =
new UPS( *( m.ups ) );
516 utm =
new UTM( *( m.utm ) );
520 strcpy( MGRSEllipsoidCode, m.MGRSEllipsoidCode );
538 ups->operator=( *m.ups );
539 utm->operator=( *m.utm );
543 strcpy( MGRSEllipsoidCode, m.MGRSEllipsoidCode );
587 double latitude = geodeticCoordinates->
latitude();
588 double longitude = geodeticCoordinates->
longitude();
610 mgrsorUSNGCoordinates = fromUTM(
611 utmCoordinates, longitude, latitude, precision );
612 delete utmCoordinates;
618 mgrsorUSNGCoordinates = fromUPS( upsCoordinates, precision );
619 delete upsCoordinates;
624 delete utmCoordinates;
625 delete upsCoordinates;
629 return mgrsorUSNGCoordinates;
651 double mgrs_northing;
658 mgrsorUSNGCoordinates->
MGRSString(), &zone, letters,
659 &mgrs_easting, &mgrs_northing, &precision );
665 utmCoordinates = toUTM(
666 zone, letters, mgrs_easting, mgrs_northing, precision );
673 delete utmCoordinates;
678 upsCoordinates = toUPS( letters, mgrs_easting, mgrs_northing );
680 delete upsCoordinates;
686 delete utmCoordinates;
687 delete upsCoordinates;
691 return geodeticCoordinates;
713 long zone = utmCoordinates->
zone();
714 char hemisphere = utmCoordinates->
hemisphere();
715 double easting = utmCoordinates->
easting();
716 double northing = utmCoordinates->
northing();
722 if ((zone < 1) || (zone > 60))
724 if ((hemisphere !=
'S') && (hemisphere !=
'N'))
740 double latitude = geodeticCoordinates->
latitude();
744 mgrsorUSNGCoordinates = fromUTM(
746 geodeticCoordinates->
longitude(), latitude, precision);
750 mgrsorUSNGCoordinates = fromUPS( upsCoordinates, precision );
755 delete upsCoordinates;
756 delete geodeticCoordinates;
760 delete upsCoordinates;
761 delete geodeticCoordinates;
763 return mgrsorUSNGCoordinates;
785 double mgrs_easting, mgrs_northing;
794 mgrsorUSNGCoordinates->
MGRSString(), &zone, letters,
795 &mgrs_easting, &mgrs_northing, &precision );
798 utmCoordinates = toUTM(
799 zone, letters, mgrs_easting, mgrs_northing, precision );
806 upsCoordinates = toUPS( letters, mgrs_easting, mgrs_northing );
813 delete utmCoordinates;
814 delete upsCoordinates;
815 delete geodeticCoordinates;
819 delete upsCoordinates;
820 delete geodeticCoordinates;
822 return utmCoordinates;
845 char hemisphere = upsCoordinates->
hemisphere();
846 double easting = upsCoordinates->
easting();
847 double northing = upsCoordinates->
northing();
853 if ((hemisphere !=
'N') && (hemisphere !=
'S'))
869 double latitude = geodeticCoordinates->
latitude();
873 mgrsorUSNGCoordinates = fromUPS( upsCoordinates, precision );
877 double longitude = geodeticCoordinates->
longitude();
878 mgrsorUSNGCoordinates = fromUTM(
879 utmCoordinates, longitude, latitude, precision );
884 delete utmCoordinates;
885 delete geodeticCoordinates;
889 delete utmCoordinates;
890 delete geodeticCoordinates;
892 return mgrsorUSNGCoordinates;
915 double mgrs_northing;
925 &zone, letters, &mgrs_easting, &mgrs_northing, &precision );
929 upsCoordinates = toUPS( letters, mgrs_easting, mgrs_northing );
936 utmCoordinates = toUTM(
937 zone, letters, mgrs_easting, mgrs_northing, precision );
949 delete utmCoordinates;
950 delete upsCoordinates;
951 delete geodeticCoordinates;
955 delete utmCoordinates;
956 delete geodeticCoordinates;
958 return upsCoordinates;
982 double pattern_offset;
983 double grid_northing;
985 long ltr2_high_value;
991 long zone = utmCoordinates->
zone();
992 char hemisphere = utmCoordinates->
hemisphere();
993 double easting = utmCoordinates->
easting();
994 double northing = utmCoordinates->
northing();
996 getLatitudeLetter( latitude, &letters[0] );
1003 natural_zone = (long)(31 + ((longitude+pad) /
_6));
1007 natural_zone = (long)(((longitude+pad) /
_6) - 29);
1010 if (natural_zone > 60)
1012 if (zone != natural_zone)
1019 utmOverride.convertFromGeodetic( &geodeticCoordinates );
1021 zone = utmCoordinatesOverride->
zone();
1022 hemisphere = utmCoordinatesOverride->
hemisphere();
1023 easting = utmCoordinatesOverride->
easting();
1024 northing = utmCoordinatesOverride->
northing();
1026 delete utmCoordinatesOverride;
1027 utmCoordinatesOverride = 0;
1033 if ((zone == 31) && (easting >=
_500000))
1038 if ((zone == 32) && (easting <
_500000))
1040 else if (((zone == 32) && (easting >=
_500000)) ||
1041 ((zone == 34) && (easting <
_500000)))
1043 else if (((zone == 34) && (easting >=
_500000)) ||
1044 ((zone == 36) && (easting <
_500000)))
1046 else if ((zone == 36) && (easting >=
_500000))
1057 utmOverride.convertFromGeodetic( &geodeticCoordinates );
1059 zone = utmCoordinatesOverride->
zone();
1060 hemisphere = utmCoordinatesOverride->
hemisphere();
1061 easting = utmCoordinatesOverride->
easting();
1062 northing = utmCoordinatesOverride->
northing();
1064 delete utmCoordinatesOverride;
1065 utmCoordinatesOverride = 0;
1070 easting = ( long )( (easting +
EPSILON2) /divisor ) * divisor;
1071 northing = ( long )( (northing+
EPSILON2) /divisor ) * divisor;
1073 if( latitude <= 0.0 && northing == 1.0e7 )
1079 getGridValues( zone, <r2_low_value, <r2_high_value, &pattern_offset );
1081 grid_northing = northing;
1083 while (grid_northing >=
TWOMIL)
1085 grid_northing = grid_northing -
TWOMIL;
1087 grid_northing = grid_northing + pattern_offset;
1088 if(grid_northing >=
TWOMIL)
1089 grid_northing = grid_northing -
TWOMIL;
1091 letters[2] = (long)(grid_northing /
ONEHT);
1093 letters[2] = letters[2] + 1;
1096 letters[2] = letters[2] + 1;
1098 letters[1] = ltr2_low_value + ((long)(easting /
ONEHT) - 1);
1100 letters[1] = letters[1] + 1;
1102 makeMGRSString( MGRSString, zone, letters, easting, northing, precision );
1131 double min_northing;
1132 double northing_offset;
1133 long ltr2_low_value;
1134 long ltr2_high_value;
1135 double pattern_offset;
1136 double grid_easting;
1137 double grid_northing;
1138 double temp_grid_northing = 0.0;
1139 double fabs_grid_northing = 0.0;
1140 double latitude = 0.0;
1141 double longitude = 0.0;
1142 double divisor = 1.0;
1145 if((letters[0] ==
LETTER_X) && ((zone == 32) || (zone == 34) || (zone == 36)))
1147 else if ((letters[0] ==
LETTER_V) && (zone == 31) && (letters[1] >
LETTER_D))
1156 getGridValues(zone, <r2_low_value, <r2_high_value, &pattern_offset);
1161 if((letters[1] < ltr2_low_value) ||
1162 (letters[1] > ltr2_high_value) ||
1166 grid_easting = (double)((letters[1]) - ltr2_low_value + 1) *
ONEHT;
1168 grid_easting = grid_easting -
ONEHT;
1170 double row_letter_northing = (double)(letters[2]) *
ONEHT;
1172 row_letter_northing = row_letter_northing -
ONEHT;
1175 row_letter_northing = row_letter_northing -
ONEHT;
1177 if (row_letter_northing >= TWOMIL)
1178 row_letter_northing = row_letter_northing -
TWOMIL;
1180 getLatitudeBandMinNorthing(letters[0], &min_northing, &northing_offset);
1182 grid_northing = row_letter_northing - pattern_offset;
1183 if(grid_northing < 0)
1186 grid_northing += northing_offset;
1188 if(grid_northing < min_northing)
1191 easting = grid_easting + easting;
1192 northing = grid_northing + northing;
1196 zone, hemisphere, easting, northing );
1206 delete utmCoordinates;
1210 double latitude = geodeticCoordinates->
latitude();
1212 delete geodeticCoordinates;
1213 geodeticCoordinates = 0;
1217 if( ! inLatitudeRange(letters[0], latitude,
PI_OVER_180/divisor) )
1220 long prevBand = letters[0] - 1;
1221 long nextBand = letters[0] + 1;
1224 prevBand = letters[0];
1227 nextBand = letters[0];
1235 if(inLatitudeRange( prevBand, latitude,
PI_OVER_180/divisor ) ||
1236 inLatitudeRange( nextBand, latitude,
PI_OVER_180/divisor ) )
1248 return utmCoordinates;
1268 double false_easting;
1269 double false_northing;
1270 double grid_easting;
1271 double grid_northing;
1272 long ltr2_low_value;
1276 char MGRSString[21];
1278 char hemisphere = upsCoordinates->
hemisphere();
1279 double easting = upsCoordinates->
easting();
1280 double northing = upsCoordinates->
northing();
1284 easting = (long)((easting +
EPSILON2) /divisor) * divisor;
1285 northing = (long)((northing+
EPSILON2) /divisor) * divisor;
1287 if (hemisphere ==
'N')
1289 if (easting >= TWOMIL)
1294 index = letters[0] - 22;
1301 if (easting >= TWOMIL)
1311 grid_northing = northing;
1312 grid_northing = grid_northing - false_northing;
1313 letters[2] = (long)(grid_northing / ONEHT);
1316 letters[2] = letters[2] + 1;
1319 letters[2] = letters[2] + 1;
1321 grid_easting = easting;
1322 grid_easting = grid_easting - false_easting;
1323 letters[1] = ltr2_low_value + ((long)(grid_easting / ONEHT));
1325 if (easting < TWOMIL)
1328 letters[1] = letters[1] + 3;
1331 letters[1] = letters[1] + 2;
1336 letters[1] = letters[1] + 2;
1339 letters[1] = letters[1] + 1;
1342 letters[1] = letters[1] + 3;
1345 makeMGRSString( MGRSString, 0, letters, easting, northing, precision );
1354 long letters[MGRS_LETTERS],
1370 long ltr2_high_value;
1371 long ltr3_high_value;
1372 long ltr2_low_value;
1373 double false_easting;
1374 double false_northing;
1375 double grid_easting;
1376 double grid_northing;
1384 index = letters[0] - 22;
1407 if ((letters[1] < ltr2_low_value) || (letters[1] > ltr2_high_value) ||
1411 ( letters[2] > ltr3_high_value))
1414 grid_northing = (double)letters[2] * ONEHT + false_northing;
1416 grid_northing = grid_northing -
ONEHT;
1419 grid_northing = grid_northing -
ONEHT;
1421 grid_easting = (double)((letters[1]) - ltr2_low_value) *ONEHT + false_easting;
1425 grid_easting = grid_easting - 300000.0;
1428 grid_easting = grid_easting - 200000.0;
1433 grid_easting = grid_easting - 200000.0;
1436 grid_easting = grid_easting -
ONEHT;
1439 grid_easting = grid_easting - 300000.0;
1442 easting = grid_easting + easting;
1443 northing = grid_northing + northing;
1450 void MGRS::getGridValues(
1452 long* ltr2_low_value,
1453 long* ltr2_high_value,
1454 double* pattern_offset )
1472 set_number = zone % 6;
1485 if ((set_number == 1) || (set_number == 4))
1490 else if ((set_number == 2) || (set_number == 5))
1495 else if ((set_number == 3) || (set_number == 6))
1504 if ((set_number % 2) == 0)
1505 *pattern_offset = 500000.0;
1507 *pattern_offset = 0.0;
1511 if ((set_number % 2) == 0)
1512 *pattern_offset = 1500000.0;
1514 *pattern_offset = 1000000.00;
1519 void MGRS::getLatitudeBandMinNorthing(
1521 double* min_northing,
1522 double* northing_offset )
1554 bool MGRS::inLatitudeRange(
long letter,
double latitude,
double border )
1565 bool result =
false;
1587 if( ((south - border) <= latitude) && (latitude <= (north + border)) )
1594 void MGRS::getLatitudeLetter(
double latitude,
int* letter )
1607 if (latitude >=
_72 && latitude <
_84_5)
1609 else if (latitude > -
_80_5 && latitude <
_72)
1611 band = (long)(((latitude +
_80) /
_8) + 1.0e-12);