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