GeographicLib  2.0
MGRS.cpp
Go to the documentation of this file.
1 /**
2  * \file MGRS.cpp
3  * \brief Implementation for GeographicLib::MGRS class
4  *
5  * Copyright (c) Charles Karney (2008-2022) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
8  **********************************************************************/
9 
10 #include <GeographicLib/MGRS.hpp>
12 
13 namespace GeographicLib {
14 
15  using namespace std;
16 
17  const char* const MGRS::hemispheres_ = "SN";
18  const char* const MGRS::utmcols_[] = { "ABCDEFGH", "JKLMNPQR", "STUVWXYZ" };
19  const char* const MGRS::utmrow_ = "ABCDEFGHJKLMNPQRSTUV";
20  const char* const MGRS::upscols_[] =
21  { "JKLPQRSTUXYZ", "ABCFGHJKLPQR", "RSTUXYZ", "ABCFGHJ" };
22  const char* const MGRS::upsrows_[] =
23  { "ABCDEFGHJKLMNPQRSTUVWXYZ", "ABCDEFGHJKLMNP" };
24  const char* const MGRS::latband_ = "CDEFGHJKLMNPQRSTUVWX";
25  const char* const MGRS::upsband_ = "ABYZ";
26  const char* const MGRS::digits_ = "0123456789";
27 
28  const int MGRS::mineasting_[] =
29  { minupsSind_, minupsNind_, minutmcol_, minutmcol_ };
30  const int MGRS::maxeasting_[] =
31  { maxupsSind_, maxupsNind_, maxutmcol_, maxutmcol_ };
32  const int MGRS::minnorthing_[] =
33  { minupsSind_, minupsNind_,
34  minutmSrow_, minutmSrow_ - (maxutmSrow_ - minutmNrow_) };
35  const int MGRS::maxnorthing_[] =
36  { maxupsSind_, maxupsNind_,
37  maxutmNrow_ + (maxutmSrow_ - minutmNrow_), maxutmNrow_ };
38 
39  void MGRS::Forward(int zone, bool northp, real x, real y, real lat,
40  int prec, std::string& mgrs) {
41  using std::isnan; // Needed for Centos 7, ubuntu 14
42  // The smallest angle s.t., 90 - angeps() < 90 (approx 50e-12 arcsec)
43  // 7 = ceil(log_2(90))
44  static const real angeps = ldexp(real(1), -(Math::digits() - 7));
45  if (zone == UTMUPS::INVALID ||
46  isnan(x) || isnan(y) || isnan(lat)) {
47  mgrs = "INVALID";
48  return;
49  }
50  bool utmp = zone != 0;
51  CheckCoords(utmp, northp, x, y);
52  if (!(zone >= UTMUPS::MINZONE && zone <= UTMUPS::MAXZONE))
53  throw GeographicErr("Zone " + Utility::str(zone) + " not in [0,60]");
54  if (!(prec >= -1 && prec <= maxprec_))
55  throw GeographicErr("MGRS precision " + Utility::str(prec)
56  + " not in [-1, "
57  + Utility::str(int(maxprec_)) + "]");
58  // Fixed char array for accumulating string. Allow space for zone, 3 block
59  // letters, easting + northing. Don't need to allow for terminating null.
60  char mgrs1[2 + 3 + 2 * maxprec_];
61  int
62  zone1 = zone - 1,
63  z = utmp ? 2 : 0,
64  mlen = z + 3 + 2 * prec;
65  if (utmp) {
66  mgrs1[0] = digits_[ zone / base_ ];
67  mgrs1[1] = digits_[ zone % base_ ];
68  // This isn't necessary...! Keep y non-neg
69  // if (!northp) y -= maxutmSrow_ * tile_;
70  }
71  // The C++ standard mandates 64 bits for long long. But
72  // check, to make sure.
73  static_assert(numeric_limits<long long>::digits >= 44,
74  "long long not wide enough to store 10e12");
75  // Guard against floor(x * mult_) being computed incorrectly on some
76  // platforms. The problem occurs when x * mult_ is held in extended
77  // precision and floor is inlined. This causes tests GeoConvert1[678] to
78  // fail. Problem reported and diagnosed by Thorkil Naur with g++ 10.2.0
79  // under Cygwin.
80  GEOGRAPHICLIB_VOLATILE real xx = x * mult_;
81  GEOGRAPHICLIB_VOLATILE real yy = y * mult_;
82  long long
83  ix = (long long)(floor(xx)),
84  iy = (long long)(floor(yy)),
85  m = (long long)(mult_) * (long long)(tile_);
86  int xh = int(ix / m), yh = int(iy / m);
87  if (utmp) {
88  int
89  // Correct fuzziness in latitude near equator
90  iband = fabs(lat) < angeps ? (northp ? 0 : -1) : LatitudeBand(lat),
91  icol = xh - minutmcol_,
92  irow = UTMRow(iband, icol, yh % utmrowperiod_);
93  if (irow != yh - (northp ? minutmNrow_ : maxutmSrow_))
94  throw GeographicErr("Latitude " + Utility::str(lat)
95  + " is inconsistent with UTM coordinates");
96  mgrs1[z++] = latband_[10 + iband];
97  mgrs1[z++] = utmcols_[zone1 % 3][icol];
98  mgrs1[z++] = utmrow_[(yh + (zone1 & 1 ? utmevenrowshift_ : 0))
99  % utmrowperiod_];
100  } else {
101  bool eastp = xh >= upseasting_;
102  int iband = (northp ? 2 : 0) + (eastp ? 1 : 0);
103  mgrs1[z++] = upsband_[iband];
104  mgrs1[z++] = upscols_[iband][xh - (eastp ? upseasting_ :
105  (northp ? minupsNind_ :
106  minupsSind_))];
107  mgrs1[z++] = upsrows_[northp][yh - (northp ? minupsNind_ : minupsSind_)];
108  }
109  if (prec > 0) {
110  ix -= m * xh; iy -= m * yh;
111  long long d = (long long)(pow(real(base_), maxprec_ - prec));
112  ix /= d; iy /= d;
113  for (int c = prec; c--;) {
114  mgrs1[z + c ] = digits_[ix % base_]; ix /= base_;
115  mgrs1[z + c + prec] = digits_[iy % base_]; iy /= base_;
116  }
117  }
118  mgrs.resize(mlen);
119  copy(mgrs1, mgrs1 + mlen, mgrs.begin());
120  }
121 
122  void MGRS::Forward(int zone, bool northp, real x, real y,
123  int prec, std::string& mgrs) {
124  real lat, lon;
125  if (zone > 0) {
126  // Does a rough estimate for latitude determine the latitude band?
127  real ys = northp ? y : y - utmNshift_;
128  // A cheap calculation of the latitude which results in an "allowed"
129  // latitude band would be
130  // lat = ApproxLatitudeBand(ys) * 8 + 4;
131  //
132  // Here we do a more careful job using the band letter corresponding to
133  // the actual latitude.
134  ys /= tile_;
135  if (fabs(ys) < 1)
136  lat = real(0.9) * ys; // accurate enough estimate near equator
137  else {
138  real
139  // The poleward bound is a fit from above of lat(x,y)
140  // for x = 500km and y = [0km, 950km]
141  latp = real(0.901) * ys + (ys > 0 ? 1 : -1) * real(0.135),
142  // The equatorward bound is a fit from below of lat(x,y)
143  // for x = 900km and y = [0km, 950km]
144  late = real(0.902) * ys * (1 - real(1.85e-6) * ys * ys);
145  if (LatitudeBand(latp) == LatitudeBand(late))
146  lat = latp;
147  else
148  // bounds straddle a band boundary so need to compute lat accurately
149  UTMUPS::Reverse(zone, northp, x, y, lat, lon);
150  }
151  } else
152  // Latitude isn't needed for UPS specs or for INVALID
153  lat = 0;
154  Forward(zone, northp, x, y, lat, prec, mgrs);
155  }
156 
157  void MGRS::Reverse(const string& mgrs,
158  int& zone, bool& northp, real& x, real& y,
159  int& prec, bool centerp) {
160  int
161  p = 0,
162  len = int(mgrs.length());
163  if (len >= 3 &&
164  toupper(mgrs[0]) == 'I' &&
165  toupper(mgrs[1]) == 'N' &&
166  toupper(mgrs[2]) == 'V') {
167  zone = UTMUPS::INVALID;
168  northp = false;
169  x = y = Math::NaN();
170  prec = -2;
171  return;
172  }
173  int zone1 = 0;
174  while (p < len) {
175  int i = Utility::lookup(digits_, mgrs[p]);
176  if (i < 0)
177  break;
178  zone1 = 10 * zone1 + i;
179  ++p;
180  }
181  if (p > 0 && !(zone1 >= UTMUPS::MINUTMZONE && zone1 <= UTMUPS::MAXUTMZONE))
182  throw GeographicErr("Zone " + Utility::str(zone1) + " not in [1,60]");
183  if (p > 2)
184  throw GeographicErr("More than 2 digits at start of MGRS "
185  + mgrs.substr(0, p));
186  if (len - p < 1)
187  throw GeographicErr("MGRS string too short " + mgrs);
188  bool utmp = zone1 != UTMUPS::UPS;
189  int zonem1 = zone1 - 1;
190  const char* band = utmp ? latband_ : upsband_;
191  int iband = Utility::lookup(band, mgrs[p++]);
192  if (iband < 0)
193  throw GeographicErr("Band letter " + Utility::str(mgrs[p-1]) + " not in "
194  + (utmp ? "UTM" : "UPS") + " set " + band);
195  bool northp1 = iband >= (utmp ? 10 : 2);
196  if (p == len) { // Grid zone only (ignore centerp)
197  // Approx length of a degree of meridian arc in units of tile.
198  real deg = real(utmNshift_) / (Math::qd * tile_);
199  zone = zone1;
200  northp = northp1;
201  if (utmp) {
202  // Pick central meridian except for 31V
203  x = ((zone == 31 && iband == 17) ? 4 : 5) * tile_;
204  // Pick center of 8deg latitude bands
205  y = floor(8 * (iband - real(9.5)) * deg + real(0.5)) * tile_
206  + (northp ? 0 : utmNshift_);
207  } else {
208  // Pick point at lat 86N or 86S
209  x = ((iband & 1 ? 1 : -1) * floor(4 * deg + real(0.5))
210  + upseasting_) * tile_;
211  // Pick point at lon 90E or 90W.
212  y = upseasting_ * tile_;
213  }
214  prec = -1;
215  return;
216  } else if (len - p < 2)
217  throw GeographicErr("Missing row letter in " + mgrs);
218  const char* col = utmp ? utmcols_[zonem1 % 3] : upscols_[iband];
219  const char* row = utmp ? utmrow_ : upsrows_[northp1];
220  int icol = Utility::lookup(col, mgrs[p++]);
221  if (icol < 0)
222  throw GeographicErr("Column letter " + Utility::str(mgrs[p-1])
223  + " not in "
224  + (utmp ? "zone " + mgrs.substr(0, p-2) :
225  "UPS band " + Utility::str(mgrs[p-2]))
226  + " set " + col );
227  int irow = Utility::lookup(row, mgrs[p++]);
228  if (irow < 0)
229  throw GeographicErr("Row letter " + Utility::str(mgrs[p-1]) + " not in "
230  + (utmp ? "UTM" :
231  "UPS " + Utility::str(hemispheres_[northp1]))
232  + " set " + row);
233  if (utmp) {
234  if (zonem1 & 1)
235  irow = (irow + utmrowperiod_ - utmevenrowshift_) % utmrowperiod_;
236  iband -= 10;
237  irow = UTMRow(iband, icol, irow);
238  if (irow == maxutmSrow_)
239  throw GeographicErr("Block " + mgrs.substr(p-2, 2)
240  + " not in zone/band " + mgrs.substr(0, p-2));
241 
242  irow = northp1 ? irow : irow + 100;
243  icol = icol + minutmcol_;
244  } else {
245  bool eastp = iband & 1;
246  icol += eastp ? upseasting_ : (northp1 ? minupsNind_ : minupsSind_);
247  irow += northp1 ? minupsNind_ : minupsSind_;
248  }
249  int prec1 = (len - p)/2;
250  real
251  unit = 1,
252  x1 = icol,
253  y1 = irow;
254  for (int i = 0; i < prec1; ++i) {
255  unit *= base_;
256  int
257  ix = Utility::lookup(digits_, mgrs[p + i]),
258  iy = Utility::lookup(digits_, mgrs[p + i + prec1]);
259  if (ix < 0 || iy < 0)
260  throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p));
261  x1 = base_ * x1 + ix;
262  y1 = base_ * y1 + iy;
263  }
264  if ((len - p) % 2) {
265  if (Utility::lookup(digits_, mgrs[len - 1]) < 0)
266  throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p));
267  else
268  throw GeographicErr("Not an even number of digits in "
269  + mgrs.substr(p));
270  }
271  if (prec1 > maxprec_)
272  throw GeographicErr("More than " + Utility::str(2*maxprec_)
273  + " digits in " + mgrs.substr(p));
274  if (centerp) {
275  unit *= 2; x1 = 2 * x1 + 1; y1 = 2 * y1 + 1;
276  }
277  zone = zone1;
278  northp = northp1;
279  x = (tile_ * x1) / unit;
280  y = (tile_ * y1) / unit;
281  prec = prec1;
282  }
283 
284  void MGRS::CheckCoords(bool utmp, bool& northp, real& x, real& y) {
285  // Limits are all multiples of 100km and are all closed on the lower end
286  // and open on the upper end -- and this is reflected in the error
287  // messages. However if a coordinate lies on the excluded upper end (e.g.,
288  // after rounding), it is shifted down by eps. This also folds UTM
289  // northings to the correct N/S hemisphere.
290 
291  // The smallest length s.t., 1.0e7 - eps() < 1.0e7 (approx 1.9 nm)
292  // 25 = ceil(log_2(2e7)) -- use half circumference here because
293  // northing 195e5 is a legal in the "southern" hemisphere.
294  static const real eps = ldexp(real(1), -(Math::digits() - 25));
295  int
296  ix = int(floor(x / tile_)),
297  iy = int(floor(y / tile_)),
298  ind = (utmp ? 2 : 0) + (northp ? 1 : 0);
299  if (! (ix >= mineasting_[ind] && ix < maxeasting_[ind]) ) {
300  if (ix == maxeasting_[ind] && x == maxeasting_[ind] * tile_)
301  x -= eps;
302  else
303  throw GeographicErr("Easting " + Utility::str(int(floor(x/1000)))
304  + "km not in MGRS/"
305  + (utmp ? "UTM" : "UPS") + " range for "
306  + (northp ? "N" : "S" ) + " hemisphere ["
307  + Utility::str(mineasting_[ind]*tile_/1000)
308  + "km, "
309  + Utility::str(maxeasting_[ind]*tile_/1000)
310  + "km)");
311  }
312  if (! (iy >= minnorthing_[ind] && iy < maxnorthing_[ind]) ) {
313  if (iy == maxnorthing_[ind] && y == maxnorthing_[ind] * tile_)
314  y -= eps;
315  else
316  throw GeographicErr("Northing " + Utility::str(int(floor(y/1000)))
317  + "km not in MGRS/"
318  + (utmp ? "UTM" : "UPS") + " range for "
319  + (northp ? "N" : "S" ) + " hemisphere ["
320  + Utility::str(minnorthing_[ind]*tile_/1000)
321  + "km, "
322  + Utility::str(maxnorthing_[ind]*tile_/1000)
323  + "km)");
324  }
325 
326  // Correct the UTM northing and hemisphere if necessary
327  if (utmp) {
328  if (northp && iy < minutmNrow_) {
329  northp = false;
330  y += utmNshift_;
331  } else if (!northp && iy >= maxutmSrow_) {
332  if (y == maxutmSrow_ * tile_)
333  // If on equator retain S hemisphere
334  y -= eps;
335  else {
336  northp = true;
337  y -= utmNshift_;
338  }
339  }
340  }
341  }
342 
343  int MGRS::UTMRow(int iband, int icol, int irow) {
344  // Input is iband = band index in [-10, 10) (as returned by LatitudeBand),
345  // icol = column index in [0,8) with origin of easting = 100km, and irow =
346  // periodic row index in [0,20) with origin = equator. Output is true row
347  // index in [-90, 95). Returns maxutmSrow_ = 100, if irow and iband are
348  // incompatible.
349 
350  // Estimate center row number for latitude band
351  // 90 deg = 100 tiles; 1 band = 8 deg = 100*8/90 tiles
352  real c = 100 * (8 * iband + 4) / real(Math::qd);
353  bool northp = iband >= 0;
354  // These are safe bounds on the rows
355  // iband minrow maxrow
356  // -10 -90 -81
357  // -9 -80 -72
358  // -8 -71 -63
359  // -7 -63 -54
360  // -6 -54 -45
361  // -5 -45 -36
362  // -4 -36 -27
363  // -3 -27 -18
364  // -2 -18 -9
365  // -1 -9 -1
366  // 0 0 8
367  // 1 8 17
368  // 2 17 26
369  // 3 26 35
370  // 4 35 44
371  // 5 44 53
372  // 6 53 62
373  // 7 62 70
374  // 8 71 79
375  // 9 80 94
376  int
377  minrow = iband > -10 ?
378  int(floor(c - real(4.3) - real(0.1) * northp)) : -90,
379  maxrow = iband < 9 ?
380  int(floor(c + real(4.4) - real(0.1) * northp)) : 94,
381  baserow = (minrow + maxrow) / 2 - utmrowperiod_ / 2;
382  // Offset irow by the multiple of utmrowperiod_ which brings it as close as
383  // possible to the center of the latitude band, (minrow + maxrow) / 2.
384  // (Add maxutmSrow_ = 5 * utmrowperiod_ to ensure operand is positive.)
385  irow = (irow - baserow + maxutmSrow_) % utmrowperiod_ + baserow;
386  if (!( irow >= minrow && irow <= maxrow )) {
387  // Outside the safe bounds, so need to check...
388  // Northing = 71e5 and 80e5 intersect band boundaries
389  // y = 71e5 in scol = 2 (x = [3e5,4e5] and x = [6e5,7e5])
390  // y = 80e5 in scol = 1 (x = [2e5,3e5] and x = [7e5,8e5])
391  // This holds for all the ellipsoids given in NGA.SIG.0012_2.0.0_UTMUPS.
392  // The following deals with these special cases.
393  int
394  // Fold [-10,-1] -> [9,0]
395  sband = iband >= 0 ? iband : -iband - 1,
396  // Fold [-90,-1] -> [89,0]
397  srow = irow >= 0 ? irow : -irow - 1,
398  // Fold [4,7] -> [3,0]
399  scol = icol < 4 ? icol : -icol + 7;
400  // For example, the safe rows for band 8 are 71 - 79. However row 70 is
401  // allowed if scol = [2,3] and row 80 is allowed if scol = [0,1].
402  if ( ! ( (srow == 70 && sband == 8 && scol >= 2) ||
403  (srow == 71 && sband == 7 && scol <= 2) ||
404  (srow == 79 && sband == 9 && scol >= 1) ||
405  (srow == 80 && sband == 8 && scol <= 1) ) )
406  irow = maxutmSrow_;
407  }
408  return irow;
409  }
410 
411  void MGRS::Check() {
412  real lat, lon, x, y, t = tile_; int zone; bool northp;
413  UTMUPS::Reverse(31, true , 1*t, 0*t, lat, lon);
414  if (!( lon < 0 ))
415  throw GeographicErr("MGRS::Check: equator coverage failure");
416  UTMUPS::Reverse(31, true , 1*t, 95*t, lat, lon);
417  if (!( lat > 84 ))
418  throw GeographicErr("MGRS::Check: UTM doesn't reach latitude = 84");
419  UTMUPS::Reverse(31, false, 1*t, 10*t, lat, lon);
420  if (!( lat < -80 ))
421  throw GeographicErr("MGRS::Check: UTM doesn't reach latitude = -80");
422  UTMUPS::Forward(56, 3, zone, northp, x, y, 32);
423  if (!( x > 1*t ))
424  throw GeographicErr("MGRS::Check: Norway exception creates a gap");
425  UTMUPS::Forward(72, 21, zone, northp, x, y, 35);
426  if (!( x > 1*t ))
427  throw GeographicErr("MGRS::Check: Svalbard exception creates a gap");
428  UTMUPS::Reverse(0, true , 20*t, 13*t, lat, lon);
429  if (!( lat < 84 ))
430  throw
431  GeographicErr("MGRS::Check: North UPS doesn't reach latitude = 84");
432  UTMUPS::Reverse(0, false, 20*t, 8*t, lat, lon);
433  if (!( lat > -80 ))
434  throw
435  GeographicErr("MGRS::Check: South UPS doesn't reach latitude = -80");
436  // Entries are [band, x, y] either side of the band boundaries. Units for
437  // x, y are t = 100km.
438  const short tab[] = {
439  0, 5, 0, 0, 9, 0, // south edge of band 0
440  0, 5, 8, 0, 9, 8, // north edge of band 0
441  1, 5, 9, 1, 9, 9, // south edge of band 1
442  1, 5, 17, 1, 9, 17, // north edge of band 1
443  2, 5, 18, 2, 9, 18, // etc.
444  2, 5, 26, 2, 9, 26,
445  3, 5, 27, 3, 9, 27,
446  3, 5, 35, 3, 9, 35,
447  4, 5, 36, 4, 9, 36,
448  4, 5, 44, 4, 9, 44,
449  5, 5, 45, 5, 9, 45,
450  5, 5, 53, 5, 9, 53,
451  6, 5, 54, 6, 9, 54,
452  6, 5, 62, 6, 9, 62,
453  7, 5, 63, 7, 9, 63,
454  7, 5, 70, 7, 7, 70, 7, 7, 71, 7, 9, 71, // y = 71t crosses boundary
455  8, 5, 71, 8, 6, 71, 8, 6, 72, 8, 9, 72, // between bands 7 and 8.
456  8, 5, 79, 8, 8, 79, 8, 8, 80, 8, 9, 80, // y = 80t crosses boundary
457  9, 5, 80, 9, 7, 80, 9, 7, 81, 9, 9, 81, // between bands 8 and 9.
458  9, 5, 95, 9, 9, 95, // north edge of band 9
459  };
460  const int bandchecks = sizeof(tab) / (3 * sizeof(short));
461  for (int i = 0; i < bandchecks; ++i) {
462  UTMUPS::Reverse(38, true, tab[3*i+1]*t, tab[3*i+2]*t, lat, lon);
463  if (!( LatitudeBand(lat) == tab[3*i+0] ))
464  throw GeographicErr("MGRS::Check: Band error, b = " +
465  Utility::str(tab[3*i+0]) + ", x = " +
466  Utility::str(tab[3*i+1]) + "00km, y = " +
467  Utility::str(tab[3*i+2]) + "00km");
468  }
469  }
470 
471 } // namespace GeographicLib
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
Header for GeographicLib::MGRS class.
#define GEOGRAPHICLIB_VOLATILE
Definition: Math.hpp:58
Header for GeographicLib::Utility class.
Exception handling for GeographicLib.
Definition: Constants.hpp:316
static void Reverse(const std::string &mgrs, int &zone, bool &northp, real &x, real &y, int &prec, bool centerp=true)
Definition: MGRS.cpp:157
static void Check()
Definition: MGRS.cpp:411
static void Forward(int zone, bool northp, real x, real y, int prec, std::string &mgrs)
Definition: MGRS.cpp:122
static int digits()
Definition: Math.cpp:26
static T NaN()
Definition: Math.cpp:250
@ qd
degrees per quarter turn
Definition: Math.hpp:141
static void Forward(real lat, real lon, int &zone, bool &northp, real &x, real &y, real &gamma, real &k, int setzone=STANDARD, bool mgrslimits=false)
Definition: UTMUPS.cpp:65
static void Reverse(int zone, bool northp, real x, real y, real &lat, real &lon, real &gamma, real &k, bool mgrslimits=false)
Definition: UTMUPS.cpp:119
static int lookup(const std::string &s, char c)
Definition: Utility.cpp:160
static std::string str(T x, int p=-1)
Definition: Utility.hpp:161
Namespace for GeographicLib.
Definition: Accumulator.cpp:12