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