115 using namespace MSP::CCS;
122 #define EPSILON 1.75e-7
152 const double PI = 3.14159265358979323e0;
167 const double _6 = (6.0 * (
PI / 180.0));
168 const double _8 = (8.0 * (
PI / 180.0));
169 const double _72 = (72.0 * (
PI / 180.0));
170 const double _80 = (80.0 * (
PI / 180.0));
174 #define _500000 500000.0
182 double northing_offset;
186 {{
LETTER_C, 1100000.0, -72.0, -80.5, 0.0},
187 {
LETTER_D, 2000000.0, -64.0, -72.0, 2000000.0},
188 {
LETTER_E, 2800000.0, -56.0, -64.0, 2000000.0},
189 {
LETTER_F, 3700000.0, -48.0, -56.0, 2000000.0},
190 {
LETTER_G, 4600000.0, -40.0, -48.0, 4000000.0},
191 {
LETTER_H, 5500000.0, -32.0, -40.0, 4000000.0},
192 {
LETTER_J, 6400000.0, -24.0, -32.0, 6000000.0},
193 {
LETTER_K, 7300000.0, -16.0, -24.0, 6000000.0},
194 {
LETTER_L, 8200000.0, -8.0, -16.0, 8000000.0},
195 {
LETTER_M, 9100000.0, 0.0, -8.0, 8000000.0},
197 {
LETTER_P, 800000.0, 16.0, 8.0, 0.0},
198 {
LETTER_Q, 1700000.0, 24.0, 16.0, 0.0},
199 {
LETTER_R, 2600000.0, 32.0, 24.0, 2000000.0},
200 {
LETTER_S, 3500000.0, 40.0, 32.0, 2000000.0},
201 {
LETTER_T, 4400000.0, 48.0, 40.0, 4000000.0},
202 {
LETTER_U, 5300000.0, 56.0, 48.0, 4000000.0},
203 {
LETTER_V, 6200000.0, 64.0, 56.0, 6000000.0},
204 {
LETTER_W, 7000000.0, 72.0, 64.0, 6000000.0},
205 {
LETTER_X, 7900000.0, 84.5, 72.0, 6000000.0}};
212 long ltr2_high_value;
213 long ltr3_high_value;
214 double false_easting;
215 double false_northing;
255 char alphabet[] =
"ABCDEFGHIJKLMNOPQRSTUVWXYZ";
259 i = sprintf (USNGString+i,
"%2.2ld",zone);
261 strncpy(USNGString,
" ", 2);
264 USNGString[i++] = alphabet[letters[j]];
265 divisor = pow (10.0, (5.0 - precision));
266 easting = fmod (easting, 100000.0);
267 if (easting >= 99999.5)
269 east = (long)(easting/divisor);
270 i += sprintf (USNGString+i,
"%*.*ld", precision, precision, east);
271 northing = fmod (northing, 100000.0);
272 if (northing >= 99999.5)
274 north = (long)(northing/divisor);
275 i += sprintf (USNGString+i,
"%*.*ld", precision, precision, north);
304 while (USNGString[i] ==
' ')
307 while (isdigit(USNGString[i]))
315 strncpy (zone_string, USNGString+j, 2);
317 sscanf (zone_string,
"%ld", zone);
318 if ((*zone < 1) || (*zone > 60))
327 while (isalpha(USNGString[i]))
330 if (num_letters == 3)
333 letters[0] = (toupper(USNGString[j]) - (long)
'A');
336 letters[1] = (toupper(USNGString[j+1]) - (long)
'A');
339 letters[2] = (toupper(USNGString[j+2]) - (long)
'A');
346 while (isdigit(USNGString[i]))
349 if ((num_digits <= 10) && (num_digits%2 == 0))
353 char north_string[6];
362 strncpy (east_string, USNGString+j, n);
364 sscanf (east_string,
"%ld", &east);
365 strncpy (north_string, USNGString+j+n, n);
367 sscanf (north_string,
"%ld", &north);
368 multiplier = pow (10.0, 5.0 - n);
369 *easting = east * multiplier;
370 *northing = north * multiplier;
388 USNG::USNG(
double ellipsoidSemiMajorAxis,
double ellipsoidFlattening,
char* ellipsoidCode ) :
403 double inv_f = 1 / ellipsoidFlattening;
405 if (ellipsoidSemiMajorAxis <= 0.0)
409 if ((inv_f < 250) || (inv_f > 350))
417 strncpy (USNGEllipsoidCode, ellipsoidCode, 2);
418 USNGEllipsoidCode[2] =
'\0';
428 ups =
new UPS( *( u.ups ) );
429 utm =
new UTM( *( u.utm ) );
433 strcpy( USNGEllipsoidCode, u.USNGEllipsoidCode );
451 ups->operator=( *u.ups );
452 utm->operator=( *u.utm );
456 strcpy( USNGEllipsoidCode, u.USNGEllipsoidCode );
500 double latitude = geodeticCoordinates->
latitude();
501 double longitude = geodeticCoordinates->
longitude();
507 if ((longitude < -
PI) || (longitude > (2*
PI)))
523 mgrsorUSNGCoordinates = fromUTM(
524 utmCoordinates, longitude, latitude, precision );
525 delete utmCoordinates;
530 mgrsorUSNGCoordinates = fromUPS( upsCoordinates, precision );
531 delete upsCoordinates;
535 delete utmCoordinates;
536 delete upsCoordinates;
540 return mgrsorUSNGCoordinates;
562 double usng_northing;
570 &zone, letters, &usng_easting, &usng_northing, &precision );
576 utmCoordinates = toUTM(
577 zone, letters, usng_easting, usng_northing, precision );
579 delete utmCoordinates;
583 upsCoordinates = toUPS( letters, usng_easting, usng_northing );
585 delete upsCoordinates;
590 delete utmCoordinates;
591 delete upsCoordinates;
595 return geodeticCoordinates;
616 long zone = utmCoordinates->
zone();
617 char hemisphere = utmCoordinates->
hemisphere();
618 double easting = utmCoordinates->
easting();
619 double northing = utmCoordinates->
northing();
621 if ((zone < 1) || (zone > 60))
623 if ((hemisphere !=
'S') && (hemisphere !=
'N'))
640 double latitude = geodeticCoordinates->
latitude();
646 mgrsorUSNGCoordinates = fromUTM(
647 utmCoordinates, geodeticCoordinates->
longitude(),
648 latitude, precision );
652 mgrsorUSNGCoordinates = fromUPS( upsCoordinates, precision );
657 delete geodeticCoordinates;
658 delete upsCoordinates;
662 delete geodeticCoordinates;
663 delete upsCoordinates;
665 return mgrsorUSNGCoordinates;
687 double usng_easting, usng_northing;
694 &zone, letters, &usng_easting, &usng_northing, &precision );
699 utmCoordinates = toUTM(
700 zone, letters, usng_easting, usng_northing, precision );
706 upsCoordinates = toUPS( letters, usng_easting, usng_northing );
713 delete utmCoordinates;
714 delete geodeticCoordinates;
715 delete upsCoordinates;
719 delete geodeticCoordinates;
720 delete upsCoordinates;
722 return utmCoordinates;
745 char hemisphere = upsCoordinates->
hemisphere();
746 double easting = upsCoordinates->
easting();
747 double northing = upsCoordinates->
northing();
749 if ((hemisphere !=
'N') && (hemisphere !=
'S'))
766 double latitude = geodeticCoordinates->
latitude();
772 mgrsorUSNGCoordinates = fromUPS( upsCoordinates, precision );
776 double longitude = geodeticCoordinates->
longitude();
777 mgrsorUSNGCoordinates = fromUTM(
778 utmCoordinates, longitude, latitude, precision );
783 delete geodeticCoordinates;
784 delete utmCoordinates;
788 delete geodeticCoordinates;
789 delete utmCoordinates;
791 return mgrsorUSNGCoordinates;
814 double usng_northing;
821 &zone, letters, &usng_easting, &usng_northing, &precision );
826 upsCoordinates = toUPS( letters, usng_easting, usng_northing );
832 utmCoordinates = toUTM(
833 zone, letters, usng_easting, usng_northing, precision );
840 delete geodeticCoordinates;
841 delete utmCoordinates;
842 delete upsCoordinates;
846 delete geodeticCoordinates;
847 delete utmCoordinates;
849 return upsCoordinates;
871 double pattern_offset;
872 double grid_northing;
874 long ltr2_high_value;
880 long zone = utmCoordinates->
zone();
881 char hemisphere = utmCoordinates->
hemisphere();
882 double easting = utmCoordinates->
easting();
883 double northing = utmCoordinates->
northing();
885 getLatitudeLetter( latitude, &letters[0] );
890 natural_zone = (long)(31 + ((longitude) /
_6));
892 natural_zone = (long)(((longitude) /
_6) - 29);
894 if (natural_zone > 60)
897 if (zone != natural_zone)
902 UTMCoordinates* utmCoordinatesOverride = utmOverride.convertFromGeodetic(
903 &geodeticCoordinates );
905 zone = utmCoordinatesOverride->
zone();
906 hemisphere = utmCoordinatesOverride->
hemisphere();
907 easting = utmCoordinatesOverride->
easting();
908 northing = utmCoordinatesOverride->
northing();
910 delete utmCoordinatesOverride;
911 utmCoordinatesOverride = 0;
917 if ((zone == 31) && (easting >=
_500000))
922 if ((zone == 32) && (easting <
_500000))
924 else if (((zone == 32) && (easting >=
_500000)) ||
925 ((zone == 34) && (easting <
_500000)))
927 else if (((zone == 34) && (easting >=
_500000)) ||
928 ((zone == 36) && (easting <
_500000)))
930 else if ((zone == 36) && (easting >=
_500000))
940 utmOverride.convertFromGeodetic( &geodeticCoordinates );
942 zone = utmCoordinatesOverride->
zone();
943 hemisphere = utmCoordinatesOverride->
hemisphere();
944 easting = utmCoordinatesOverride->
easting();
945 northing = utmCoordinatesOverride->
northing();
947 delete utmCoordinatesOverride;
948 utmCoordinatesOverride = 0;
952 double divisor = pow (10.0, (5.0 - precision));
953 easting = (long)(easting/divisor) * divisor;
954 northing = (long)(northing/divisor) * divisor;
956 if( latitude <= 0.0 && northing == 1.0e7)
962 getGridValues( zone, <r2_low_value, <r2_high_value, &pattern_offset );
964 grid_northing = northing;
966 while (grid_northing >=
TWOMIL)
968 grid_northing = grid_northing -
TWOMIL;
970 grid_northing = grid_northing + pattern_offset;
971 if(grid_northing >=
TWOMIL)
972 grid_northing = grid_northing -
TWOMIL;
974 letters[2] = (long)(grid_northing /
ONEHT);
976 letters[2] = letters[2] + 1;
979 letters[2] = letters[2] + 1;
981 letters[1] = ltr2_low_value + ((long)(easting /
ONEHT) -1);
983 letters[1] = letters[1] + 1;
985 makeUSNGString( USNGString, zone, letters, easting, northing, precision );
1012 double min_northing;
1013 double northing_offset;
1014 long ltr2_low_value;
1015 long ltr2_high_value;
1016 double pattern_offset;
1017 double upper_lat_limit;
1018 double lower_lat_limit;
1019 double grid_easting;
1020 double grid_northing;
1021 double temp_grid_northing = 0.0;
1022 double fabs_grid_northing = 0.0;
1023 double latitude = 0.0;
1024 double longitude = 0.0;
1025 double divisor = 1.0;
1028 if((letters[0] ==
LETTER_X) && ((zone == 32) || (zone == 34) || (zone == 36)))
1030 else if ((letters[0] ==
LETTER_V) && (zone == 31) && (letters[1] >
LETTER_D))
1039 getGridValues(zone, <r2_low_value, <r2_high_value, &pattern_offset);
1044 if((letters[1] < ltr2_low_value) ||
1045 (letters[1] > ltr2_high_value) ||
1049 grid_easting = (double)((letters[1]) - ltr2_low_value + 1) *
ONEHT;
1051 grid_easting = grid_easting -
ONEHT;
1053 double row_letter_northing = (double)(letters[2]) *
ONEHT;
1055 row_letter_northing = row_letter_northing -
ONEHT;
1058 row_letter_northing = row_letter_northing -
ONEHT;
1060 if (row_letter_northing >= TWOMIL)
1061 row_letter_northing = row_letter_northing -
TWOMIL;
1063 getLatitudeBandMinNorthing(letters[0], &min_northing, &northing_offset);
1065 grid_northing = row_letter_northing - pattern_offset;
1066 if(grid_northing < 0)
1069 grid_northing += northing_offset;
1071 if(grid_northing < min_northing)
1074 easting = grid_easting + easting;
1075 northing = grid_northing + northing;
1079 zone, hemisphere, easting, northing );
1089 delete utmCoordinates;
1093 divisor = pow (10.0, (
double)precision);
1094 getLatitudeRange(letters[0], &upper_lat_limit, &lower_lat_limit);
1096 double latitude = geodeticCoordinates->
latitude();
1098 delete geodeticCoordinates;
1100 if(!(((lower_lat_limit -
PI_OVER_180/divisor) <= latitude) &&
1101 (latitude <= (upper_lat_limit +
PI_OVER_180/divisor))))
1105 return utmCoordinates;
1125 double false_easting;
1126 double false_northing;
1127 double grid_easting;
1128 double grid_northing;
1129 long ltr2_low_value;
1133 char USNGString[21];
1135 char hemisphere = upsCoordinates->
hemisphere();
1136 double easting = upsCoordinates->
easting();
1137 double northing = upsCoordinates->
northing();
1139 divisor = pow (10.0, (5.0 - precision));
1140 easting = (long)(easting/divisor) * divisor;
1141 northing = (long)(northing/divisor) * divisor;
1143 if (hemisphere ==
'N')
1145 if (easting >= TWOMIL)
1150 index = letters[0] - 22;
1157 if (easting >= TWOMIL)
1167 grid_northing = northing;
1168 grid_northing = grid_northing - false_northing;
1169 letters[2] = (long)(grid_northing / ONEHT);
1172 letters[2] = letters[2] + 1;
1175 letters[2] = letters[2] + 1;
1177 grid_easting = easting;
1178 grid_easting = grid_easting - false_easting;
1179 letters[1] = ltr2_low_value + ((long)(grid_easting / ONEHT));
1181 if (easting < TWOMIL)
1184 letters[1] = letters[1] + 3;
1187 letters[1] = letters[1] + 2;
1192 letters[1] = letters[1] + 2;
1195 letters[1] = letters[1] + 1;
1198 letters[1] = letters[1] + 3;
1201 makeUSNGString( USNGString, 0, letters, easting, northing, precision );
1208 long letters[USNG_LETTERS],
1224 long ltr2_high_value;
1225 long ltr3_high_value;
1226 long ltr2_low_value;
1227 double false_easting;
1228 double false_northing;
1229 double grid_easting;
1230 double grid_northing;
1238 index = letters[0] - 22;
1261 if ((letters[1] < ltr2_low_value) || (letters[1] > ltr2_high_value) ||
1265 (letters[2] > ltr3_high_value))
1268 grid_northing = (double)letters[2] * ONEHT + false_northing;
1270 grid_northing = grid_northing -
ONEHT;
1273 grid_northing = grid_northing -
ONEHT;
1275 grid_easting = (double)((letters[1]) - ltr2_low_value) * ONEHT +false_easting;
1279 grid_easting = grid_easting - 300000.0;
1282 grid_easting = grid_easting - 200000.0;
1287 grid_easting = grid_easting - 200000.0;
1290 grid_easting = grid_easting -
ONEHT;
1293 grid_easting = grid_easting - 300000.0;
1296 easting = grid_easting + easting;
1297 northing = grid_northing + northing;
1301 hemisphere, easting, northing);
1305 void USNG::getGridValues(
1307 long* ltr2_low_value,
1308 long* ltr2_high_value,
1309 double *pattern_offset )
1326 set_number = zone % 6;
1331 if ((set_number == 1) || (set_number == 4))
1336 else if ((set_number == 2) || (set_number == 5))
1341 else if ((set_number == 3) || (set_number == 6))
1348 if ((set_number % 2) == 0)
1349 *pattern_offset = 500000.0;
1351 *pattern_offset = 0.0;
1355 void USNG::getLatitudeBandMinNorthing(
1357 double* min_northing,
1358 double* northing_offset )
1390 void USNG::getLatitudeRange(
long letter,
double* north,
double* south )
1422 void USNG::getLatitudeLetter(
double latitude,
int* letter )
1435 if (latitude >=
_72 && latitude <
_84_5)
1437 else if (latitude > -
_80_5 && latitude <
_72)
1439 band = (long)(((latitude +
_80) /
_8) + 1.0e-12);