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