GeographicLib  1.38
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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-2014) <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 MGRS (periodic) row index and output is true row index. Band
343  // index is in [-10, 10) (as returned by LatitudeBand). Column index
344  // origin is easting = 100km. Returns maxutmSrow_ if irow and iband are
345  // incompatible. Row index origin is equator.
346 
347  // Estimate center row number for latitude band
348  // 90 deg = 100 tiles; 1 band = 8 deg = 100*8/90 tiles
349  real c = 100 * (8 * iband + 4)/real(90);
350  bool northp = iband >= 0;
351  // iband minrow maxrow
352  // -10 -90 -81
353  // -9 -80 -72
354  // -8 -71 -63
355  // -7 -63 -54
356  // -6 -54 -45
357  // -5 -45 -36
358  // -4 -36 -27
359  // -3 -27 -18
360  // -2 -18 -9
361  // -1 -9 -1
362  // 0 0 8
363  // 1 8 17
364  // 2 17 26
365  // 3 26 35
366  // 4 35 44
367  // 5 44 53
368  // 6 53 62
369  // 7 62 70
370  // 8 71 79
371  // 9 80 94
372  int
373  minrow = iband > -10 ?
374  int(floor(c - real(4.3) - real(0.1) * northp)) : -90,
375  maxrow = iband < 9 ?
376  int(floor(c + real(4.4) - real(0.1) * northp)) : 94,
377  baserow = (minrow + maxrow) / 2 - utmrowperiod_ / 2;
378  // Add maxutmSrow_ = 5 * utmrowperiod_ to ensure operand is positive
379  irow = (irow - baserow + maxutmSrow_) % utmrowperiod_ + baserow;
380  if (irow < minrow || irow > maxrow) {
381  // Northing = 71*100km and 80*100km intersect band boundaries
382  // The following deals with these special cases.
383  int
384  // Fold [-10,-1] -> [9,0]
385  sband = iband >= 0 ? iband : -iband - 1,
386  // Fold [-90,-1] -> [89,0]
387  srow = irow >= 0 ? irow : -irow - 1,
388  // Fold [4,7] -> [3,0]
389  scol = icol < 4 ? icol : -icol + 7;
390  if ( ! ( (srow == 70 && sband == 8 && scol >= 2) ||
391  (srow == 71 && sband == 7 && scol <= 2) ||
392  (srow == 79 && sband == 9 && scol >= 1) ||
393  (srow == 80 && sband == 8 && scol <= 1) ) )
394  irow = maxutmSrow_;
395  }
396  return irow;
397  }
398 
399 } // namespace GeographicLib
static T NaN()
Definition: Math.hpp:460
GeographicLib::Math::real real
Definition: GeodSolve.cpp:40
Header for GeographicLib::Utility class.
static bool isnan(T x)
Definition: Math.hpp:477
Header for GeographicLib::MGRS class.
static std::string str(T x, int p=-1)
Definition: Utility.hpp:266
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:362
static int lookup(const std::string &s, char c)
Definition: Utility.hpp:382
static void Forward(int zone, bool northp, real x, real y, int prec, std::string &mgrs)
Definition: MGRS.cpp:123