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