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