GeographicLib  1.36
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-2012) <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 Math::real MGRS::eps_ =
18  // 25 = ceil(log_2(2e7)) -- use half circumference here because northing
19  // 195e5 is a legal in the "southern" hemisphere.
20  pow(real(0.5), numeric_limits<real>::digits - 25);
21  const Math::real MGRS::angeps_ =
22  // 7 = ceil(log_2(90))
23  pow(real(0.5), numeric_limits<real>::digits - 7);
24  const string MGRS::hemispheres_ = "SN";
25  const string MGRS::utmcols_[3] = { "ABCDEFGH", "JKLMNPQR", "STUVWXYZ" };
26  const string MGRS::utmrow_ = "ABCDEFGHJKLMNPQRSTUV";
27  const string MGRS::upscols_[4] =
28  { "JKLPQRSTUXYZ", "ABCFGHJKLPQR", "RSTUXYZ", "ABCFGHJ" };
29  const string MGRS::upsrows_[2] =
30  { "ABCDEFGHJKLMNPQRSTUVWXYZ", "ABCDEFGHJKLMNP" };
31  const string MGRS::latband_ = "CDEFGHJKLMNPQRSTUVWX";
32  const string MGRS::upsband_ = "ABYZ";
33  const string MGRS::digits_ = "0123456789";
34 
35  const int MGRS::mineasting_[4] =
36  { minupsSind_, minupsNind_, minutmcol_, minutmcol_ };
37  const int MGRS::maxeasting_[4] =
38  { maxupsSind_, maxupsNind_, maxutmcol_, maxutmcol_ };
39  const int MGRS::minnorthing_[4] =
40  { minupsSind_, minupsNind_,
41  minutmSrow_, minutmSrow_ - (maxutmSrow_ - minutmNrow_) };
42  const int MGRS::maxnorthing_[4] =
43  { maxupsSind_, maxupsNind_,
44  maxutmNrow_ + (maxutmSrow_ - minutmNrow_), maxutmNrow_ };
45 
46  void MGRS::Forward(int zone, bool northp, real x, real y, real lat,
47  int prec, std::string& mgrs) {
48  if (zone == UTMUPS::INVALID ||
49  Math::isnan(x) || Math::isnan(y) || Math::isnan(lat)) {
50  mgrs = "INVALID";
51  return;
52  }
53  bool utmp = zone != 0;
54  CheckCoords(utmp, northp, x, y);
55  if (!(zone >= UTMUPS::MINZONE && zone <= UTMUPS::MAXZONE))
56  throw GeographicErr("Zone " + Utility::str(zone) + " not in [0,60]");
57  if (!(prec >= -1 && prec <= maxprec_))
58  throw GeographicErr("MGRS precision " + Utility::str(prec)
59  + " not in [-1, "
60  + Utility::str(int(maxprec_)) + "]");
61  // Fixed char array for accumulating string. Allow space for zone, 3 block
62  // letters, easting + northing. Don't need to allow for terminating null.
63  char mgrs1[2 + 3 + 2 * maxprec_];
64  int
65  zone1 = zone - 1,
66  z = utmp ? 2 : 0,
67  mlen = z + 3 + 2 * prec;
68  if (utmp) {
69  mgrs1[0] = digits_[ zone / base_ ];
70  mgrs1[1] = digits_[ zone % base_ ];
71  // This isn't necessary...! Keep y non-neg
72  // if (!northp) y -= maxutmSrow_ * tile_;
73  }
74  int
75  xh = int(floor(x)) / tile_,
76  yh = int(floor(y)) / tile_;
77  real
78  xf = x - tile_ * xh,
79  yf = y - tile_ * yh;
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_ : minupsSind_))];
99  mgrs1[z++] = upsrows_[northp][yh - (northp ? minupsNind_ : minupsSind_)];
100  }
101  if (prec > 0) {
102  real mult = pow(real(base_), max(tilelevel_ - prec, 0));
103  int
104  ix = int(floor(xf / mult)),
105  iy = int(floor(yf / mult));
106  for (int c = min(prec, int(tilelevel_)); c--;) {
107  mgrs1[z + c] = digits_[ ix % base_ ];
108  ix /= base_;
109  mgrs1[z + c + prec] = digits_[ iy % base_ ];
110  iy /= base_;
111  }
112  if (prec > tilelevel_) {
113  xf -= floor(xf / mult);
114  yf -= floor(yf / mult);
115  mult = pow(real(base_), prec - tilelevel_);
116  ix = int(floor(xf * mult));
117  iy = int(floor(yf * mult));
118  for (int c = prec - tilelevel_; c--;) {
119  mgrs1[z + c + tilelevel_] = digits_[ ix % base_ ];
120  ix /= base_;
121  mgrs1[z + c + tilelevel_ + prec] = digits_[ iy % base_ ];
122  iy /= base_;
123  }
124  }
125  }
126  mgrs.resize(mlen);
127  copy(mgrs1, mgrs1 + mlen, mgrs.begin());
128  }
129 
130  void MGRS::Forward(int zone, bool northp, real x, real y,
131  int prec, std::string& mgrs) {
132  real lat, lon;
133  if (zone > 0) {
134  // Does a rough estimate for latitude determine the latitude band?
135  real ys = northp ? y : y - utmNshift_;
136  // A cheap calculation of the latitude which results in an "allowed"
137  // latitude band would be
138  // lat = ApproxLatitudeBand(ys) * 8 + 4;
139  //
140  // Here we do a more careful job using the band letter corresponding to
141  // the actual latitude.
142  ys /= tile_;
143  if (abs(ys) < 1)
144  lat = 0.9 * ys; // accurate enough estimate near equator
145  else {
146  real
147  // The poleward bound is a fit from above of lat(x,y)
148  // for x = 500km and y = [0km, 950km]
149  latp = real(0.901) * ys + (ys > 0 ? 1 : -1) * real(0.135),
150  // The equatorward bound is a fit from below of lat(x,y)
151  // for x = 900km and y = [0km, 950km]
152  late = real(0.902) * ys * (1 - real(1.85e-6) * ys * ys);
153  if (LatitudeBand(latp) == LatitudeBand(late))
154  lat = latp;
155  else
156  // bounds straddle a band boundary so need to compute lat accurately
157  UTMUPS::Reverse(zone, northp, x, y, lat, lon);
158  }
159  } else
160  // Latitude isn't needed for UPS specs or for INVALID
161  lat = 0;
162  Forward(zone, northp, x, y, lat, prec, mgrs);
163  }
164 
165  void MGRS::Reverse(const std::string& mgrs,
166  int& zone, bool& northp, real& x, real& y,
167  int& prec, bool centerp) {
168  int
169  p = 0,
170  len = int(mgrs.size());
171  if (len >= 3 &&
172  toupper(mgrs[0]) == 'I' &&
173  toupper(mgrs[1]) == 'N' &&
174  toupper(mgrs[2]) == 'V') {
175  zone = UTMUPS::INVALID;
176  northp = false;
177  x = y = Math::NaN<real>();
178  prec = -2;
179  return;
180  }
181  int zone1 = 0;
182  while (p < len) {
183  int i = Utility::lookup(digits_, mgrs[p]);
184  if (i < 0)
185  break;
186  zone1 = 10 * zone1 + i;
187  ++p;
188  }
189  if (p > 0 && !(zone1 >= UTMUPS::MINUTMZONE && zone1 <= UTMUPS::MAXUTMZONE))
190  throw GeographicErr("Zone " + Utility::str(zone1) + " not in [1,60]");
191  if (p > 2)
192  throw GeographicErr("More than 2 digits_ at start of MGRS "
193  + mgrs.substr(0, p));
194  if (len - p < 1)
195  throw GeographicErr("MGRS string too short " + mgrs);
196  bool utmp = zone1 != UTMUPS::UPS;
197  int zonem1 = zone1 - 1;
198  const string& band = utmp ? latband_ : upsband_;
199  int iband = Utility::lookup(band, mgrs[p++]);
200  if (iband < 0)
201  throw GeographicErr("Band letter " + Utility::str(mgrs[p-1]) + " not in "
202  + (utmp ? "UTM" : "UPS") + " set " + band);
203  bool northp1 = iband >= (utmp ? 10 : 2);
204  if (p == len) { // Grid zone only (ignore centerp)
205  // Approx length of a degree of meridian arc in units of tile.
206  real deg = real(utmNshift_) / (90 * tile_);
207  zone = zone1;
208  northp = northp1;
209  if (utmp) {
210  // Pick central meridian except for 31V
211  x = ((zone == 31 && iband == 17) ? 4 : 5) * tile_;
212  // Pick center of 8deg latitude bands
213  y = floor(8 * (iband - real(9.5)) * deg + real(0.5)) * tile_
214  + (northp ? 0 : utmNshift_);
215  } else {
216  // Pick point at lat 86N or 86S
217  x = ((iband & 1 ? 1 : -1) * floor(4 * deg + real(0.5))
218  + upseasting_) * tile_;
219  // Pick point at lon 90E or 90W.
220  y = upseasting_ * tile_;
221  }
222  prec = -1;
223  return;
224  } else if (len - p < 2)
225  throw GeographicErr("Missing row letter in " + mgrs);
226  const string& col = utmp ? utmcols_[zonem1 % 3] : upscols_[iband];
227  const string& row = utmp ? utmrow_ : upsrows_[northp1];
228  int icol = Utility::lookup(col, mgrs[p++]);
229  if (icol < 0)
230  throw GeographicErr("Column letter " + Utility::str(mgrs[p-1])
231  + " not in "
232  + (utmp ? "zone " + mgrs.substr(0, p-2) :
233  "UPS band " + Utility::str(mgrs[p-2]))
234  + " set " + col );
235  int irow = Utility::lookup(row, mgrs[p++]);
236  if (irow < 0)
237  throw GeographicErr("Row letter " + Utility::str(mgrs[p-1]) + " not in "
238  + (utmp ? "UTM" :
239  "UPS " + Utility::str(hemispheres_[northp1]))
240  + " set " + row);
241  if (utmp) {
242  if (zonem1 & 1)
243  irow = (irow + utmrowperiod_ - utmevenrowshift_) % utmrowperiod_;
244  iband -= 10;
245  irow = UTMRow(iband, icol, irow);
246  if (irow == maxutmSrow_)
247  throw GeographicErr("Block " + mgrs.substr(p-2, 2)
248  + " not in zone/band " + mgrs.substr(0, p-2));
249 
250  irow = northp1 ? irow : irow + 100;
251  icol = icol + minutmcol_;
252  } else {
253  bool eastp = iband & 1;
254  icol += eastp ? upseasting_ : (northp1 ? minupsNind_ : minupsSind_);
255  irow += northp1 ? minupsNind_ : minupsSind_;
256  }
257  int prec1 = (len - p)/2;
258  real
259  unit = tile_,
260  x1 = unit * icol,
261  y1 = unit * irow;
262  for (int i = 0; i < prec1; ++i) {
263  unit /= base_;
264  int
265  ix = Utility::lookup(digits_, mgrs[p + i]),
266  iy = Utility::lookup(digits_, mgrs[p + i + prec1]);
267  if (ix < 0 || iy < 0)
268  throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p));
269  x1 += unit * ix;
270  y1 += unit * iy;
271  }
272  if ((len - p) % 2) {
273  if (Utility::lookup(digits_, mgrs[len - 1]) < 0)
274  throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p));
275  else
276  throw GeographicErr("Not an even number of digits_ in "
277  + mgrs.substr(p));
278  }
279  if (prec1 > maxprec_)
280  throw GeographicErr("More than " + Utility::str(2*maxprec_)
281  + " digits_ in "
282  + mgrs.substr(p));
283  if (centerp) {
284  x1 += unit/2;
285  y1 += unit/2;
286  }
287  zone = zone1;
288  northp = northp1;
289  x = x1;
290  y = y1;
291  prec = prec1;
292  }
293 
294  void MGRS::CheckCoords(bool utmp, bool& northp, real& x, real& y) {
295  // Limits are all multiples of 100km and are all closed on the lower end
296  // and open on the upper end -- and this is reflected in the error
297  // messages. However if a coordinate lies on the excluded upper end (e.g.,
298  // after rounding), it is shifted down by eps_. This also folds UTM
299  // northings to the correct N/S hemisphere.
300  int
301  ix = int(floor(x / tile_)),
302  iy = int(floor(y / tile_)),
303  ind = (utmp ? 2 : 0) + (northp ? 1 : 0);
304  if (! (ix >= mineasting_[ind] && ix < maxeasting_[ind]) ) {
305  if (ix == maxeasting_[ind] && x == maxeasting_[ind] * tile_)
306  x -= eps_;
307  else
308  throw GeographicErr("Easting " + Utility::str(int(floor(x/1000)))
309  + "km not in MGRS/"
310  + (utmp ? "UTM" : "UPS") + " range for "
311  + (northp ? "N" : "S" ) + " hemisphere ["
312  + Utility::str(mineasting_[ind]*tile_/1000)
313  + "km, "
314  + Utility::str(maxeasting_[ind]*tile_/1000)
315  + "km)");
316  }
317  if (! (iy >= minnorthing_[ind] && iy < maxnorthing_[ind]) ) {
318  if (iy == maxnorthing_[ind] && y == maxnorthing_[ind] * tile_)
319  y -= eps_;
320  else
321  throw GeographicErr("Northing " + Utility::str(int(floor(y/1000)))
322  + "km not in MGRS/"
323  + (utmp ? "UTM" : "UPS") + " range for "
324  + (northp ? "N" : "S" ) + " hemisphere ["
325  + Utility::str(minnorthing_[ind]*tile_/1000)
326  + "km, "
327  + Utility::str(maxnorthing_[ind]*tile_/1000)
328  + "km)");
329  }
330 
331  // Correct the UTM northing and hemisphere if necessary
332  if (utmp) {
333  if (northp && iy < minutmNrow_) {
334  northp = false;
335  y += utmNshift_;
336  } else if (!northp && iy >= maxutmSrow_) {
337  if (y == maxutmSrow_ * tile_)
338  // If on equator retain S hemisphere
339  y -= eps_;
340  else {
341  northp = true;
342  y -= utmNshift_;
343  }
344  }
345  }
346  }
347 
348  int MGRS::UTMRow(int iband, int icol, int irow) {
349  // Input is MGRS (periodic) row index and output is true row index. Band
350  // index is in [-10, 10) (as returned by LatitudeBand). Column index
351  // origin is easting = 100km. Returns maxutmSrow_ if irow and iband are
352  // incompatible. Row index origin is equator.
353 
354  // Estimate center row number for latitude band
355  // 90 deg = 100 tiles; 1 band = 8 deg = 100*8/90 tiles
356  real c = 100 * (8 * iband + 4)/real(90);
357  bool northp = iband >= 0;
358  // iband minrow maxrow
359  // -10 -90 -81
360  // -9 -80 -72
361  // -8 -71 -63
362  // -7 -63 -54
363  // -6 -54 -45
364  // -5 -45 -36
365  // -4 -36 -27
366  // -3 -27 -18
367  // -2 -18 -9
368  // -1 -9 -1
369  // 0 0 8
370  // 1 8 17
371  // 2 17 26
372  // 3 26 35
373  // 4 35 44
374  // 5 44 53
375  // 6 53 62
376  // 7 62 70
377  // 8 71 79
378  // 9 80 94
379  int
380  minrow = iband > -10 ?
381  int(floor(c - real(4.3) - real(0.1) * northp)) : -90,
382  maxrow = iband < 9 ?
383  int(floor(c + real(4.4) - real(0.1) * northp)) : 94,
384  baserow = (minrow + maxrow) / 2 - utmrowperiod_ / 2;
385  // Add maxutmSrow_ = 5 * utmrowperiod_ to ensure operand is positive
386  irow = (irow - baserow + maxutmSrow_) % utmrowperiod_ + baserow;
387  if (irow < minrow || irow > maxrow) {
388  // Northing = 71*100km and 80*100km intersect band boundaries
389  // The following deals with these special cases.
390  int
391  // Fold [-10,-1] -> [9,0]
392  sband = iband >= 0 ? iband : -iband - 1,
393  // Fold [-90,-1] -> [89,0]
394  srow = irow >= 0 ? irow : -irow - 1,
395  // Fold [4,7] -> [3,0]
396  scol = icol < 4 ? icol : -icol + 7;
397  if ( ! ( (srow == 70 && sband == 8 && scol >= 2) ||
398  (srow == 71 && sband == 7 && scol <= 2) ||
399  (srow == 79 && sband == 9 && scol >= 1) ||
400  (srow == 80 && sband == 8 && scol <= 1) ) )
401  irow = maxutmSrow_;
402  }
403  return irow;
404  }
405 
406 } // namespace GeographicLib