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