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