GeographicLib  1.30
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>
11 #include <GeographicLib/Utility.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  UTMUPS::Reverse(zone, northp, x, y, lat, lon);
134  else
135  // Latitude isn't needed for UPS specs or for INVALID
136  lat = 0;
137  Forward(zone, northp, x, y, lat, prec, mgrs);
138  }
139 
140  void MGRS::Reverse(const std::string& mgrs,
141  int& zone, bool& northp, real& x, real& y,
142  int& prec, bool centerp) {
143  int
144  p = 0,
145  len = int(mgrs.size());
146  if (len >= 3 &&
147  toupper(mgrs[0]) == 'I' &&
148  toupper(mgrs[1]) == 'N' &&
149  toupper(mgrs[2]) == 'V') {
150  zone = UTMUPS::INVALID;
151  northp = false;
152  x = y = Math::NaN<real>();
153  prec = -1;
154  return;
155  }
156  int zone1 = 0;
157  while (p < len) {
158  int i = Utility::lookup(digits_, mgrs[p]);
159  if (i < 0)
160  break;
161  zone1 = 10 * zone1 + i;
162  ++p;
163  }
164  if (p > 0 && !(zone1 >= UTMUPS::MINUTMZONE && zone1 <= UTMUPS::MAXUTMZONE))
165  throw GeographicErr("Zone " + Utility::str(zone1) + " not in [1,60]");
166  if (p > 2)
167  throw GeographicErr("More than 2 digits_ at start of MGRS "
168  + mgrs.substr(0, p));
169  if (len - p < 3)
170  throw GeographicErr("MGRS string too short " + mgrs);
171  bool utmp = zone1 != UTMUPS::UPS;
172  int zonem1 = zone1 - 1;
173  const string& band = utmp ? latband_ : upsband_;
174  int iband = Utility::lookup(band, mgrs[p++]);
175  if (iband < 0)
176  throw GeographicErr("Band letter " + Utility::str(mgrs[p-1]) + " not in "
177  + (utmp ? "UTM" : "UPS") + " set " + band);
178  bool northp1 = iband >= (utmp ? 10 : 2);
179  const string& col = utmp ? utmcols_[zonem1 % 3] : upscols_[iband];
180  const string& row = utmp ? utmrow_ : upsrows_[northp1];
181  int icol = Utility::lookup(col, mgrs[p++]);
182  if (icol < 0)
183  throw GeographicErr("Column letter " + Utility::str(mgrs[p-1])
184  + " not in "
185  + (utmp ? "zone " + mgrs.substr(0, p-2) :
186  "UPS band " + Utility::str(mgrs[p-2]))
187  + " set " + col );
188  int irow = Utility::lookup(row, mgrs[p++]);
189  if (irow < 0)
190  throw GeographicErr("Row letter " + Utility::str(mgrs[p-1]) + " not in "
191  + (utmp ? "UTM" :
192  "UPS " + Utility::str(hemispheres_[northp1]))
193  + " set " + row);
194  if (utmp) {
195  if (zonem1 & 1)
196  irow = (irow + utmrowperiod_ - utmevenrowshift_) % utmrowperiod_;
197  iband -= 10;
198  irow = UTMRow(iband, icol, irow);
199  if (irow == maxutmSrow_)
200  throw GeographicErr("Block " + mgrs.substr(p-2, 2)
201  + " not in zone/band " + mgrs.substr(0, p-2));
202 
203  irow = northp1 ? irow : irow + 100;
204  icol = icol + minutmcol_;
205  } else {
206  bool eastp = iband & 1;
207  icol += eastp ? upseasting_ : (northp1 ? minupsNind_ : minupsSind_);
208  irow += northp1 ? minupsNind_ : minupsSind_;
209  }
210  int prec1 = (len - p)/2;
211  real
212  unit = tile_,
213  x1 = unit * icol,
214  y1 = unit * irow;
215  for (int i = 0; i < prec1; ++i) {
216  unit /= base_;
217  int
218  ix = Utility::lookup(digits_, mgrs[p + i]),
219  iy = Utility::lookup(digits_, mgrs[p + i + prec1]);
220  if (ix < 0 || iy < 0)
221  throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p));
222  x1 += unit * ix;
223  y1 += unit * iy;
224  }
225  if ((len - p) % 2) {
226  if (Utility::lookup(digits_, mgrs[len - 1]) < 0)
227  throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p));
228  else
229  throw GeographicErr("Not an even number of digits_ in "
230  + mgrs.substr(p));
231  }
232  if (prec1 > maxprec_)
233  throw GeographicErr("More than " + Utility::str(2*maxprec_)
234  + " digits_ in "
235  + mgrs.substr(p));
236  if (centerp) {
237  x1 += unit/2;
238  y1 += unit/2;
239  }
240  zone = zone1;
241  northp = northp1;
242  x = x1;
243  y = y1;
244  prec = prec1;
245  }
246 
247  void MGRS::CheckCoords(bool utmp, bool& northp, real& x, real& y) {
248  // Limits are all multiples of 100km and are all closed on the lower end
249  // and open on the upper end -- and this is reflected in the error
250  // messages. However if a coordinate lies on the excluded upper end (e.g.,
251  // after rounding), it is shifted down by eps_. This also folds UTM
252  // northings to the correct N/S hemisphere.
253  int
254  ix = int(floor(x / tile_)),
255  iy = int(floor(y / tile_)),
256  ind = (utmp ? 2 : 0) + (northp ? 1 : 0);
257  if (! (ix >= mineasting_[ind] && ix < maxeasting_[ind]) ) {
258  if (ix == maxeasting_[ind] && x == maxeasting_[ind] * tile_)
259  x -= eps_;
260  else
261  throw GeographicErr("Easting " + Utility::str(int(floor(x/1000)))
262  + "km not in MGRS/"
263  + (utmp ? "UTM" : "UPS") + " range for "
264  + (northp ? "N" : "S" ) + " hemisphere ["
265  + Utility::str(mineasting_[ind]*tile_/1000)
266  + "km, "
267  + Utility::str(maxeasting_[ind]*tile_/1000)
268  + "km)");
269  }
270  if (! (iy >= minnorthing_[ind] && iy < maxnorthing_[ind]) ) {
271  if (iy == maxnorthing_[ind] && y == maxnorthing_[ind] * tile_)
272  y -= eps_;
273  else
274  throw GeographicErr("Northing " + Utility::str(int(floor(y/1000)))
275  + "km not in MGRS/"
276  + (utmp ? "UTM" : "UPS") + " range for "
277  + (northp ? "N" : "S" ) + " hemisphere ["
278  + Utility::str(minnorthing_[ind]*tile_/1000)
279  + "km, "
280  + Utility::str(maxnorthing_[ind]*tile_/1000)
281  + "km)");
282  }
283 
284  // Correct the UTM northing and hemisphere if necessary
285  if (utmp) {
286  if (northp && iy < minutmNrow_) {
287  northp = false;
288  y += utmNshift_;
289  } else if (!northp && iy >= maxutmSrow_) {
290  if (y == maxutmSrow_ * tile_)
291  // If on equator retain S hemisphere
292  y -= eps_;
293  else {
294  northp = true;
295  y -= utmNshift_;
296  }
297  }
298  }
299  }
300 
301  int MGRS::UTMRow(int iband, int icol, int irow) throw() {
302  // Input is MGRS (periodic) row index and output is true row index. Band
303  // index is in [-10, 10) (as returned by LatitudeBand). Column index
304  // origin is easting = 100km. Returns maxutmSrow_ if irow and iband are
305  // incompatible. Row index origin is equator.
306 
307  // Estimate center row number for latitude band
308  // 90 deg = 100 tiles; 1 band = 8 deg = 100*8/90 tiles
309  real c = 100 * (8 * iband + 4)/real(90);
310  bool northp = iband >= 0;
311  int
312  minrow = iband > -10 ?
313  int(floor(c - real(4.3) - real(0.1) * northp)) : -90,
314  maxrow = iband < 9 ?
315  int(floor(c + real(4.4) - real(0.1) * northp)) : 94,
316  baserow = (minrow + maxrow) / 2 - utmrowperiod_ / 2;
317  // Add maxutmSrow_ = 5 * utmrowperiod_ to ensure operand is positive
318  irow = (irow - baserow + maxutmSrow_) % utmrowperiod_ + baserow;
319  if (irow < minrow || irow > maxrow) {
320  // Northing = 71*100km and 80*100km intersect band boundaries
321  // The following deals with these special cases.
322  int
323  // Fold [-10,-1] -> [9,0]
324  sband = iband >= 0 ? iband : -iband - 1,
325  // Fold [-90,-1] -> [89,0]
326  srow = irow >= 0 ? irow : -irow - 1,
327  // Fold [4,7] -> [3,0]
328  scol = icol < 4 ? icol : -icol + 7;
329  if ( ! ( (srow == 70 && sband == 8 && scol >= 2) ||
330  (srow == 71 && sband == 7 && scol <= 2) ||
331  (srow == 79 && sband == 9 && scol >= 1) ||
332  (srow == 80 && sband == 8 && scol <= 1) ) )
333  irow = maxutmSrow_;
334  }
335  return irow;
336  }
337 
338 } // namespace GeographicLib
Generated by   doxygen 1.8.2