1 /***************************************************************************
2  * $Id: mgrs.c 10645 2007-01-18 02:22:39Z warmerdam $
3  *
4  * Project:  MGRS Converter
5  * Purpose:  Geotrans code for MGRS translation (slightly adapted)
6  * Author:   Unknown (NIMA)
7  *
8  ***************************************************************************
9  ***************************************************************************
10  * RSC IDENTIFIER:  MGRS
11  *
12  * ABSTRACT
13  *
14  *    This component converts between geodetic coordinates (latitude and
15  *    longitude) and Military Grid Reference System (MGRS) coordinates.
16  *
17  * ERROR HANDLING
18  *
19  *    This component checks parameters for valid values.  If an invalid value
20  *    is found, the error code is combined with the current error code using
21  *    the bitwise or.  This combining allows multiple error codes to be
22  *    returned. The possible error codes are:
23  *
24  *          MGRS_NO_ERROR          : No errors occurred in function
25  *          MGRS_LAT_ERROR         : Latitude outside of valid range
26  *                                    (-90 to 90 degrees)
27  *          MGRS_LON_ERROR         : Longitude outside of valid range
28  *                                    (-180 to 360 degrees)
29  *          MGRS_STR_ERROR         : An MGRS string error: string too long,
30  *                                    too short, or badly formed
31  *          MGRS_PRECISION_ERROR   : The precision must be between 0 and 5
32  *                                    inclusive.
33  *          MGRS_A_ERROR           : Semi-major axis less than or equal to zero
34  *          MGRS_INV_F_ERROR       : Inverse flattening outside of valid range
35  *									                  (250 to 350)
36  *          MGRS_EASTING_ERROR     : Easting outside of valid range
37  *                                    (100,000 to 900,000 meters for UTM)
38  *                                    (0 to 4,000,000 meters for UPS)
39  *          MGRS_NORTHING_ERROR    : Northing outside of valid range
40  *                                    (0 to 10,000,000 meters for UTM)
41  *                                    (0 to 4,000,000 meters for UPS)
42  *          MGRS_ZONE_ERROR        : Zone outside of valid range (1 to 60)
43  *          MGRS_HEMISPHERE_ERROR  : Invalid hemisphere ('N' or 'S')
44  *
45  * REUSE NOTES
46  *
47  *    MGRS is intended for reuse by any application that does conversions
48  *    between geodetic coordinates and MGRS coordinates.
49  *
50  * REFERENCES
51  *
52  *    Further information on MGRS can be found in the Reuse Manual.
53  *
54  *    MGRS originated from : U.S. Army Topographic Engineering Center
55  *                           Geospatial Information Division
56  *                           7701 Telegraph Road
57  *                           Alexandria, VA  22310-3864
58  *
59  * LICENSES
60  *
61  *    None apply to this component.
62  *
63  * RESTRICTIONS
64  *
65  *
66  * ENVIRONMENT
67  *
68  *    MGRS was tested and certified in the following environments:
69  *
70  *    1. Solaris 2.5 with GCC version 2.8.1
71  *    2. Windows 95 with MS Visual C++ version 6
72  *
73  * MODIFICATIONS
74  *
75  *    Date              Description
76  *    ----              -----------
77  *    16-11-94          Original Code
78  *    15-09-99          Reengineered upper layers
79  *    02-05-03          Corrected latitude band bug in GRID_UTM
80  *    08-20-03          Reengineered lower layers
81  */
82 
83 
84 /***************************************************************************/
85 /*
86  *                               INCLUDES
87  */
88 #include <ctype.h>
89 #include <math.h>
90 #include <stdio.h>
91 #include <string.h>
92 #include "mgrs.h"
93 
94 /*
95  *      ctype.h     - Standard C character handling library
96  *      math.h      - Standard C math library
97  *      stdio.h     - Standard C input/output library
98  *      string.h    - Standard C string handling library
99  *      ups.h       - Universal Polar Stereographic (UPS) projection
100  *      utm.h       - Universal Transverse Mercator (UTM) projection
101  *      mgrs.h      - function prototype error checking
102  */
103 
104 
105 /***************************************************************************/
106 /*
107  *                              GLOBAL DECLARATIONS
108  */
109 #define DEG_TO_RAD       0.017453292519943295 /* PI/180                      */
110 #define RAD_TO_DEG       57.29577951308232087 /* 180/PI                      */
111 #define LETTER_A               0   /* ARRAY INDEX FOR LETTER A               */
112 #define LETTER_B               1   /* ARRAY INDEX FOR LETTER B               */
113 #define LETTER_C               2   /* ARRAY INDEX FOR LETTER C               */
114 #define LETTER_D               3   /* ARRAY INDEX FOR LETTER D               */
115 #define LETTER_E               4   /* ARRAY INDEX FOR LETTER E               */
116 #define LETTER_F               5   /* ARRAY INDEX FOR LETTER E               */
117 #define LETTER_G               6   /* ARRAY INDEX FOR LETTER H               */
118 #define LETTER_H               7   /* ARRAY INDEX FOR LETTER H               */
119 #define LETTER_I               8   /* ARRAY INDEX FOR LETTER I               */
120 #define LETTER_J               9   /* ARRAY INDEX FOR LETTER J               */
121 #define LETTER_K              10   /* ARRAY INDEX FOR LETTER J               */
122 #define LETTER_L              11   /* ARRAY INDEX FOR LETTER L               */
123 #define LETTER_M              12   /* ARRAY INDEX FOR LETTER M               */
124 #define LETTER_N              13   /* ARRAY INDEX FOR LETTER N               */
125 #define LETTER_O              14   /* ARRAY INDEX FOR LETTER O               */
126 #define LETTER_P              15   /* ARRAY INDEX FOR LETTER P               */
127 #define LETTER_Q              16   /* ARRAY INDEX FOR LETTER Q               */
128 #define LETTER_R              17   /* ARRAY INDEX FOR LETTER R               */
129 #define LETTER_S              18   /* ARRAY INDEX FOR LETTER S               */
130 #define LETTER_T              19   /* ARRAY INDEX FOR LETTER S               */
131 #define LETTER_U              20   /* ARRAY INDEX FOR LETTER U               */
132 #define LETTER_V              21   /* ARRAY INDEX FOR LETTER V               */
133 #define LETTER_W              22   /* ARRAY INDEX FOR LETTER W               */
134 #define LETTER_X              23   /* ARRAY INDEX FOR LETTER X               */
135 #define LETTER_Y              24   /* ARRAY INDEX FOR LETTER Y               */
136 #define LETTER_Z              25   /* ARRAY INDEX FOR LETTER Z               */
137 #define MGRS_LETTERS            3  /* NUMBER OF LETTERS IN MGRS              */
138 #define ONEHT          100000.e0    /* ONE HUNDRED THOUSAND                  */
139 #define TWOMIL        2000000.e0    /* TWO MILLION                           */
140 #define TRUE                      1  /* CONSTANT VALUE FOR TRUE VALUE  */
141 #define FALSE                     0  /* CONSTANT VALUE FOR FALSE VALUE */
142 #define PI    3.14159265358979323e0  /* PI                             */
143 #define PI_OVER_2  (PI / 2.0e0)
144 
145 #define MIN_EASTING  100000
146 #define MAX_EASTING  900000
147 #define MIN_NORTHING 0
148 #define MAX_NORTHING 10000000
149 #define MAX_PRECISION           5   /* Maximum precision of easting & northing */
150 #define MIN_UTM_LAT      ( (-80 * PI) / 180.0 ) /* -80 degrees in radians    */
151 #define MAX_UTM_LAT      ( (84 * PI) / 180.0 )  /* 84 degrees in radians     */
152 
153 #define MIN_EAST_NORTH 0
154 #define MAX_EAST_NORTH 4000000
155 
156 
157 /* Ellipsoid parameters, default to WGS 84 */
158 double MGRS_a = 6378137.0;    /* Semi-major axis of ellipsoid in meters */
159 double MGRS_f = 1 / 298.257223563; /* Flattening of ellipsoid           */
160 double MGRS_recpf = 298.257223563;
161 char   MGRS_Ellipsoid_Code[3] = {'W','E',0};
162 
163 
164 /*
165  *    CLARKE_1866 : Ellipsoid code for CLARKE_1866
166  *    CLARKE_1880 : Ellipsoid code for CLARKE_1880
167  *    BESSEL_1841 : Ellipsoid code for BESSEL_1841
168  *    BESSEL_1841_NAMIBIA : Ellipsoid code for BESSEL 1841 (NAMIBIA)
169  */
170 const char* CLARKE_1866 = "CC";
171 const char* CLARKE_1880 = "CD";
172 const char* BESSEL_1841 = "BR";
173 const char* BESSEL_1841_NAMIBIA = "BN";
174 
175 
176 typedef struct Latitude_Band_Value
177 {
178   long letter;            /* letter representing latitude band  */
179   double min_northing;    /* minimum northing for latitude band */
180   double north;           /* upper latitude for latitude band   */
181   double south;           /* lower latitude for latitude band   */
182 } Latitude_Band;
183 
184 static const Latitude_Band Latitude_Band_Table[20] =
185   {{LETTER_C, 1100000.0, -72.0, -80.5},
186   {LETTER_D, 2000000.0, -64.0, -72.0},
187   {LETTER_E, 2800000.0, -56.0, -64.0},
188   {LETTER_F, 3700000.0, -48.0, -56.0},
189   {LETTER_G, 4600000.0, -40.0, -48.0},
190   {LETTER_H, 5500000.0, -32.0, -40.0},
191   {LETTER_J, 6400000.0, -24.0, -32.0},
192   {LETTER_K, 7300000.0, -16.0, -24.0},
193   {LETTER_L, 8200000.0, -8.0, -16.0},
194   {LETTER_M, 9100000.0, 0.0, -8.0},
195   {LETTER_N, 0.0, 8.0, 0.0},
196   {LETTER_P, 800000.0, 16.0, 8.0},
197   {LETTER_Q, 1700000.0, 24.0, 16.0},
198   {LETTER_R, 2600000.0, 32.0, 24.0},
199   {LETTER_S, 3500000.0, 40.0, 32.0},
200   {LETTER_T, 4400000.0, 48.0, 40.0},
201   {LETTER_U, 5300000.0, 56.0, 48.0},
202   {LETTER_V, 6200000.0, 64.0, 56.0},
203   {LETTER_W, 7000000.0, 72.0, 64.0},
204   {LETTER_X, 7900000.0, 84.5, 72.0}};
205 
206 
207 typedef struct UPS_Constant_Value
208 {
209   long letter;            /* letter representing latitude band      */
210   long ltr2_low_value;    /* 2nd letter range - high number         */
211   long ltr2_high_value;   /* 2nd letter range - low number          */
212   long ltr3_high_value;   /* 3rd letter range - high number (UPS)   */
213   double false_easting;   /* False easting based on 2nd letter      */
214   double false_northing;  /* False northing based on 3rd letter     */
215 } UPS_Constant;
216 
217 static const UPS_Constant UPS_Constant_Table[4] =
218   {{LETTER_A, LETTER_J, LETTER_Z, LETTER_Z, 800000.0, 800000.0},
219   {LETTER_B, LETTER_A, LETTER_R, LETTER_Z, 2000000.0, 800000.0},
220   {LETTER_Y, LETTER_J, LETTER_Z, LETTER_P, 800000.0, 1300000.0},
221   {LETTER_Z, LETTER_A, LETTER_J, LETTER_P, 2000000.0, 1300000.0}};
222 
223 /***************************************************************************/
224 /*
225  *                              FUNCTIONS
226  */
227 
Get_Latitude_Band_Min_Northing(long letter,double * min_northing)228 long Get_Latitude_Band_Min_Northing(long letter, double* min_northing)
229 /*
230  * The function Get_Latitude_Band_Min_Northing receives a latitude band letter
231  * and uses the Latitude_Band_Table to determine the minimum northing for that
232  * latitude band letter.
233  *
234  *   letter        : Latitude band letter             (input)
235  *   min_northing  : Minimum northing for that letter	(output)
236  */
237 { /* Get_Latitude_Band_Min_Northing */
238   long error_code = MGRS_NO_ERROR;
239 
240   if ((letter >= LETTER_C) && (letter <= LETTER_H))
241     *min_northing = Latitude_Band_Table[letter-2].min_northing;
242   else if ((letter >= LETTER_J) && (letter <= LETTER_N))
243     *min_northing = Latitude_Band_Table[letter-3].min_northing;
244   else if ((letter >= LETTER_P) && (letter <= LETTER_X))
245     *min_northing = Latitude_Band_Table[letter-4].min_northing;
246   else
247     error_code |= MGRS_STRING_ERROR;
248 
249   return error_code;
250 } /* Get_Latitude_Band_Min_Northing */
251 
252 
Get_Latitude_Range(long letter,double * north,double * south)253 long Get_Latitude_Range(long letter, double* north, double* south)
254 /*
255  * The function Get_Latitude_Range receives a latitude band letter
256  * and uses the Latitude_Band_Table to determine the latitude band
257  * boundaries for that latitude band letter.
258  *
259  *   letter   : Latitude band letter                        (input)
260  *   north    : Northern latitude boundary for that letter	(output)
261  *   north    : Southern latitude boundary for that letter	(output)
262  */
263 { /* Get_Latitude_Range */
264   long error_code = MGRS_NO_ERROR;
265 
266   if ((letter >= LETTER_C) && (letter <= LETTER_H))
267   {
268     *north = Latitude_Band_Table[letter-2].north * DEG_TO_RAD;
269     *south = Latitude_Band_Table[letter-2].south * DEG_TO_RAD;
270   }
271   else if ((letter >= LETTER_J) && (letter <= LETTER_N))
272   {
273     *north = Latitude_Band_Table[letter-3].north * DEG_TO_RAD;
274     *south = Latitude_Band_Table[letter-3].south * DEG_TO_RAD;
275   }
276   else if ((letter >= LETTER_P) && (letter <= LETTER_X))
277   {
278     *north = Latitude_Band_Table[letter-4].north * DEG_TO_RAD;
279     *south = Latitude_Band_Table[letter-4].south * DEG_TO_RAD;
280   }
281   else
282     error_code |= MGRS_STRING_ERROR;
283 
284   return error_code;
285 } /* Get_Latitude_Range */
286 
287 
Get_Latitude_Letter(double latitude,int * letter)288 long Get_Latitude_Letter(double latitude, int* letter)
289 /*
290  * The function Get_Latitude_Letter receives a latitude value
291  * and uses the Latitude_Band_Table to determine the latitude band
292  * letter for that latitude.
293  *
294  *   latitude   : Latitude              (input)
295  *   letter     : Latitude band letter  (output)
296  */
297 { /* Get_Latitude_Letter */
298   double temp = 0.0;
299   long error_code = MGRS_NO_ERROR;
300   double lat_deg = latitude * RAD_TO_DEG;
301 
302   if (lat_deg >= 72 && lat_deg < 84.5)
303     *letter = LETTER_X;
304   else if (lat_deg > -80.5 && lat_deg < 72)
305   {
306     temp = ((latitude + (80.0 * DEG_TO_RAD)) / (8.0 * DEG_TO_RAD)) + 1.0e-12;
307     *letter = Latitude_Band_Table[(int)temp].letter;
308   }
309   else
310     error_code |= MGRS_LAT_ERROR;
311 
312   return error_code;
313 } /* Get_Latitude_Letter */
314 
315 
Check_Zone(char * MGRS,long * zone_exists)316 long Check_Zone(char* MGRS, long* zone_exists)
317 /*
318  * The function Check_Zone receives an MGRS coordinate string.
319  * If a zone is given, TRUE is returned. Otherwise, FALSE
320  * is returned.
321  *
322  *   MGRS           : MGRS coordinate string        (input)
323  *   zone_exists    : TRUE if a zone is given,
324  *                    FALSE if a zone is not given  (output)
325  */
326 { /* Check_Zone */
327   int i = 0;
328   int j = 0;
329   int num_digits = 0;
330   long error_code = MGRS_NO_ERROR;
331 
332   /* skip any leading blanks */
333   while (MGRS[i] == ' ')
334     i++;
335   j = i;
336   while (isdigit(MGRS[i]))
337     i++;
338   num_digits = i - j;
339   if (num_digits <= 2)
340     if (num_digits > 0)
341       *zone_exists = TRUE;
342     else
343       *zone_exists = FALSE;
344   else
345     error_code |= MGRS_STRING_ERROR;
346 
347   return error_code;
348 } /* Check_Zone */
349 
350 
Round_MGRS(double value)351 long Round_MGRS (double value)
352 /*
353  * The function Round_MGRS rounds the input value to the
354  * nearest integer, using the standard engineering rule.
355  * The rounded integer value is then returned.
356  *
357  *   value           : Value to be rounded  (input)
358  */
359 { /* Round_MGRS */
360   double ivalue;
361   long ival;
362   double fraction = modf (value, &ivalue);
363   ival = (long)(ivalue);
364   if ((fraction > 0.5) || ((fraction == 0.5) && (ival%2 == 1)))
365     ival++;
366   return (ival);
367 } /* Round_MGRS */
368 
369 
Make_MGRS_String(char * MGRS,long Zone,int Letters[MGRS_LETTERS],double Easting,double Northing,long Precision)370 long Make_MGRS_String (char* MGRS,
371                        long Zone,
372                        int Letters[MGRS_LETTERS],
373                        double Easting,
374                        double Northing,
375                        long Precision)
376 /*
377  * The function Make_MGRS_String constructs an MGRS string
378  * from its component parts.
379  *
380  *   MGRS           : MGRS coordinate string          (output)
381  *   Zone           : UTM Zone                        (input)
382  *   Letters        : MGRS coordinate string letters  (input)
383  *   Easting        : Easting value                   (input)
384  *   Northing       : Northing value                  (input)
385  *   Precision      : Precision level of MGRS string  (input)
386  */
387 { /* Make_MGRS_String */
388   long i;
389   long j;
390   double divisor;
391   long east;
392   long north;
393   char alphabet[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
394   long error_code = MGRS_NO_ERROR;
395 
396   i = 0;
397   if (Zone)
398     i = sprintf (MGRS+i,"%2.2ld",Zone);
399   else
400     strncpy(MGRS, "  ", 2);  // 2 spaces
401 
402   for (j=0;j<3;j++)
403     MGRS[i++] = alphabet[Letters[j]];
404   divisor = pow (10.0, (5 - Precision));
405   Easting = fmod (Easting, 100000.0);
406   if (Easting >= 99999.5)
407     Easting = 99999.0;
408   east = (long)(Easting/divisor);
409   i += sprintf (MGRS+i, "%*.*ld", (int) Precision, (int) Precision, east);
410   Northing = fmod (Northing, 100000.0);
411   if (Northing >= 99999.5)
412     Northing = 99999.0;
413   north = (long)(Northing/divisor);
414   i += sprintf (MGRS+i, "%*.*ld", (int) Precision, (int) Precision, north);
415   return (error_code);
416 } /* Make_MGRS_String */
417 
418 
Break_MGRS_String(char * MGRS,long * Zone,long Letters[MGRS_LETTERS],double * Easting,double * Northing,long * Precision)419 long Break_MGRS_String (char* MGRS,
420                         long* Zone,
421                         long Letters[MGRS_LETTERS],
422                         double* Easting,
423                         double* Northing,
424                         long* Precision)
425 /*
426  * The function Break_MGRS_String breaks down an MGRS
427  * coordinate string into its component parts.
428  *
429  *   MGRS           : MGRS coordinate string          (input)
430  *   Zone           : UTM Zone                        (output)
431  *   Letters        : MGRS coordinate string letters  (output)
432  *   Easting        : Easting value                   (output)
433  *   Northing       : Northing value                  (output)
434  *   Precision      : Precision level of MGRS string  (output)
435  */
436 { /* Break_MGRS_String */
437   long num_digits;
438   long num_letters;
439   long i = 0;
440   long j = 0;
441   long error_code = MGRS_NO_ERROR;
442 
443   while (MGRS[i] == ' ')
444     i++;  /* skip any leading blanks */
445   j = i;
446   while (isdigit(MGRS[i]))
447     i++;
448   num_digits = i - j;
449   if (num_digits <= 2)
450     if (num_digits > 0)
451     {
452       char zone_string[3];
453       /* get zone */
454       strncpy (zone_string, MGRS+j, 2);
455       zone_string[2] = 0;
456       sscanf (zone_string, "%ld", Zone);
457       if ((*Zone < 1) || (*Zone > 60))
458         error_code |= MGRS_STRING_ERROR;
459     }
460     else
461       *Zone = 0;
462   else
463     error_code |= MGRS_STRING_ERROR;
464   j = i;
465 
466   while (isalpha(MGRS[i]))
467     i++;
468   num_letters = i - j;
469   if (num_letters == 3)
470   {
471     /* get letters */
472     Letters[0] = (toupper(MGRS[j]) - (long)'A');
473     if ((Letters[0] == LETTER_I) || (Letters[0] == LETTER_O))
474       error_code |= MGRS_STRING_ERROR;
475     Letters[1] = (toupper(MGRS[j+1]) - (long)'A');
476     if ((Letters[1] == LETTER_I) || (Letters[1] == LETTER_O))
477       error_code |= MGRS_STRING_ERROR;
478     Letters[2] = (toupper(MGRS[j+2]) - (long)'A');
479     if ((Letters[2] == LETTER_I) || (Letters[2] == LETTER_O))
480       error_code |= MGRS_STRING_ERROR;
481   }
482   else
483     error_code |= MGRS_STRING_ERROR;
484   j = i;
485   while (isdigit(MGRS[i]))
486     i++;
487   num_digits = i - j;
488   if ((num_digits <= 10) && (num_digits%2 == 0))
489   {
490     long n;
491     char east_string[6];
492     char north_string[6];
493     long east;
494     long north;
495     double multiplier;
496     /* get easting & northing */
497     n = num_digits/2;
498     *Precision = n;
499     if (n > 0)
500     {
501       strncpy (east_string, MGRS+j, n);
502       east_string[n] = 0;
503       sscanf (east_string, "%ld", &east);
504       strncpy (north_string, MGRS+j+n, n);
505       north_string[n] = 0;
506       sscanf (north_string, "%ld", &north);
507       multiplier = pow (10.0, 5 - n);
508       *Easting = east * multiplier;
509       *Northing = north * multiplier;
510     }
511     else
512     {
513       *Easting = 0.0;
514       *Northing = 0.0;
515     }
516   }
517   else
518     error_code |= MGRS_STRING_ERROR;
519 
520   return (error_code);
521 } /* Break_MGRS_String */
522 
523 
Get_Grid_Values(long zone,long * ltr2_low_value,long * ltr2_high_value,double * false_northing)524 void Get_Grid_Values (long zone,
525                       long* ltr2_low_value,
526                       long* ltr2_high_value,
527                       double *false_northing)
528 /*
529  * The function Get_Grid_Values sets the letter range used for
530  * the 2nd letter in the MGRS coordinate string, based on the set
531  * number of the utm zone. It also sets the false northing using a
532  * value of A for the second letter of the grid square, based on
533  * the grid pattern and set number of the utm zone.
534  *
535  *    zone            : Zone number             (input)
536  *    ltr2_low_value  : 2nd letter low number   (output)
537  *    ltr2_high_value : 2nd letter high number  (output)
538  *    false_northing  : False northing          (output)
539  */
540 { /* BEGIN Get_Grid_Values */
541   long set_number;    /* Set number (1-6) based on UTM zone number */
542   long aa_pattern;    /* Pattern based on ellipsoid code */
543 
544   set_number = zone % 6;
545 
546   if (!set_number)
547     set_number = 6;
548 
549   if (!strcmp(MGRS_Ellipsoid_Code,CLARKE_1866) || !strcmp(MGRS_Ellipsoid_Code, CLARKE_1880) ||
550       !strcmp(MGRS_Ellipsoid_Code,BESSEL_1841) || !strcmp(MGRS_Ellipsoid_Code,BESSEL_1841_NAMIBIA))
551     aa_pattern = FALSE;
552   else
553     aa_pattern = TRUE;
554 
555   if ((set_number == 1) || (set_number == 4))
556   {
557     *ltr2_low_value = LETTER_A;
558     *ltr2_high_value = LETTER_H;
559   }
560   else if ((set_number == 2) || (set_number == 5))
561   {
562     *ltr2_low_value = LETTER_J;
563     *ltr2_high_value = LETTER_R;
564   }
565   else if ((set_number == 3) || (set_number == 6))
566   {
567     *ltr2_low_value = LETTER_S;
568     *ltr2_high_value = LETTER_Z;
569   }
570 
571   /* False northing at A for second letter of grid square */
572   if (aa_pattern)
573   {
574     if ((set_number % 2) ==  0)
575       *false_northing = 1500000.0;
576     else
577       *false_northing = 0.0;
578   }
579   else
580   {
581     if ((set_number % 2) == 0)
582       *false_northing =  500000.0;
583     else
584       *false_northing = 1000000.00;
585   }
586 } /* END OF Get_Grid_Values */
587 
588 
UTM_To_MGRS(long Zone,double Latitude,double Easting,double Northing,long Precision,char * MGRS)589 long UTM_To_MGRS (long Zone,
590                   double Latitude,
591                   double Easting,
592                   double Northing,
593                   long Precision,
594                   char *MGRS)
595 /*
596  * The function UTM_To_MGRS calculates an MGRS coordinate string
597  * based on the zone, latitude, easting and northing.
598  *
599  *    Zone      : Zone number             (input)
600  *    Latitude  : Latitude in radians     (input)
601  *    Easting   : Easting                 (input)
602  *    Northing  : Northing                (input)
603  *    Precision : Precision               (input)
604  *    MGRS      : MGRS coordinate string  (output)
605  */
606 { /* BEGIN UTM_To_MGRS */
607   double false_northing;      /* False northing for 3rd letter               */
608   double grid_easting;        /* Easting used to derive 2nd letter of MGRS   */
609   double grid_northing;       /* Northing used to derive 3rd letter of MGRS  */
610   long ltr2_low_value;        /* 2nd letter range - low number               */
611   long ltr2_high_value;       /* 2nd letter range - high number              */
612   int letters[MGRS_LETTERS];  /* Number location of 3 letters in alphabet    */
613   double divisor;
614   long error_code = MGRS_NO_ERROR;
615 
616 	/* Round easting and northing values */
617 	divisor = pow (10.0, (5 - Precision));
618 	Easting = Round_MGRS (Easting/divisor) * divisor;
619 	Northing = Round_MGRS (Northing/divisor) * divisor;
620 
621   Get_Grid_Values(Zone, &ltr2_low_value, &ltr2_high_value, &false_northing);
622 
623   error_code = Get_Latitude_Letter(Latitude, &letters[0]);
624 
625   if (!error_code)
626   {
627     grid_northing = Northing;
628     if (grid_northing == 1.e7)
629       grid_northing = grid_northing - 1.0;
630 
631     while (grid_northing >= TWOMIL)
632     {
633       grid_northing = grid_northing - TWOMIL;
634     }
635     grid_northing = grid_northing - false_northing;
636 
637     if (grid_northing < 0.0)
638       grid_northing = grid_northing + TWOMIL;
639 
640     letters[2] = (long)(grid_northing / ONEHT);
641     if (letters[2] > LETTER_H)
642       letters[2] = letters[2] + 1;
643 
644     if (letters[2] > LETTER_N)
645       letters[2] = letters[2] + 1;
646 
647     grid_easting = Easting;
648     if (((letters[0] == LETTER_V) && (Zone == 31)) && (grid_easting == 500000.0))
649       grid_easting = grid_easting - 1.0; /* SUBTRACT 1 METER */
650 
651     letters[1] = ltr2_low_value + ((long)(grid_easting / ONEHT) -1);
652     if ((ltr2_low_value == LETTER_J) && (letters[1] > LETTER_N))
653       letters[1] = letters[1] + 1;
654 
655     Make_MGRS_String (MGRS, Zone, letters, Easting, Northing, Precision);
656   }
657   return error_code;
658 } /* END UTM_To_MGRS */
659 
660 
Set_MGRS_Parameters(double a,double f,char * Ellipsoid_Code)661 long Set_MGRS_Parameters (double a,
662                           double f,
663                           char   *Ellipsoid_Code)
664 /*
665  * The function SET_MGRS_PARAMETERS receives the ellipsoid parameters and sets
666  * the corresponding state variables. If any errors occur, the error code(s)
667  * are returned by the function, otherwise MGRS_NO_ERROR is returned.
668  *
669  *   a                : Semi-major axis of ellipsoid in meters  (input)
670  *   f                : Flattening of ellipsoid					        (input)
671  *   Ellipsoid_Code   : 2-letter code for ellipsoid             (input)
672  */
673 { /* Set_MGRS_Parameters  */
674 
675   double inv_f = 1 / f;
676   long Error_Code = MGRS_NO_ERROR;
677 
678   if (a <= 0.0)
679   { /* Semi-major axis must be greater than zero */
680     Error_Code |= MGRS_A_ERROR;
681   }
682   if ((inv_f < 250) || (inv_f > 350))
683   { /* Inverse flattening must be between 250 and 350 */
684     Error_Code |= MGRS_INV_F_ERROR;
685   }
686   if (!Error_Code)
687   { /* no errors */
688     MGRS_a = a;
689     MGRS_f = f;
690     MGRS_recpf = inv_f;
691     strcpy (MGRS_Ellipsoid_Code, Ellipsoid_Code);
692   }
693   return (Error_Code);
694 }  /* Set_MGRS_Parameters  */
695 
696 
Get_MGRS_Parameters(double * a,double * f,char * Ellipsoid_Code)697 void Get_MGRS_Parameters (double *a,
698                           double *f,
699                           char* Ellipsoid_Code)
700 /*
701  * The function Get_MGRS_Parameters returns the current ellipsoid
702  * parameters.
703  *
704  *  a                : Semi-major axis of ellipsoid, in meters (output)
705  *  f                : Flattening of ellipsoid					       (output)
706  *  Ellipsoid_Code   : 2-letter code for ellipsoid             (output)
707  */
708 { /* Get_MGRS_Parameters */
709   *a = MGRS_a;
710   *f = MGRS_f;
711   strcpy (Ellipsoid_Code, MGRS_Ellipsoid_Code);
712   return;
713 } /* Get_MGRS_Parameters */
714 
715 #ifdef notdef
Convert_UTM_To_MGRS(long Zone,char Hemisphere,double Easting,double Northing,long Precision,char * MGRS)716 long Convert_UTM_To_MGRS (long Zone,
717                           char Hemisphere,
718                           double Easting,
719                           double Northing,
720                           long Precision,
721                           char* MGRS)
722 /*
723  * The function Convert_UTM_To_MGRS converts UTM (zone, easting, and
724  * northing) coordinates to an MGRS coordinate string, according to the
725  * current ellipsoid parameters.  If any errors occur, the error code(s)
726  * are returned by the function, otherwise MGRS_NO_ERROR is returned.
727  *
728  *    Zone       : UTM zone                         (input)
729  *    Hemisphere : North or South hemisphere        (input)
730  *    Easting    : Easting (X) in meters            (input)
731  *    Northing   : Northing (Y) in meters           (input)
732  *    Precision  : Precision level of MGRS string   (input)
733  *    MGRS       : MGRS coordinate string           (output)
734  */
735 { /* Convert_UTM_To_MGRS */
736   double latitude;           /* Latitude of UTM point */
737   double longitude;          /* Longitude of UTM point */
738   long temp_error = MGRS_NO_ERROR;
739   long error_code = MGRS_NO_ERROR;
740 
741   if ((Zone < 1) || (Zone > 60))
742     error_code |= MGRS_ZONE_ERROR;
743   if ((Hemisphere != 'S') && (Hemisphere != 'N'))
744     error_code |= MGRS_HEMISPHERE_ERROR;
745   if ((Easting < MIN_EASTING) || (Easting > MAX_EASTING))
746     error_code |= MGRS_EASTING_ERROR;
747   if ((Northing < MIN_NORTHING) || (Northing > MAX_NORTHING))
748     error_code |= MGRS_NORTHING_ERROR;
749   if ((Precision < 0) || (Precision > MAX_PRECISION))
750     error_code |= MGRS_PRECISION_ERROR;
751   if (!error_code)
752   {
753     Set_UTM_Parameters (MGRS_a, MGRS_f, 0);
754     temp_error = Convert_UTM_To_Geodetic (Zone, Hemisphere, Easting, Northing, &latitude, &longitude);
755 
756 	  /* Special check for rounding to (truncated) eastern edge of zone 31V */
757 	  if ((Zone == 31) && (latitude >= 56.0 * DEG_TO_RAD) && (latitude < 64.0 * DEG_TO_RAD) &&
758         (longitude >= 3.0 * DEG_TO_RAD))
759 	  { /* Reconvert to UTM zone 32 */
760       Set_UTM_Parameters (MGRS_a, MGRS_f, 32);
761       temp_error = Convert_Geodetic_To_UTM (latitude, longitude, &Zone, &Hemisphere, &Easting, &Northing);
762 	  }
763 
764 	  error_code = UTM_To_MGRS (Zone, latitude, Easting, Northing, Precision, MGRS);
765   }
766   return (error_code);
767 } /* Convert_UTM_To_MGRS */
768 #endif
769 
Convert_MGRS_To_UTM(char * MGRS,long * Zone,char * Hemisphere,double * Easting,double * Northing)770 long Convert_MGRS_To_UTM (char   *MGRS,
771                           long   *Zone,
772                           char   *Hemisphere,
773                           double *Easting,
774                           double *Northing)
775 /*
776  * The function Convert_MGRS_To_UTM converts an MGRS coordinate string
777  * to UTM projection (zone, hemisphere, easting and northing) coordinates
778  * according to the current ellipsoid parameters.  If any errors occur,
779  * the error code(s) are returned by the function, otherwise UTM_NO_ERROR
780  * is returned.
781  *
782  *    MGRS       : MGRS coordinate string           (input)
783  *    Zone       : UTM zone                         (output)
784  *    Hemisphere : North or South hemisphere        (output)
785  *    Easting    : Easting (X) in meters            (output)
786  *    Northing   : Northing (Y) in meters           (output)
787  */
788 { /* Convert_MGRS_To_UTM */
789   double scaled_min_northing;
790   double min_northing;
791   long ltr2_low_value;
792   long ltr2_high_value;
793   double false_northing;
794   double grid_easting;        /* Easting for 100,000 meter grid square      */
795   double grid_northing;       /* Northing for 100,000 meter grid square     */
796   long letters[MGRS_LETTERS];
797   long in_precision;
798 #ifdef notdef
799   double upper_lat_limit;     /* North latitude limits based on 1st letter  */
800   double lower_lat_limit;     /* South latitude limits based on 1st letter  */
801   double latitude = 0.0;
802   double longitude = 0.0;
803   double divisor = 1.0;
804 #endif
805   long error_code = MGRS_NO_ERROR;
806 
807   error_code = Break_MGRS_String (MGRS, Zone, letters, Easting, Northing, &in_precision);
808   if (!*Zone)
809     error_code |= MGRS_STRING_ERROR;
810   else
811   {
812     if (!error_code)
813     {
814       if ((letters[0] == LETTER_X) && ((*Zone == 32) || (*Zone == 34) || (*Zone == 36)))
815         error_code |= MGRS_STRING_ERROR;
816       else
817       {
818         if (letters[0] < LETTER_N)
819           *Hemisphere = 'S';
820         else
821           *Hemisphere = 'N';
822 
823         Get_Grid_Values(*Zone, &ltr2_low_value, &ltr2_high_value, &false_northing);
824 
825         /* Check that the second letter of the MGRS string is within
826          * the range of valid second letter values
827          * Also check that the third letter is valid */
828         if ((letters[1] < ltr2_low_value) || (letters[1] > ltr2_high_value) || (letters[2] > LETTER_V))
829           error_code |= MGRS_STRING_ERROR;
830 
831         if (!error_code)
832         {
833           grid_northing = (double)(letters[2]) * ONEHT + false_northing;
834           grid_easting = (double)((letters[1]) - ltr2_low_value + 1) * ONEHT;
835           if ((ltr2_low_value == LETTER_J) && (letters[1] > LETTER_O))
836             grid_easting = grid_easting - ONEHT;
837 
838           if (letters[2] > LETTER_O)
839             grid_northing = grid_northing - ONEHT;
840 
841           if (letters[2] > LETTER_I)
842             grid_northing = grid_northing - ONEHT;
843 
844           if (grid_northing >= TWOMIL)
845             grid_northing = grid_northing - TWOMIL;
846 
847           error_code = Get_Latitude_Band_Min_Northing(letters[0], &min_northing);
848           if (!error_code)
849           {
850             scaled_min_northing = min_northing;
851             while (scaled_min_northing >= TWOMIL)
852             {
853               scaled_min_northing = scaled_min_northing - TWOMIL;
854             }
855 
856             grid_northing = grid_northing - scaled_min_northing;
857             if (grid_northing < 0.0)
858               grid_northing = grid_northing + TWOMIL;
859 
860             grid_northing = min_northing + grid_northing;
861 
862             *Easting = grid_easting + *Easting;
863             *Northing = grid_northing + *Northing;
864 #ifdef notdef
865             /* check that point is within Zone Letter bounds */
866             error_code = Set_UTM_Parameters(MGRS_a,MGRS_f,*Zone);
867             if (!error_code)
868             {
869               error_code = Convert_UTM_To_Geodetic(*Zone,*Hemisphere,*Easting,*Northing,&latitude,&longitude);
870               if (!error_code)
871               {
872                 divisor = pow (10.0, in_precision);
873                 error_code = Get_Latitude_Range(letters[0], &upper_lat_limit, &lower_lat_limit);
874                 if (!error_code)
875                 {
876                   if (!(((lower_lat_limit - DEG_TO_RAD/divisor) <= latitude) && (latitude <= (upper_lat_limit + DEG_TO_RAD/divisor))))
877                     error_code |= MGRS_LAT_ERROR;
878                 }
879               }
880             }
881 #endif /* notdef */
882           }
883         }
884       }
885     }
886   }
887   return (error_code);
888 } /* Convert_MGRS_To_UTM */
889 
890 
Convert_UPS_To_MGRS(char Hemisphere,double Easting,double Northing,long Precision,char * MGRS)891 long Convert_UPS_To_MGRS (char   Hemisphere,
892                           double Easting,
893                           double Northing,
894                           long   Precision,
895                           char*  MGRS)
896 /*
897  *  The function Convert_UPS_To_MGRS converts UPS (hemisphere, easting,
898  *  and northing) coordinates to an MGRS coordinate string according to
899  *  the current ellipsoid parameters.  If any errors occur, the error
900  *  code(s) are returned by the function, otherwise UPS_NO_ERROR is
901  *  returned.
902  *
903  *    Hemisphere    : Hemisphere either 'N' or 'S'     (input)
904  *    Easting       : Easting/X in meters              (input)
905  *    Northing      : Northing/Y in meters             (input)
906  *    Precision     : Precision level of MGRS string   (input)
907  *    MGRS          : MGRS coordinate string           (output)
908  */
909 { /* Convert_UPS_To_MGRS */
910   double false_easting;       /* False easting for 2nd letter                 */
911   double false_northing;      /* False northing for 3rd letter                */
912   double grid_easting;        /* Easting used to derive 2nd letter of MGRS    */
913   double grid_northing;       /* Northing used to derive 3rd letter of MGRS   */
914   long ltr2_low_value;        /* 2nd letter range - low number                */
915   int letters[MGRS_LETTERS];  /* Number location of 3 letters in alphabet     */
916   double divisor;
917   int index = 0;
918   long error_code = MGRS_NO_ERROR;
919 
920   if ((Hemisphere != 'N') && (Hemisphere != 'S'))
921     error_code |= MGRS_HEMISPHERE_ERROR;
922   if ((Easting < MIN_EAST_NORTH) || (Easting > MAX_EAST_NORTH))
923     error_code |= MGRS_EASTING_ERROR;
924   if ((Northing < MIN_EAST_NORTH) || (Northing > MAX_EAST_NORTH))
925     error_code |= MGRS_NORTHING_ERROR;
926   if ((Precision < 0) || (Precision > MAX_PRECISION))
927     error_code |= MGRS_PRECISION_ERROR;
928   if (!error_code)
929   {
930     divisor = pow (10.0, (5 - Precision));
931     Easting = Round_MGRS (Easting/divisor) * divisor;
932     Northing = Round_MGRS (Northing/divisor) * divisor;
933 
934     if (Hemisphere == 'N')
935     {
936       if (Easting >= TWOMIL)
937         letters[0] = LETTER_Z;
938       else
939         letters[0] = LETTER_Y;
940 
941       index = letters[0] - 22;
942       ltr2_low_value = UPS_Constant_Table[index].ltr2_low_value;
943       false_easting = UPS_Constant_Table[index].false_easting;
944       false_northing = UPS_Constant_Table[index].false_northing;
945     }
946     else
947     {
948       if (Easting >= TWOMIL)
949         letters[0] = LETTER_B;
950       else
951         letters[0] = LETTER_A;
952 
953       ltr2_low_value = UPS_Constant_Table[letters[0]].ltr2_low_value;
954       false_easting = UPS_Constant_Table[letters[0]].false_easting;
955       false_northing = UPS_Constant_Table[letters[0]].false_northing;
956     }
957 
958     grid_northing = Northing;
959     grid_northing = grid_northing - false_northing;
960     letters[2] = (long)(grid_northing / ONEHT);
961 
962     if (letters[2] > LETTER_H)
963       letters[2] = letters[2] + 1;
964 
965     if (letters[2] > LETTER_N)
966       letters[2] = letters[2] + 1;
967 
968     grid_easting = Easting;
969     grid_easting = grid_easting - false_easting;
970     letters[1] = ltr2_low_value + ((long)(grid_easting / ONEHT));
971 
972     if (Easting < TWOMIL)
973     {
974       if (letters[1] > LETTER_L)
975         letters[1] = letters[1] + 3;
976 
977       if (letters[1] > LETTER_U)
978         letters[1] = letters[1] + 2;
979     }
980     else
981     {
982       if (letters[1] > LETTER_C)
983         letters[1] = letters[1] + 2;
984 
985       if (letters[1] > LETTER_H)
986         letters[1] = letters[1] + 1;
987 
988       if (letters[1] > LETTER_L)
989         letters[1] = letters[1] + 3;
990     }
991 
992     Make_MGRS_String (MGRS, 0, letters, Easting, Northing, Precision);
993   }
994   return (error_code);
995 } /* Convert_UPS_To_MGRS */
996 
997 
Convert_MGRS_To_UPS(char * MGRS,char * Hemisphere,double * Easting,double * Northing)998 long Convert_MGRS_To_UPS ( char   *MGRS,
999                            char   *Hemisphere,
1000                            double *Easting,
1001                            double *Northing)
1002 /*
1003  *  The function Convert_MGRS_To_UPS converts an MGRS coordinate string
1004  *  to UPS (hemisphere, easting, and northing) coordinates, according
1005  *  to the current ellipsoid parameters. If any errors occur, the error
1006  *  code(s) are returned by the function, otherwide UPS_NO_ERROR is returned.
1007  *
1008  *    MGRS          : MGRS coordinate string           (input)
1009  *    Hemisphere    : Hemisphere either 'N' or 'S'     (output)
1010  *    Easting       : Easting/X in meters              (output)
1011  *    Northing      : Northing/Y in meters             (output)
1012  */
1013 { /* Convert_MGRS_To_UPS */
1014   long ltr2_high_value;       /* 2nd letter range - high number             */
1015   long ltr3_high_value;       /* 3rd letter range - high number (UPS)       */
1016   long ltr2_low_value;        /* 2nd letter range - low number              */
1017   double false_easting;       /* False easting for 2nd letter               */
1018   double false_northing;      /* False northing for 3rd letter              */
1019   double grid_easting;        /* easting for 100,000 meter grid square      */
1020   double grid_northing;       /* northing for 100,000 meter grid square     */
1021   long zone;
1022   long letters[MGRS_LETTERS];
1023   long in_precision;
1024   int index = 0;
1025   long error_code = MGRS_NO_ERROR;
1026 
1027   error_code = Break_MGRS_String (MGRS, &zone, letters, Easting, Northing, &in_precision);
1028   if (zone)
1029     error_code |= MGRS_STRING_ERROR;
1030   else
1031   {
1032     if (!error_code)
1033     {
1034       if (letters[0] >= LETTER_Y)
1035       {
1036         *Hemisphere = 'N';
1037 
1038         index = letters[0] - 22;
1039         ltr2_low_value = UPS_Constant_Table[index].ltr2_low_value;
1040         ltr2_high_value = UPS_Constant_Table[index].ltr2_high_value;
1041         ltr3_high_value = UPS_Constant_Table[index].ltr3_high_value;
1042         false_easting = UPS_Constant_Table[index].false_easting;
1043         false_northing = UPS_Constant_Table[index].false_northing;
1044       }
1045       else
1046       {
1047         *Hemisphere = 'S';
1048 
1049         ltr2_low_value = UPS_Constant_Table[letters[0]].ltr2_low_value;
1050         ltr2_high_value = UPS_Constant_Table[letters[0]].ltr2_high_value;
1051         ltr3_high_value = UPS_Constant_Table[letters[0]].ltr3_high_value;
1052         false_easting = UPS_Constant_Table[letters[0]].false_easting;
1053         false_northing = UPS_Constant_Table[letters[0]].false_northing;
1054       }
1055 
1056       /* Check that the second letter of the MGRS string is within
1057        * the range of valid second letter values
1058        * Also check that the third letter is valid */
1059       if ((letters[1] < ltr2_low_value) || (letters[1] > ltr2_high_value) ||
1060           ((letters[1] == LETTER_D) || (letters[1] == LETTER_E) ||
1061           (letters[1] == LETTER_M) || (letters[1] == LETTER_N) ||
1062           (letters[1] == LETTER_V) || (letters[1] == LETTER_W)) ||
1063           (letters[2] > ltr3_high_value))
1064           error_code = MGRS_STRING_ERROR;
1065 
1066       if (!error_code)
1067       {
1068         grid_northing = (double)letters[2] * ONEHT + false_northing;
1069         if (letters[2] > LETTER_I)
1070           grid_northing = grid_northing - ONEHT;
1071 
1072         if (letters[2] > LETTER_O)
1073           grid_northing = grid_northing - ONEHT;
1074 
1075         grid_easting = (double)((letters[1]) - ltr2_low_value) * ONEHT + false_easting;
1076         if (ltr2_low_value != LETTER_A)
1077         {
1078           if (letters[1] > LETTER_L)
1079             grid_easting = grid_easting - 300000.0;
1080 
1081           if (letters[1] > LETTER_U)
1082             grid_easting = grid_easting - 200000.0;
1083         }
1084         else
1085         {
1086           if (letters[1] > LETTER_C)
1087             grid_easting = grid_easting - 200000.0;
1088 
1089           if (letters[1] > LETTER_I)
1090             grid_easting = grid_easting - ONEHT;
1091 
1092           if (letters[1] > LETTER_L)
1093             grid_easting = grid_easting - 300000.0;
1094         }
1095 
1096         *Easting = grid_easting + *Easting;
1097         *Northing = grid_northing + *Northing;
1098       }
1099     }
1100   }
1101   return (error_code);
1102 } /* Convert_MGRS_To_UPS */
1103 
1104 
1105 
1106