1
// CLASSIFICATION: UNCLASSIFIED
3
/***************************************************************************/
4
/* RSC IDENTIFIER: USNG
8
* This component converts between geodetic coordinates (latitude and
9
* longitude) and United States National Grid (USNG) coordinates.
13
* This component checks parameters for valid values. If an invalid value
14
* is found, the error code is combined with the current error code using
15
* the bitwise or. This combining allows multiple error codes to be
16
* returned. The possible error codes are:
18
* USNG_NO_ERROR : No errors occurred in function
19
* USNG_LAT_ERROR : Latitude outside of valid range
21
* USNG_LON_ERROR : Longitude outside of valid range
22
* (-180 to 360 degrees)
23
* USNG_STR_ERROR : An USNG string error: string too long,
24
* too short, or badly formed
25
* USNG_PRECISION_ERROR : The precision must be between 0 and 5
27
* USNG_A_ERROR : Semi-major axis less than or equal to zero
28
* USNG_INV_F_ERROR : Inverse flattening outside of valid range
30
* USNG_EASTING_ERROR : Easting outside of valid range
31
* (100,000 to 900,000 meters for UTM)
32
* (0 to 4,000,000 meters for UPS)
33
* USNG_NORTHING_ERROR : Northing outside of valid range
34
* (0 to 10,000,000 meters for UTM)
35
* (0 to 4,000,000 meters for UPS)
36
* USNG_ZONE_ERROR : Zone outside of valid range (1 to 60)
37
* USNG_HEMISPHERE_ERROR : Invalid hemisphere ('N' or 'S')
41
* USNG is intended for reuse by any application that does conversions
42
* between geodetic coordinates and USNG coordinates.
46
* Further information on USNG can be found in the Reuse Manual.
48
* USNG originated from : Federal Geographic Data Committee
50
* 12201 Sunrise Valley Drive
55
* None apply to this component.
62
* USNG was tested and certified in the following environments:
64
* 1. Solaris 2.5 with GCC version 2.8.1
65
* 2. Windows XP with MS Visual C++ version 6
71
* 3-1-07 Original Code (cloned from MGRS)
75
/***************************************************************************/
86
#include "EllipsoidParameters.h"
87
#include "MGRSorUSNGCoordinates.h"
88
#include "GeodeticCoordinates.h"
89
#include "UPSCoordinates.h"
90
#include "UTMCoordinates.h"
91
#include "CoordinateConversionException.h"
92
#include "ErrorMessages.h"
93
#include "WarningMessages.h"
96
* ctype.h - Standard C character handling library
97
* math.h - Standard C math library
98
* stdio.h - Standard C input/output library
99
* string.h - Standard C string handling library
100
* UPS.h - Universal Polar Stereographic (UPS) projection
101
* UTM.h - Universal Transverse Mercator (UTM) projection
102
* USNG.h - function prototype error checking
103
* MGRSorUSNGCoordinates.h - defines mgrs coordinates
104
* GeodeticCoordinates.h - defines geodetic coordinates
105
* UPSCoordinates.h - defines ups coordinates
106
* UTMCoordinates.h - defines utm coordinates
107
* CoordinateConversionException.h - Exception handler
108
* ErrorMessages.h - Contains exception messages
109
* WarningMessages.h - Contains warning messages
113
using namespace MSP::CCS;
116
/************************************************************************/
121
#define EPSILON 1.75e-7 /* approx 1.0e-5 degrees (~1 meter) in radians */
123
const int LETTER_A = 0; /* ARRAY INDEX FOR LETTER A */
124
const int LETTER_B = 1; /* ARRAY INDEX FOR LETTER B */
125
const int LETTER_C = 2; /* ARRAY INDEX FOR LETTER C */
126
const int LETTER_D = 3; /* ARRAY INDEX FOR LETTER D */
127
const int LETTER_E = 4; /* ARRAY INDEX FOR LETTER E */
128
const int LETTER_F = 5; /* ARRAY INDEX FOR LETTER F */
129
const int LETTER_G = 6; /* ARRAY INDEX FOR LETTER G */
130
const int LETTER_H = 7; /* ARRAY INDEX FOR LETTER H */
131
const int LETTER_I = 8; /* ARRAY INDEX FOR LETTER I */
132
const int LETTER_J = 9; /* ARRAY INDEX FOR LETTER J */
133
const int LETTER_K = 10; /* ARRAY INDEX FOR LETTER K */
134
const int LETTER_L = 11; /* ARRAY INDEX FOR LETTER L */
135
const int LETTER_M = 12; /* ARRAY INDEX FOR LETTER M */
136
const int LETTER_N = 13; /* ARRAY INDEX FOR LETTER N */
137
const int LETTER_O = 14; /* ARRAY INDEX FOR LETTER O */
138
const int LETTER_P = 15; /* ARRAY INDEX FOR LETTER P */
139
const int LETTER_Q = 16; /* ARRAY INDEX FOR LETTER Q */
140
const int LETTER_R = 17; /* ARRAY INDEX FOR LETTER R */
141
const int LETTER_S = 18; /* ARRAY INDEX FOR LETTER S */
142
const int LETTER_T = 19; /* ARRAY INDEX FOR LETTER T */
143
const int LETTER_U = 20; /* ARRAY INDEX FOR LETTER U */
144
const int LETTER_V = 21; /* ARRAY INDEX FOR LETTER V */
145
const int LETTER_W = 22; /* ARRAY INDEX FOR LETTER W */
146
const int LETTER_X = 23; /* ARRAY INDEX FOR LETTER X */
147
const int LETTER_Y = 24; /* ARRAY INDEX FOR LETTER Y */
148
const int LETTER_Z = 25; /* ARRAY INDEX FOR LETTER Z */
149
const double ONEHT = 100000.e0; /* ONE HUNDRED THOUSAND */
150
const double TWOMIL = 2000000.e0; /* TWO MILLION */
151
const int TRUE = 1; /* CONSTANT VALUE FOR TRUE VALUE */
152
const int FALSE = 0; /* CONSTANT VALUE FOR FALSE VALUE */
153
const double PI = 3.14159265358979323e0; /* PI */
154
const double PI_OVER_2 = (PI / 2.0e0);
155
const double PI_OVER_180 = (PI / 180.0e0);
157
const double MIN_EASTING = 100000.0;
158
const double MAX_EASTING = 900000.0;
159
const double MIN_NORTHING = 0.0;
160
const double MAX_NORTHING = 10000000.0;
161
const int MAX_PRECISION = 5; /* Maximum precision of easting & northing */
162
const double MIN_USNG_NON_POLAR_LAT = -80.0 * ( PI / 180.0 ); /* -80 degrees in radians */
163
const double MAX_USNG_NON_POLAR_LAT = 84.0 * ( PI / 180.0 ); /* 84 degrees in radians */
165
const double MIN_EAST_NORTH = 0.0;
166
const double MAX_EAST_NORTH = 3999999.0;
168
const double _6 = (6.0 * (PI / 180.0));
169
const double _8 = (8.0 * (PI / 180.0));
170
const double _72 = (72.0 * (PI / 180.0));
171
const double _80 = (80.0 * (PI / 180.0));
172
const double _80_5 = (80.5 * (PI / 180.0));
173
const double _84_5 = (84.5 * (PI / 180.0));
175
#define _500000 500000.0
179
long letter; /* letter representing latitude band */
180
double min_northing; /* minimum northing for latitude band */
181
double north; /* upper latitude for latitude band */
182
double south; /* lower latitude for latitude band */
183
double northing_offset; /* latitude band northing offset */
186
const Latitude_Band Latitude_Band_Table[20] =
187
{{LETTER_C, 1100000.0, -72.0, -80.5, 0.0},
188
{LETTER_D, 2000000.0, -64.0, -72.0, 2000000.0},
189
{LETTER_E, 2800000.0, -56.0, -64.0, 2000000.0},
190
{LETTER_F, 3700000.0, -48.0, -56.0, 2000000.0},
191
{LETTER_G, 4600000.0, -40.0, -48.0, 4000000.0},
192
{LETTER_H, 5500000.0, -32.0, -40.0, 4000000.0},
193
{LETTER_J, 6400000.0, -24.0, -32.0, 6000000.0},
194
{LETTER_K, 7300000.0, -16.0, -24.0, 6000000.0},
195
{LETTER_L, 8200000.0, -8.0, -16.0, 8000000.0},
196
{LETTER_M, 9100000.0, 0.0, -8.0, 8000000.0},
197
{LETTER_N, 0.0, 8.0, 0.0, 0.0},
198
{LETTER_P, 800000.0, 16.0, 8.0, 0.0},
199
{LETTER_Q, 1700000.0, 24.0, 16.0, 0.0},
200
{LETTER_R, 2600000.0, 32.0, 24.0, 2000000.0},
201
{LETTER_S, 3500000.0, 40.0, 32.0, 2000000.0},
202
{LETTER_T, 4400000.0, 48.0, 40.0, 4000000.0},
203
{LETTER_U, 5300000.0, 56.0, 48.0, 4000000.0},
204
{LETTER_V, 6200000.0, 64.0, 56.0, 6000000.0},
205
{LETTER_W, 7000000.0, 72.0, 64.0, 6000000.0},
206
{LETTER_X, 7900000.0, 84.5, 72.0, 6000000.0}};
211
long letter; /* letter representing latitude band */
212
long ltr2_low_value; /* 2nd letter range - low number */
213
long ltr2_high_value; /* 2nd letter range - high number */
214
long ltr3_high_value; /* 3rd letter range - high number (UPS) */
215
double false_easting; /* False easting based on 2nd letter */
216
double false_northing; /* False northing based on 3rd letter */
219
const UPS_Constant UPS_Constant_Table[4] =
220
{{LETTER_A, LETTER_J, LETTER_Z, LETTER_Z, 800000.0, 800000.0},
221
{LETTER_B, LETTER_A, LETTER_R, LETTER_Z, 2000000.0, 800000.0},
222
{LETTER_Y, LETTER_J, LETTER_Z, LETTER_P, 800000.0, 1300000.0},
223
{LETTER_Z, LETTER_A, LETTER_J, LETTER_P, 2000000.0, 1300000.0}};
226
/************************************************************************/
231
long roundUSNG( double value )
234
* The function roundUSNG rounds the input value to the
235
* nearest integer, using the standard engineering rule.
236
* The rounded integer value is then returned.
238
* value : Value to be rounded (input)
243
double fraction = modf (value, &ivalue);
244
ival = (long)(ivalue);
245
if ((fraction > 0.5) || ((fraction == 0.5) && (ival%2 == 1)))
251
void makeUSNGString( char* USNGString, long zone, int letters[USNG_LETTERS], double easting, double northing, long precision )
254
* The function makeUSNGString constructs an USNG string
255
* from its component parts.
257
* USNGString : USNG coordinate string (output)
258
* zone : UTM Zone (input)
259
* letters : USNG coordinate string letters (input)
260
* easting : Easting value (input)
261
* northing : Northing value (input)
262
* precision : Precision level of USNG string (input)
270
char alphabet[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
274
i = sprintf (USNGString+i,"%2.2ld",zone);
276
strncpy(USNGString, " ", 2); // 2 spaces
279
USNGString[i++] = alphabet[letters[j]];
280
divisor = pow (10.0, (5.0 - precision));
281
easting = fmod (easting, 100000.0);
282
if (easting >= 99999.5)
284
east = (long)(easting/divisor);
285
i += sprintf (USNGString+i, "%*.*ld", precision, precision, east);
286
northing = fmod (northing, 100000.0);
287
if (northing >= 99999.5)
289
north = (long)(northing/divisor);
290
i += sprintf (USNGString+i, "%*.*ld", precision, precision, north);
294
void breakUSNGString( char* USNGString, long* zone, long letters[USNG_LETTERS], double* easting, double* northing, long* precision )
297
* The function breakUSNGString breaks down an USNG
298
* coordinate string into its component parts.
300
* USNG : USNG coordinate string (input)
301
* zone : UTM Zone (output)
302
* letters : USNG coordinate string letters (output)
303
* easting : Easting value (output)
304
* northing : Northing value (output)
305
* precision : Precision level of USNG string (output)
313
while (USNGString[i] == ' ')
314
i++; /* skip any leading blanks */
316
while (isdigit(USNGString[i]))
324
strncpy (zone_string, USNGString+j, 2);
326
sscanf (zone_string, "%ld", zone);
327
if ((*zone < 1) || (*zone > 60))
328
throw CoordinateConversionException( ErrorMessages::usngString );
333
throw CoordinateConversionException( ErrorMessages::usngString );
336
while (isalpha(USNGString[i]))
339
if (num_letters == 3)
342
letters[0] = (toupper(USNGString[j]) - (long)'A');
343
if ((letters[0] == LETTER_I) || (letters[0] == LETTER_O))
344
throw CoordinateConversionException( ErrorMessages::usngString );
345
letters[1] = (toupper(USNGString[j+1]) - (long)'A');
346
if ((letters[1] == LETTER_I) || (letters[1] == LETTER_O))
347
throw CoordinateConversionException( ErrorMessages::usngString );
348
letters[2] = (toupper(USNGString[j+2]) - (long)'A');
349
if ((letters[2] == LETTER_I) || (letters[2] == LETTER_O))
350
throw CoordinateConversionException( ErrorMessages::usngString );
353
throw CoordinateConversionException( ErrorMessages::usngString );
355
while (isdigit(USNGString[i]))
358
if ((num_digits <= 10) && (num_digits%2 == 0))
362
char north_string[6];
366
/* get easting & northing */
371
strncpy (east_string, USNGString+j, n);
373
sscanf (east_string, "%ld", &east);
374
strncpy (north_string, USNGString+j+n, n);
376
sscanf (north_string, "%ld", &north);
377
multiplier = pow (10.0, 5.0 - n);
378
*easting = east * multiplier;
379
*northing = north * multiplier;
388
throw CoordinateConversionException( ErrorMessages::usngString );
392
/************************************************************************/
397
USNG::USNG( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, char* ellipsoidCode ) :
403
* The constructor receives the ellipsoid parameters and sets
404
* the corresponding state variables. If any errors occur, an exception is thrown with a description
407
* ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid in meters (input)
408
* ellipsoidFlattening : Flattening of ellipsoid (input)
409
* ellipsoid_Code : 2-letter code for ellipsoid (input)
412
double inv_f = 1 / ellipsoidFlattening;
413
char errorStatus[500] = "";
415
if (ellipsoidSemiMajorAxis <= 0.0)
416
{ /* Semi-major axis must be greater than zero */
417
strcat( errorStatus, ErrorMessages::semiMajorAxis );
419
if ((inv_f < 250) || (inv_f > 350))
420
{ /* Inverse flattening must be between 250 and 350 */
421
strcat( errorStatus, ErrorMessages::ellipsoidFlattening );
424
if( strlen( errorStatus ) > 0)
425
throw CoordinateConversionException( errorStatus );
427
semiMajorAxis = ellipsoidSemiMajorAxis;
428
flattening = ellipsoidFlattening;
430
strncpy (USNGEllipsoidCode, ellipsoidCode, 2);
431
USNGEllipsoidCode[2] = '\0';
433
ups = new UPS( semiMajorAxis, flattening );
435
utm = new UTM( semiMajorAxis, flattening, 0 );
439
USNG::USNG( const USNG &u )
441
ups = new UPS( *( u.ups ) );
442
utm = new UTM( *( u.utm ) );
444
semiMajorAxis = u.semiMajorAxis;
445
flattening = u.flattening;
446
strcpy( USNGEllipsoidCode, u.USNGEllipsoidCode );
460
USNG& USNG::operator=( const USNG &u )
464
ups->operator=( *u.ups );
465
utm->operator=( *u.utm );
467
semiMajorAxis = u.semiMajorAxis;
468
flattening = u.flattening;
469
strcpy( USNGEllipsoidCode, u.USNGEllipsoidCode );
476
EllipsoidParameters* USNG::getParameters() const
479
* The function getParameters returns the current ellipsoid
482
* ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
483
* ellipsoidFlattening : Flattening of ellipsoid (output)
484
* ellipsoidCode : 2-letter code for ellipsoid (output)
487
return new EllipsoidParameters( semiMajorAxis, flattening, (char*)USNGEllipsoidCode );
491
MSP::CCS::MGRSorUSNGCoordinates* USNG::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates, long precision )
494
* The function convertFromGeodetic converts Geodetic (latitude and
495
* longitude) coordinates to an USNG coordinate string, according to the
496
* current ellipsoid parameters. If any errors occur, an exception is thrown with a description
499
* latitude : Latitude in radians (input)
500
* longitude : Longitude in radians (input)
501
* precision : Precision level of USNG string (input)
502
* USNGString : USNG coordinate string (output)
506
MGRSorUSNGCoordinates* mgrsorUSNGCoordinates = 0;
507
char errorStatus[50] = "";
509
double latitude = geodeticCoordinates->latitude();
510
double longitude = geodeticCoordinates->longitude();
512
if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
513
{ /* latitude out of range */
514
strcat( errorStatus, ErrorMessages::latitude );
516
if ((longitude < -PI) || (longitude > (2*PI)))
517
{ /* longitude out of range */
518
strcat( errorStatus, ErrorMessages::longitude );
520
if ((precision < 0) || (precision > MAX_PRECISION))
521
strcat( errorStatus, ErrorMessages::precision );
523
if( strlen( errorStatus ) > 0)
524
throw CoordinateConversionException( errorStatus );
526
// If the latitude is within the valid usng non polar range [-80, 84),
527
// convert to usng using the utm path, otherwise convert to usng using the ups path
528
if ((latitude >= MIN_USNG_NON_POLAR_LAT) && (latitude < MAX_USNG_NON_POLAR_LAT))
530
UTMCoordinates* utmCoordinates = utm->convertFromGeodetic( geodeticCoordinates );
531
mgrsorUSNGCoordinates = fromUTM( utmCoordinates, longitude, latitude, precision );
532
delete utmCoordinates;
537
UPSCoordinates* upsCoordinates = ups->convertFromGeodetic( geodeticCoordinates );
538
mgrsorUSNGCoordinates = fromUPS( upsCoordinates, precision );
539
delete upsCoordinates;
543
return mgrsorUSNGCoordinates;
547
MSP::CCS::GeodeticCoordinates* USNG::convertToGeodetic( MSP::CCS::MGRSorUSNGCoordinates* usngCoordinates )
550
* The function convertToGeodetic converts an USNG coordinate string
551
* to Geodetic (latitude and longitude) coordinates
552
* according to the current ellipsoid parameters. If any errors occur, an exception is thrown with a description
555
* USNG : USNG coordinate string (input)
556
* latitude : Latitude in radians (output)
557
* longitude : Longitude in radians (output)
562
long letters[USNG_LETTERS];
564
double usng_northing;
566
GeodeticCoordinates* geodeticCoordinates = 0;
568
breakUSNGString( usngCoordinates->MGRSString(), &zone, letters, &usng_easting, &usng_northing, &precision );
572
UTMCoordinates* utmCoordinates = toUTM( zone, letters, usng_easting, usng_northing, precision );
573
geodeticCoordinates = utm->convertToGeodetic( utmCoordinates );
574
delete utmCoordinates;
579
UPSCoordinates* upsCoordinates = toUPS( letters, usng_easting, usng_northing );
580
geodeticCoordinates = ups->convertToGeodetic( upsCoordinates );
581
delete upsCoordinates;
585
return geodeticCoordinates;
589
MSP::CCS::MGRSorUSNGCoordinates* USNG::convertFromUTM ( UTMCoordinates* utmCoordinates, long precision )
592
* The function convertFromUTM converts UTM (zone, easting, and
593
* northing) coordinates to an USNG coordinate string, according to the
594
* current ellipsoid parameters. If any errors occur, an exception is thrown with a description
597
* zone : UTM zone (input)
598
* hemisphere : North or South hemisphere (input)
599
* easting : Easting (X) in meters (input)
600
* northing : Northing (Y) in meters (input)
601
* precision : Precision level of USNG string (input)
602
* USNGString : USNG coordinate string (output)
605
char errorStatus[50] = "";
607
long zone = utmCoordinates->zone();
608
char hemisphere = utmCoordinates->hemisphere();
609
double easting = utmCoordinates->easting();
610
double northing = utmCoordinates->northing();
612
if ((zone < 1) || (zone > 60))
613
strcat( errorStatus, ErrorMessages::zone );
614
if ((hemisphere != 'S') && (hemisphere != 'N'))
615
strcat( errorStatus, ErrorMessages::hemisphere );
616
if ((easting < MIN_EASTING) || (easting > MAX_EASTING))
617
strcat( errorStatus, ErrorMessages::easting );
618
if ((northing < MIN_NORTHING) || (northing > MAX_NORTHING))
619
strcat( errorStatus, ErrorMessages::northing );
620
if ((precision < 0) || (precision > MAX_PRECISION))
621
strcat( errorStatus, ErrorMessages::precision );
623
if( strlen( errorStatus ) > 0)
624
throw CoordinateConversionException( errorStatus );
626
GeodeticCoordinates* geodeticCoordinates = utm->convertToGeodetic( utmCoordinates );
628
// If the latitude is within the valid mgrs non polar range [-80, 84),
629
// convert to mgrs using the utm path, otherwise convert to mgrs using the ups path
630
MGRSorUSNGCoordinates* mgrsorUSNGCoordinates = 0;
631
double latitude = geodeticCoordinates->latitude();
633
if ((latitude >= (MIN_USNG_NON_POLAR_LAT - EPSILON)) && (latitude < (MAX_USNG_NON_POLAR_LAT + EPSILON)))
634
mgrsorUSNGCoordinates = fromUTM( utmCoordinates, geodeticCoordinates->longitude(), latitude, precision );
637
UPSCoordinates* upsCoordinates = ups->convertFromGeodetic( geodeticCoordinates );
638
mgrsorUSNGCoordinates = fromUPS( upsCoordinates, precision );
639
delete upsCoordinates;
643
delete geodeticCoordinates;
644
geodeticCoordinates = 0;
646
return mgrsorUSNGCoordinates;
650
MSP::CCS::UTMCoordinates* USNG::convertToUTM( MSP::CCS::MGRSorUSNGCoordinates* mgrsorUSNGCoordinates )
653
* The function convertToUTM converts an USNG coordinate string
654
* to UTM projection (zone, hemisphere, easting and northing) coordinates
655
* according to the current ellipsoid parameters. If any errors occur,
656
* an exception is thrown with a description of the error.
658
* USNGString : USNG coordinate string (input)
659
* zone : UTM zone (output)
660
* hemisphere : North or South hemisphere (output)
661
* easting : Easting (X) in meters (output)
662
* northing : Northing (Y) in meters (output)
666
long letters[USNG_LETTERS];
667
double usng_easting, usng_northing;
669
UTMCoordinates* utmCoordinates = 0;
670
GeodeticCoordinates* geodeticCoordinates = 0;
671
char errorStatus[50] = "";
673
breakUSNGString( mgrsorUSNGCoordinates->MGRSString(), &zone, letters, &usng_easting, &usng_northing, &precision );
676
utmCoordinates = toUTM( zone, letters, usng_easting, usng_northing, precision );
677
// Convert to geodetic to make sure the coordinates are in the valid utm range
678
geodeticCoordinates = utm->convertToGeodetic( utmCoordinates );
682
UPSCoordinates* upsCoordinates = toUPS( letters, usng_easting, usng_northing );
683
geodeticCoordinates = ups->convertToGeodetic( upsCoordinates );
684
delete upsCoordinates;
686
utmCoordinates = utm->convertFromGeodetic( geodeticCoordinates );
689
delete geodeticCoordinates;
690
geodeticCoordinates = 0;
692
return utmCoordinates;
696
MSP::CCS::MGRSorUSNGCoordinates* USNG::convertFromUPS( MSP::CCS::UPSCoordinates* upsCoordinates, long precision )
699
* The function convertFromUPS converts UPS (hemisphere, easting,
700
* and northing) coordinates to an USNG coordinate string according to
701
* the current ellipsoid parameters. If any errors occur, an
702
* exception is thrown with a description of the error.
704
* hemisphere : Hemisphere either 'N' or 'S' (input)
705
* easting : Easting/X in meters (input)
706
* northing : Northing/Y in meters (input)
707
* precision : Precision level of USNG string (input)
708
* USNGString : USNG coordinate string (output)
712
char errorStatus[50] = "";
714
char hemisphere = upsCoordinates->hemisphere();
715
double easting = upsCoordinates->easting();
716
double northing = upsCoordinates->northing();
718
if ((hemisphere != 'N') && (hemisphere != 'S'))
719
strcat( errorStatus, ErrorMessages::hemisphere );
720
if ((easting < MIN_EAST_NORTH) || (easting > MAX_EAST_NORTH))
721
strcat( errorStatus, ErrorMessages::easting );
722
if ((northing < MIN_EAST_NORTH) || (northing > MAX_EAST_NORTH))
723
strcat( errorStatus, ErrorMessages::northing );
724
if ((precision < 0) || (precision > MAX_PRECISION))
725
strcat( errorStatus, ErrorMessages::precision );
727
if( strlen( errorStatus ) > 0)
728
throw CoordinateConversionException( errorStatus );
730
GeodeticCoordinates* geodeticCoordinates = ups->convertToGeodetic( upsCoordinates );
732
// If the latitude is within the valid mgrs polar range [-90, -80) or [84, 90],
733
// convert to mgrs using the ups path, otherwise convert to mgrs using the utm path
734
MGRSorUSNGCoordinates* mgrsorUSNGCoordinates = 0;
735
double latitude = geodeticCoordinates->latitude();
737
if ((latitude < (MIN_USNG_NON_POLAR_LAT + EPSILON)) || (latitude >= (MAX_USNG_NON_POLAR_LAT - EPSILON)))
738
mgrsorUSNGCoordinates = fromUPS( upsCoordinates, precision );
741
UTMCoordinates* utmCoordinates = utm->convertFromGeodetic( geodeticCoordinates );
743
double longitude = geodeticCoordinates->longitude();
744
mgrsorUSNGCoordinates = fromUTM( utmCoordinates, longitude, latitude, precision );
745
delete utmCoordinates;
749
delete geodeticCoordinates;
750
geodeticCoordinates = 0;
752
return mgrsorUSNGCoordinates;
756
MSP::CCS::UPSCoordinates* USNG::convertToUPS( MSP::CCS::MGRSorUSNGCoordinates* mgrsorUSNGCoordinates )
759
* The function convertToUPS converts an USNG coordinate string
760
* to UPS (hemisphere, easting, and northing) coordinates, according
761
* to the current ellipsoid parameters. If any errors occur, an
762
* exception is thrown with a description of the error.
764
* USNGString : USNG coordinate string (input)
765
* hemisphere : Hemisphere either 'N' or 'S' (output)
766
* easting : Easting/X in meters (output)
767
* northing : Northing/Y in meters (output)
771
long letters[USNG_LETTERS];
774
double usng_northing;
776
UPSCoordinates* upsCoordinates = 0;
777
GeodeticCoordinates* geodeticCoordinates = 0;
779
breakUSNGString( mgrsorUSNGCoordinates->MGRSString(), &zone, letters, &usng_easting, &usng_northing, &precision );
782
upsCoordinates = toUPS( letters, usng_easting, usng_northing );
783
geodeticCoordinates = ups->convertToGeodetic( upsCoordinates );
787
UTMCoordinates* utmCoordinates = toUTM( zone, letters, usng_easting, usng_northing, precision );
788
geodeticCoordinates = utm->convertToGeodetic( utmCoordinates );
789
delete utmCoordinates;
791
upsCoordinates = ups->convertFromGeodetic( geodeticCoordinates );
794
delete geodeticCoordinates;
795
geodeticCoordinates = 0;
797
return upsCoordinates;
801
MSP::CCS::MGRSorUSNGCoordinates* USNG::fromUTM( MSP::CCS::UTMCoordinates* utmCoordinates, double longitude, double latitude, long precision )
804
* The function fromUTM calculates an USNG coordinate string
805
* based on the zone, latitude, easting and northing.
807
* zone : Zone number (input)
808
* latitude : Latitude in radians (input)
809
* easting : Easting (input)
810
* northing : Northing (input)
811
* precision : Precision (input)
812
* USNGString : USNG coordinate string (output)
815
double pattern_offset; /* Pattern offset for 3rd letter */
816
double grid_northing; /* Northing used to derive 3rd letter of USNG */
817
long ltr2_low_value; /* 2nd letter range - low number */
818
long ltr2_high_value; /* 2nd letter range - high number */
819
int letters[USNG_LETTERS]; /* Number location of 3 letters in alphabet */
824
long zone = utmCoordinates->zone();
825
char hemisphere = utmCoordinates->hemisphere();
826
double easting = utmCoordinates->easting();
827
double northing = utmCoordinates->northing();
829
getLatitudeLetter( latitude, &letters[0] );
831
easting = roundUSNG(easting);
833
// Check if the point is within it's natural zone
834
// If it is not, put it there
836
natural_zone = (long)(31 + ((longitude) / _6));
838
natural_zone = (long)(((longitude) / _6) - 29);
840
if (natural_zone > 60)
842
if (zone != natural_zone)
843
{ // reconvert to override zone
844
UTM utmOverride( semiMajorAxis, flattening, natural_zone );
845
GeodeticCoordinates geodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
846
UTMCoordinates* utmCoordinatesOverride = utmOverride.convertFromGeodetic( &geodeticCoordinates );
848
zone = utmCoordinatesOverride->zone();
849
hemisphere = utmCoordinatesOverride->hemisphere();
850
easting = utmCoordinatesOverride->easting();
851
northing = utmCoordinatesOverride->northing();
853
delete utmCoordinatesOverride;
854
utmCoordinatesOverride = 0;
857
easting = roundUSNG(easting);
859
/* UTM special cases */
860
if (letters[0] == LETTER_V) // V latitude band
862
if ((zone == 31) && (easting >= _500000))
863
override = 32; // extension of zone 32V
865
else if (letters[0] == LETTER_X)
867
if ((zone == 32) && (easting < _500000)) // extension of zone 31X
869
else if (((zone == 32) && (easting >= _500000)) || // western extension of zone 33X
870
((zone == 34) && (easting < _500000))) // eastern extension of zone 33X
872
else if (((zone == 34) && (easting >= _500000)) || // western extension of zone 35X
873
((zone == 36) && (easting < _500000))) // eastern extension of zone 35X
875
else if ((zone == 36) && (easting >= _500000)) // western extension of zone 37X
880
{ // reconvert to override zone
881
UTM utmOverride( semiMajorAxis, flattening, override );
882
GeodeticCoordinates geodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
883
UTMCoordinates* utmCoordinatesOverride = utmOverride.convertFromGeodetic( &geodeticCoordinates );
885
zone = utmCoordinatesOverride->zone();
886
hemisphere = utmCoordinatesOverride->hemisphere();
887
easting = utmCoordinatesOverride->easting();
888
northing = utmCoordinatesOverride->northing();
890
delete utmCoordinatesOverride;
891
utmCoordinatesOverride = 0;
894
easting = roundUSNG(easting);
895
northing = roundUSNG(northing);
896
/* Truncate easting and northing values */
897
double divisor = pow (10.0, (5.0 - precision));
898
easting = (long)(easting/divisor) * divisor;
899
northing = (long)(northing/divisor) * divisor;
901
if( latitude <= 0.0 && northing == 1.0e7)
907
getGridValues( zone, <r2_low_value, <r2_high_value, &pattern_offset );
909
grid_northing = northing;
911
while (grid_northing >= TWOMIL)
913
grid_northing = grid_northing - TWOMIL;
915
grid_northing = grid_northing + pattern_offset;
916
if(grid_northing >= TWOMIL)
917
grid_northing = grid_northing - TWOMIL;
919
letters[2] = (long)(grid_northing / ONEHT);
920
if (letters[2] > LETTER_H)
921
letters[2] = letters[2] + 1;
923
if (letters[2] > LETTER_N)
924
letters[2] = letters[2] + 1;
926
letters[1] = ltr2_low_value + ((long)(easting / ONEHT) -1);
927
if ((ltr2_low_value == LETTER_J) && (letters[1] > LETTER_N))
928
letters[1] = letters[1] + 1;
930
makeUSNGString( USNGString, zone, letters, easting, northing, precision );
932
return new MGRSorUSNGCoordinates( CoordinateType::usNationalGrid, USNGString );
936
MSP::CCS::UTMCoordinates* USNG::toUTM( long zone, long letters[USNG_LETTERS], double easting, double northing, long precision )
939
* The function toUTM converts an USNG coordinate string
940
* to UTM projection (zone, hemisphere, easting and northing) coordinates
941
* according to the current ellipsoid parameters. If any errors occur,
942
* an exception is thrown with a description of the error.
944
* USNGString : USNG coordinate string (input)
945
* zone : UTM zone (output)
946
* hemisphere : North or South hemisphere (output)
947
* easting : Easting (X) in meters (output)
948
* northing : Northing (Y) in meters (output)
953
double northing_offset;
955
long ltr2_high_value;
956
double pattern_offset;
957
double upper_lat_limit; /* North latitude limits based on 1st letter */
958
double lower_lat_limit; /* South latitude limits based on 1st letter */
959
double grid_easting; /* Easting for 100,000 meter grid square */
960
double grid_northing; /* Northing for 100,000 meter grid square */
961
double temp_grid_northing = 0.0;
962
double fabs_grid_northing = 0.0;
963
double latitude = 0.0;
964
double longitude = 0.0;
965
double divisor = 1.0;
966
UTMCoordinates* utmCoordinates = 0;
967
char errorStatus[50] = "";
969
if ((letters[0] == LETTER_X) && ((zone == 32) || (zone == 34) || (zone == 36)))
970
throw CoordinateConversionException( ErrorMessages::usngString );
971
else if ((letters[0] == LETTER_V) && (zone == 31) && (letters[1] > LETTER_D))
972
throw CoordinateConversionException( ErrorMessages::usngString );
975
if (letters[0] < LETTER_N)
980
getGridValues(zone, <r2_low_value, <r2_high_value, &pattern_offset);
982
/* Check that the second letter of the USNG string is within
983
* the range of valid second letter values
984
* Also check that the third letter is valid */
985
if ((letters[1] < ltr2_low_value) || (letters[1] > ltr2_high_value) || (letters[2] > LETTER_V))
986
throw CoordinateConversionException( ErrorMessages::usngString );
988
grid_easting = (double)((letters[1]) - ltr2_low_value + 1) * ONEHT;
989
if ((ltr2_low_value == LETTER_J) && (letters[1] > LETTER_O))
990
grid_easting = grid_easting - ONEHT;
992
double row_letter_northing = (double)(letters[2]) * ONEHT;
993
if (letters[2] > LETTER_O)
994
row_letter_northing = row_letter_northing - ONEHT;
996
if (letters[2] > LETTER_I)
997
row_letter_northing = row_letter_northing - ONEHT;
999
if (row_letter_northing >= TWOMIL)
1000
row_letter_northing = row_letter_northing - TWOMIL;
1002
getLatitudeBandMinNorthing(letters[0], &min_northing, &northing_offset);
1004
grid_northing = row_letter_northing - pattern_offset;
1005
if(grid_northing < 0)
1006
grid_northing += TWOMIL;
1008
grid_northing += northing_offset;
1010
if(grid_northing < min_northing)
1011
grid_northing += TWOMIL;
1013
easting = grid_easting + easting;
1014
northing = grid_northing + northing;
1016
utmCoordinates = new UTMCoordinates( CoordinateType::universalTransverseMercator, zone, hemisphere, easting, northing );
1018
/* check that point is within Zone Letter bounds */
1019
GeodeticCoordinates* geodeticCoordinates = utm->convertToGeodetic( utmCoordinates );
1021
divisor = pow (10.0, (double)precision);
1022
getLatitudeRange(letters[0], &upper_lat_limit, &lower_lat_limit);
1024
double latitude = geodeticCoordinates->latitude();
1026
delete geodeticCoordinates;
1027
geodeticCoordinates = 0;
1029
if (!(((lower_lat_limit - PI_OVER_180/divisor) <= latitude) && (latitude <= (upper_lat_limit + PI_OVER_180/divisor))))
1030
utmCoordinates->setWarningMessage( MSP::CCS::WarningMessages::latitude );
1033
return utmCoordinates;
1037
MSP::CCS::MGRSorUSNGCoordinates* USNG::fromUPS( MSP::CCS::UPSCoordinates* upsCoordinates, long precision )
1040
* The function fromUPS converts UPS (hemisphere, easting,
1041
* and northing) coordinates to an USNG coordinate string according to
1042
* the current ellipsoid parameters.
1044
* hemisphere : Hemisphere either 'N' or 'S' (input)
1045
* easting : Easting/X in meters (input)
1046
* northing : Northing/Y in meters (input)
1047
* precision : Precision level of USNG string (input)
1048
* USNGString : USNG coordinate string (output)
1051
double false_easting; /* False easting for 2nd letter */
1052
double false_northing; /* False northing for 3rd letter */
1053
double grid_easting; /* Easting used to derive 2nd letter of USNG */
1054
double grid_northing; /* Northing used to derive 3rd letter of USNG */
1055
long ltr2_low_value; /* 2nd letter range - low number */
1056
int letters[USNG_LETTERS]; /* Number location of 3 letters in alphabet */
1059
char USNGString[21];
1061
char hemisphere = upsCoordinates->hemisphere();
1062
double easting = upsCoordinates->easting();
1063
double northing = upsCoordinates->northing();
1065
easting = roundUSNG(easting);
1066
northing = roundUSNG(northing);
1067
divisor = pow (10.0, (5.0 - precision));
1068
easting = (long)(easting/divisor) * divisor;
1069
northing = (long)(northing/divisor) * divisor;
1071
if (hemisphere == 'N')
1073
if (easting >= TWOMIL)
1074
letters[0] = LETTER_Z;
1076
letters[0] = LETTER_Y;
1078
index = letters[0] - 22;
1079
ltr2_low_value = UPS_Constant_Table[index].ltr2_low_value;
1080
false_easting = UPS_Constant_Table[index].false_easting;
1081
false_northing = UPS_Constant_Table[index].false_northing;
1085
if (easting >= TWOMIL)
1086
letters[0] = LETTER_B;
1088
letters[0] = LETTER_A;
1090
ltr2_low_value = UPS_Constant_Table[letters[0]].ltr2_low_value;
1091
false_easting = UPS_Constant_Table[letters[0]].false_easting;
1092
false_northing = UPS_Constant_Table[letters[0]].false_northing;
1095
grid_northing = northing;
1096
grid_northing = grid_northing - false_northing;
1097
letters[2] = (long)(grid_northing / ONEHT);
1099
if (letters[2] > LETTER_H)
1100
letters[2] = letters[2] + 1;
1102
if (letters[2] > LETTER_N)
1103
letters[2] = letters[2] + 1;
1105
grid_easting = easting;
1106
grid_easting = grid_easting - false_easting;
1107
letters[1] = ltr2_low_value + ((long)(grid_easting / ONEHT));
1109
if (easting < TWOMIL)
1111
if (letters[1] > LETTER_L)
1112
letters[1] = letters[1] + 3;
1114
if (letters[1] > LETTER_U)
1115
letters[1] = letters[1] + 2;
1119
if (letters[1] > LETTER_C)
1120
letters[1] = letters[1] + 2;
1122
if (letters[1] > LETTER_H)
1123
letters[1] = letters[1] + 1;
1125
if (letters[1] > LETTER_L)
1126
letters[1] = letters[1] + 3;
1129
makeUSNGString( USNGString, 0, letters, easting, northing, precision );
1131
return new MGRSorUSNGCoordinates( CoordinateType::usNationalGrid, USNGString );
1135
MSP::CCS::UPSCoordinates* USNG::toUPS( long letters[USNG_LETTERS], double easting, double northing )
1138
* The function toUPS converts an USNG coordinate string
1139
* to UPS (hemisphere, easting, and northing) coordinates, according
1140
* to the current ellipsoid parameters. If any errors occur, an
1141
* exception is thrown with a description of the error.
1143
* USNGString : USNG coordinate string (input)
1144
* hemisphere : Hemisphere either 'N' or 'S' (output)
1145
* easting : Easting/X in meters (output)
1146
* northing : Northing/Y in meters (output)
1149
long ltr2_high_value; /* 2nd letter range - high number */
1150
long ltr3_high_value; /* 3rd letter range - high number (UPS) */
1151
long ltr2_low_value; /* 2nd letter range - low number */
1152
double false_easting; /* False easting for 2nd letter */
1153
double false_northing; /* False northing for 3rd letter */
1154
double grid_easting; /* easting for 100,000 meter grid square */
1155
double grid_northing; /* northing for 100,000 meter grid square */
1159
if ((letters[0] == LETTER_Y) || (letters[0] == LETTER_Z))
1163
index = letters[0] - 22;
1164
ltr2_low_value = UPS_Constant_Table[index].ltr2_low_value;
1165
ltr2_high_value = UPS_Constant_Table[index].ltr2_high_value;
1166
ltr3_high_value = UPS_Constant_Table[index].ltr3_high_value;
1167
false_easting = UPS_Constant_Table[index].false_easting;
1168
false_northing = UPS_Constant_Table[index].false_northing;
1170
else if ((letters[0] == LETTER_A) || (letters[0] == LETTER_B))
1174
ltr2_low_value = UPS_Constant_Table[letters[0]].ltr2_low_value;
1175
ltr2_high_value = UPS_Constant_Table[letters[0]].ltr2_high_value;
1176
ltr3_high_value = UPS_Constant_Table[letters[0]].ltr3_high_value;
1177
false_easting = UPS_Constant_Table[letters[0]].false_easting;
1178
false_northing = UPS_Constant_Table[letters[0]].false_northing;
1181
throw CoordinateConversionException( ErrorMessages::usngString );
1183
/* Check that the second letter of the USNG string is within
1184
* the range of valid second letter values
1185
* Also check that the third letter is valid */
1186
if ((letters[1] < ltr2_low_value) || (letters[1] > ltr2_high_value) ||
1187
((letters[1] == LETTER_D) || (letters[1] == LETTER_E) ||
1188
(letters[1] == LETTER_M) || (letters[1] == LETTER_N) ||
1189
(letters[1] == LETTER_V) || (letters[1] == LETTER_W)) ||
1190
(letters[2] > ltr3_high_value))
1191
throw CoordinateConversionException( ErrorMessages::usngString );
1193
grid_northing = (double)letters[2] * ONEHT + false_northing;
1194
if (letters[2] > LETTER_I)
1195
grid_northing = grid_northing - ONEHT;
1197
if (letters[2] > LETTER_O)
1198
grid_northing = grid_northing - ONEHT;
1200
grid_easting = (double)((letters[1]) - ltr2_low_value) * ONEHT + false_easting;
1201
if (ltr2_low_value != LETTER_A)
1203
if (letters[1] > LETTER_L)
1204
grid_easting = grid_easting - 300000.0;
1206
if (letters[1] > LETTER_U)
1207
grid_easting = grid_easting - 200000.0;
1211
if (letters[1] > LETTER_C)
1212
grid_easting = grid_easting - 200000.0;
1214
if (letters[1] > LETTER_I)
1215
grid_easting = grid_easting - ONEHT;
1217
if (letters[1] > LETTER_L)
1218
grid_easting = grid_easting - 300000.0;
1221
easting = grid_easting + easting;
1222
northing = grid_northing + northing;
1224
return new UPSCoordinates( CoordinateType::universalPolarStereographic, hemisphere, easting, northing );
1228
void USNG::getGridValues( long zone, long* ltr2_low_value, long* ltr2_high_value, double *pattern_offset )
1231
* The function getGridValues sets the letter range used for
1232
* the 2nd letter in the USNG coordinate string, based on the set
1233
* number of the utm zone. It also sets the pattern offset using a
1234
* value of A for the second letter of the grid square, based on
1235
* the grid pattern and set number of the utm zone.
1237
* zone : Zone number (input)
1238
* ltr2_low_value : 2nd letter low number (output)
1239
* ltr2_high_value : 2nd letter high number (output)
1240
* pattern_offset : Pattern offset (output)
1243
long set_number; /* Set number (1-6) based on UTM zone number */
1245
set_number = zone % 6;
1250
if ((set_number == 1) || (set_number == 4))
1252
*ltr2_low_value = LETTER_A;
1253
*ltr2_high_value = LETTER_H;
1255
else if ((set_number == 2) || (set_number == 5))
1257
*ltr2_low_value = LETTER_J;
1258
*ltr2_high_value = LETTER_R;
1260
else if ((set_number == 3) || (set_number == 6))
1262
*ltr2_low_value = LETTER_S;
1263
*ltr2_high_value = LETTER_Z;
1266
/* False northing at A for second letter of grid square */
1267
if ((set_number % 2) == 0)
1268
*pattern_offset = 500000.0;
1270
*pattern_offset = 0.0;
1274
void USNG::getLatitudeBandMinNorthing( long letter, double* min_northing, double* northing_offset )
1277
* The function getLatitudeBandMinNorthing receives a latitude band letter
1278
* and uses the Latitude_Band_Table to determine the minimum northing and northing offset
1279
* for that latitude band letter.
1281
* letter : Latitude band letter (input)
1282
* min_northing : Minimum northing for that letter (output)
1283
* northing_offset : Latitude band northing offset (output)
1286
if ((letter >= LETTER_C) && (letter <= LETTER_H))
1288
*min_northing = Latitude_Band_Table[letter-2].min_northing;
1289
*northing_offset = Latitude_Band_Table[letter-2].northing_offset;
1291
else if ((letter >= LETTER_J) && (letter <= LETTER_N))
1293
*min_northing = Latitude_Band_Table[letter-3].min_northing;
1294
*northing_offset = Latitude_Band_Table[letter-3].northing_offset;
1296
else if ((letter >= LETTER_P) && (letter <= LETTER_X))
1298
*min_northing = Latitude_Band_Table[letter-4].min_northing;
1299
*northing_offset = Latitude_Band_Table[letter-4].northing_offset;
1302
throw CoordinateConversionException( ErrorMessages::usngString );
1306
void USNG::getLatitudeRange( long letter, double* north, double* south )
1309
* The function getLatitudeRange receives a latitude band letter
1310
* and uses the Latitude_Band_Table to determine the latitude band
1311
* boundaries for that latitude band letter.
1313
* letter : Latitude band letter (input)
1314
* north : Northern latitude boundary for that letter (output)
1315
* north : Southern latitude boundary for that letter (output)
1318
if ((letter >= LETTER_C) && (letter <= LETTER_H))
1320
*north = Latitude_Band_Table[letter-2].north * PI_OVER_180;
1321
*south = Latitude_Band_Table[letter-2].south * PI_OVER_180;
1323
else if ((letter >= LETTER_J) && (letter <= LETTER_N))
1325
*north = Latitude_Band_Table[letter-3].north * PI_OVER_180;
1326
*south = Latitude_Band_Table[letter-3].south * PI_OVER_180;
1328
else if ((letter >= LETTER_P) && (letter <= LETTER_X))
1330
*north = Latitude_Band_Table[letter-4].north * PI_OVER_180;
1331
*south = Latitude_Band_Table[letter-4].south * PI_OVER_180;
1334
throw CoordinateConversionException( ErrorMessages::usngString );
1338
void USNG::getLatitudeLetter( double latitude, int* letter )
1341
* The function getLatitudeLetter receives a latitude value
1342
* and uses the Latitude_Band_Table to determine the latitude band
1343
* letter for that latitude.
1345
* latitude : Latitude (input)
1346
* letter : Latitude band letter (output)
1351
if (latitude >= _72 && latitude < _84_5)
1353
else if (latitude > -_80_5 && latitude < _72)
1355
band = (long)(((latitude + _80) / _8) + 1.0e-12);
1358
*letter = Latitude_Band_Table[band].letter;
1361
throw CoordinateConversionException( ErrorMessages::latitude );
1368
// CLASSIFICATION: UNCLASSIFIED
added
removed
.pc/002-gxx-4-changes.patch/CCS/src/dtcc/CoordinateSystems/usng/USNG.cpp
