GeographicLib 2.1.1
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-2022) <charles@karney.com> and licensed
6 * under the MIT/X11 License. For more information, see
7 * https://geographiclib.sourceforge.io/
8 **********************************************************************/
9
12
13#if defined(_MSC_VER)
14// Squelch warnings about enum-float expressions and mixing enums
15# pragma warning (disable: 5055 5054)
16#endif
17
18namespace GeographicLib {
19
20 using namespace std;
21
22 const char* const MGRS::hemispheres_ = "SN";
23 const char* const MGRS::utmcols_[] = { "ABCDEFGH", "JKLMNPQR", "STUVWXYZ" };
24 const char* const MGRS::utmrow_ = "ABCDEFGHJKLMNPQRSTUV";
25 const char* const MGRS::upscols_[] =
26 { "JKLPQRSTUXYZ", "ABCFGHJKLPQR", "RSTUXYZ", "ABCFGHJ" };
27 const char* const MGRS::upsrows_[] =
28 { "ABCDEFGHJKLMNPQRSTUVWXYZ", "ABCDEFGHJKLMNP" };
29 const char* const MGRS::latband_ = "CDEFGHJKLMNPQRSTUVWX";
30 const char* const MGRS::upsband_ = "ABYZ";
31 const char* const MGRS::digits_ = "0123456789";
32
33 const int MGRS::mineasting_[] =
34 { minupsSind_, minupsNind_, minutmcol_, minutmcol_ };
35 const int MGRS::maxeasting_[] =
36 { maxupsSind_, maxupsNind_, maxutmcol_, maxutmcol_ };
37 const int MGRS::minnorthing_[] =
38 { minupsSind_, minupsNind_,
39 minutmSrow_, minutmSrow_ - (maxutmSrow_ - minutmNrow_) };
40 const int MGRS::maxnorthing_[] =
41 { maxupsSind_, maxupsNind_,
42 maxutmNrow_ + (maxutmSrow_ - minutmNrow_), maxutmNrow_ };
43
44 void MGRS::Forward(int zone, bool northp, real x, real y, real lat,
45 int prec, std::string& mgrs) {
46 using std::isnan; // Needed for Centos 7, ubuntu 14
47 // The smallest angle s.t., 90 - angeps() < 90 (approx 50e-12 arcsec)
48 // 7 = ceil(log_2(90))
49 static const real angeps = ldexp(real(1), -(Math::digits() - 7));
50 if (zone == UTMUPS::INVALID ||
51 isnan(x) || isnan(y) || isnan(lat)) {
52 mgrs = "INVALID";
53 return;
54 }
55 bool utmp = zone != 0;
56 CheckCoords(utmp, northp, x, y);
57 if (!(zone >= UTMUPS::MINZONE && zone <= UTMUPS::MAXZONE))
58 throw GeographicErr("Zone " + Utility::str(zone) + " not in [0,60]");
59 if (!(prec >= -1 && prec <= maxprec_))
60 throw GeographicErr("MGRS precision " + Utility::str(prec)
61 + " not in [-1, "
62 + Utility::str(int(maxprec_)) + "]");
63 // Fixed char array for accumulating string. Allow space for zone, 3 block
64 // letters, easting + northing. Don't need to allow for terminating null.
65 char mgrs1[2 + 3 + 2 * maxprec_];
66 int
67 zone1 = zone - 1,
68 z = utmp ? 2 : 0,
69 mlen = z + 3 + 2 * prec;
70 if (utmp) {
71 mgrs1[0] = digits_[ zone / base_ ];
72 mgrs1[1] = digits_[ zone % base_ ];
73 // This isn't necessary...! Keep y non-neg
74 // if (!northp) y -= maxutmSrow_ * tile_;
75 }
76 // The C++ standard mandates 64 bits for long long. But
77 // check, to make sure.
78 static_assert(numeric_limits<long long>::digits >= 44,
79 "long long not wide enough to store 10e12");
80 // Guard against floor(x * mult_) being computed incorrectly on some
81 // platforms. The problem occurs when x * mult_ is held in extended
82 // precision and floor is inlined. This causes tests GeoConvert1[678] to
83 // fail. Problem reported and diagnosed by Thorkil Naur with g++ 10.2.0
84 // under Cygwin.
85 GEOGRAPHICLIB_VOLATILE real xx = x * mult_;
86 GEOGRAPHICLIB_VOLATILE real yy = y * mult_;
87 long long
88 ix = (long long)(floor(xx)),
89 iy = (long long)(floor(yy)),
90 m = (long long)(mult_) * (long long)(tile_);
91 int xh = int(ix / m), yh = int(iy / m);
92 if (utmp) {
93 int
94 // Correct fuzziness in latitude near equator
95 iband = fabs(lat) < angeps ? (northp ? 0 : -1) : LatitudeBand(lat),
96 icol = xh - minutmcol_,
97 irow = UTMRow(iband, icol, yh % utmrowperiod_);
98 if (irow != yh - (northp ? minutmNrow_ : maxutmSrow_))
99 throw GeographicErr("Latitude " + Utility::str(lat)
100 + " is inconsistent with UTM coordinates");
101 mgrs1[z++] = latband_[10 + iband];
102 mgrs1[z++] = utmcols_[zone1 % 3][icol];
103 mgrs1[z++] = utmrow_[(yh + (zone1 & 1 ? utmevenrowshift_ : 0))
104 % utmrowperiod_];
105 } else {
106 bool eastp = xh >= upseasting_;
107 int iband = (northp ? 2 : 0) + (eastp ? 1 : 0);
108 mgrs1[z++] = upsband_[iband];
109 mgrs1[z++] = upscols_[iband][xh - (eastp ? upseasting_ :
110 (northp ? minupsNind_ :
111 minupsSind_))];
112 mgrs1[z++] = upsrows_[northp][yh - (northp ? minupsNind_ : minupsSind_)];
113 }
114 if (prec > 0) {
115 ix -= m * xh; iy -= m * yh;
116 long long d = (long long)(pow(real(base_), maxprec_ - prec));
117 ix /= d; iy /= d;
118 for (int c = prec; c--;) {
119 mgrs1[z + c ] = digits_[ix % base_]; ix /= base_;
120 mgrs1[z + c + prec] = digits_[iy % base_]; iy /= base_;
121 }
122 }
123 mgrs.resize(mlen);
124 copy(mgrs1, mgrs1 + mlen, mgrs.begin());
125 }
126
127 void MGRS::Forward(int zone, bool northp, real x, real y,
128 int prec, std::string& mgrs) {
129 real lat, lon;
130 if (zone > 0) {
131 // Does a rough estimate for latitude determine the latitude band?
132 real ys = northp ? y : y - utmNshift_;
133 // A cheap calculation of the latitude which results in an "allowed"
134 // latitude band would be
135 // lat = ApproxLatitudeBand(ys) * 8 + 4;
136 //
137 // Here we do a more careful job using the band letter corresponding to
138 // the actual latitude.
139 ys /= tile_;
140 if (fabs(ys) < 1)
141 lat = real(0.9) * ys; // accurate enough estimate near equator
142 else {
143 real
144 // The poleward bound is a fit from above of lat(x,y)
145 // for x = 500km and y = [0km, 950km]
146 latp = real(0.901) * ys + (ys > 0 ? 1 : -1) * real(0.135),
147 // The equatorward bound is a fit from below of lat(x,y)
148 // for x = 900km and y = [0km, 950km]
149 late = real(0.902) * ys * (1 - real(1.85e-6) * ys * ys);
150 if (LatitudeBand(latp) == LatitudeBand(late))
151 lat = latp;
152 else
153 // bounds straddle a band boundary so need to compute lat accurately
154 UTMUPS::Reverse(zone, northp, x, y, lat, lon);
155 }
156 } else
157 // Latitude isn't needed for UPS specs or for INVALID
158 lat = 0;
159 Forward(zone, northp, x, y, lat, prec, mgrs);
160 }
161
162 void MGRS::Reverse(const string& mgrs,
163 int& zone, bool& northp, real& x, real& y,
164 int& prec, bool centerp) {
165 int
166 p = 0,
167 len = int(mgrs.length());
168 if (len >= 3 &&
169 toupper(mgrs[0]) == 'I' &&
170 toupper(mgrs[1]) == 'N' &&
171 toupper(mgrs[2]) == 'V') {
172 zone = UTMUPS::INVALID;
173 northp = false;
174 x = y = Math::NaN();
175 prec = -2;
176 return;
177 }
178 int zone1 = 0;
179 while (p < len) {
180 int i = Utility::lookup(digits_, mgrs[p]);
181 if (i < 0)
182 break;
183 zone1 = 10 * zone1 + i;
184 ++p;
185 }
186 if (p > 0 && !(zone1 >= UTMUPS::MINUTMZONE && zone1 <= UTMUPS::MAXUTMZONE))
187 throw GeographicErr("Zone " + Utility::str(zone1) + " not in [1,60]");
188 if (p > 2)
189 throw GeographicErr("More than 2 digits at start of MGRS "
190 + mgrs.substr(0, p));
191 if (len - p < 1)
192 throw GeographicErr("MGRS string too short " + mgrs);
193 bool utmp = zone1 != UTMUPS::UPS;
194 int zonem1 = zone1 - 1;
195 const char* band = utmp ? latband_ : upsband_;
196 int iband = Utility::lookup(band, mgrs[p++]);
197 if (iband < 0)
198 throw GeographicErr("Band letter " + Utility::str(mgrs[p-1]) + " not in "
199 + (utmp ? "UTM" : "UPS") + " set " + band);
200 bool northp1 = iband >= (utmp ? 10 : 2);
201 if (p == len) { // Grid zone only (ignore centerp)
202 // Approx length of a degree of meridian arc in units of tile.
203 real deg = real(utmNshift_) / (Math::qd * tile_);
204 zone = zone1;
205 northp = northp1;
206 if (utmp) {
207 // Pick central meridian except for 31V
208 x = ((zone == 31 && iband == 17) ? 4 : 5) * tile_;
209 // Pick center of 8deg latitude bands
210 y = floor(8 * (iband - real(9.5)) * deg + real(0.5)) * tile_
211 + (northp ? 0 : utmNshift_);
212 } else {
213 // Pick point at lat 86N or 86S
214 x = ((iband & 1 ? 1 : -1) * floor(4 * deg + real(0.5))
215 + upseasting_) * tile_;
216 // Pick point at lon 90E or 90W.
217 y = upseasting_ * tile_;
218 }
219 prec = -1;
220 return;
221 } else if (len - p < 2)
222 throw GeographicErr("Missing row letter in " + mgrs);
223 const char* col = utmp ? utmcols_[zonem1 % 3] : upscols_[iband];
224 const char* row = utmp ? utmrow_ : upsrows_[northp1];
225 int icol = Utility::lookup(col, mgrs[p++]);
226 if (icol < 0)
227 throw GeographicErr("Column letter " + Utility::str(mgrs[p-1])
228 + " not in "
229 + (utmp ? "zone " + mgrs.substr(0, p-2) :
230 "UPS band " + Utility::str(mgrs[p-2]))
231 + " set " + col );
232 int irow = Utility::lookup(row, mgrs[p++]);
233 if (irow < 0)
234 throw GeographicErr("Row letter " + Utility::str(mgrs[p-1]) + " not in "
235 + (utmp ? "UTM" :
236 "UPS " + Utility::str(hemispheres_[northp1]))
237 + " set " + row);
238 if (utmp) {
239 if (zonem1 & 1)
240 irow = (irow + utmrowperiod_ - utmevenrowshift_) % utmrowperiod_;
241 iband -= 10;
242 irow = UTMRow(iband, icol, irow);
243 if (irow == maxutmSrow_)
244 throw GeographicErr("Block " + mgrs.substr(p-2, 2)
245 + " not in zone/band " + mgrs.substr(0, p-2));
246
247 irow = northp1 ? irow : irow + 100;
248 icol = icol + minutmcol_;
249 } else {
250 bool eastp = iband & 1;
251 icol += eastp ? upseasting_ : (northp1 ? minupsNind_ : minupsSind_);
252 irow += northp1 ? minupsNind_ : minupsSind_;
253 }
254 int prec1 = (len - p)/2;
255 real
256 unit = 1,
257 x1 = icol,
258 y1 = irow;
259 for (int i = 0; i < prec1; ++i) {
260 unit *= base_;
261 int
262 ix = Utility::lookup(digits_, mgrs[p + i]),
263 iy = Utility::lookup(digits_, mgrs[p + i + prec1]);
264 if (ix < 0 || iy < 0)
265 throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p));
266 x1 = base_ * x1 + ix;
267 y1 = base_ * y1 + iy;
268 }
269 if ((len - p) % 2) {
270 if (Utility::lookup(digits_, mgrs[len - 1]) < 0)
271 throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p));
272 else
273 throw GeographicErr("Not an even number of digits in "
274 + mgrs.substr(p));
275 }
276 if (prec1 > maxprec_)
277 throw GeographicErr("More than " + Utility::str(2*maxprec_)
278 + " digits in " + mgrs.substr(p));
279 if (centerp) {
280 unit *= 2; x1 = 2 * x1 + 1; y1 = 2 * y1 + 1;
281 }
282 zone = zone1;
283 northp = northp1;
284 x = (tile_ * x1) / unit;
285 y = (tile_ * y1) / unit;
286 prec = prec1;
287 }
288
289 void MGRS::CheckCoords(bool utmp, bool& northp, real& x, real& y) {
290 // Limits are all multiples of 100km and are all closed on the lower end
291 // and open on the upper end -- and this is reflected in the error
292 // messages. However if a coordinate lies on the excluded upper end (e.g.,
293 // after rounding), it is shifted down by eps. This also folds UTM
294 // northings to the correct N/S hemisphere.
295
296 // The smallest length s.t., 1.0e7 - eps() < 1.0e7 (approx 1.9 nm)
297 // 25 = ceil(log_2(2e7)) -- use half circumference here because
298 // northing 195e5 is a legal in the "southern" hemisphere.
299 static const real eps = ldexp(real(1), -(Math::digits() - 25));
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 iband = band index in [-10, 10) (as returned by LatitudeBand),
350 // icol = column index in [0,8) with origin of easting = 100km, and irow =
351 // periodic row index in [0,20) with origin = equator. Output is true row
352 // index in [-90, 95). Returns maxutmSrow_ = 100, if irow and iband are
353 // incompatible.
354
355 // Estimate center row number for latitude band
356 // 90 deg = 100 tiles; 1 band = 8 deg = 100*8/90 tiles
357 real c = 100 * (8 * iband + 4) / real(Math::qd);
358 bool northp = iband >= 0;
359 // These are safe bounds on the rows
360 // iband minrow maxrow
361 // -10 -90 -81
362 // -9 -80 -72
363 // -8 -71 -63
364 // -7 -63 -54
365 // -6 -54 -45
366 // -5 -45 -36
367 // -4 -36 -27
368 // -3 -27 -18
369 // -2 -18 -9
370 // -1 -9 -1
371 // 0 0 8
372 // 1 8 17
373 // 2 17 26
374 // 3 26 35
375 // 4 35 44
376 // 5 44 53
377 // 6 53 62
378 // 7 62 70
379 // 8 71 79
380 // 9 80 94
381 int
382 minrow = iband > -10 ?
383 int(floor(c - real(4.3) - real(0.1) * northp)) : -90,
384 maxrow = iband < 9 ?
385 int(floor(c + real(4.4) - real(0.1) * northp)) : 94,
386 baserow = (minrow + maxrow) / 2 - utmrowperiod_ / 2;
387 // Offset irow by the multiple of utmrowperiod_ which brings it as close as
388 // possible to the center of the latitude band, (minrow + maxrow) / 2.
389 // (Add maxutmSrow_ = 5 * utmrowperiod_ to ensure operand is positive.)
390 irow = (irow - baserow + maxutmSrow_) % utmrowperiod_ + baserow;
391 if (!( irow >= minrow && irow <= maxrow )) {
392 // Outside the safe bounds, so need to check...
393 // Northing = 71e5 and 80e5 intersect band boundaries
394 // y = 71e5 in scol = 2 (x = [3e5,4e5] and x = [6e5,7e5])
395 // y = 80e5 in scol = 1 (x = [2e5,3e5] and x = [7e5,8e5])
396 // This holds for all the ellipsoids given in NGA.SIG.0012_2.0.0_UTMUPS.
397 // The following deals with these special cases.
398 int
399 // Fold [-10,-1] -> [9,0]
400 sband = iband >= 0 ? iband : -iband - 1,
401 // Fold [-90,-1] -> [89,0]
402 srow = irow >= 0 ? irow : -irow - 1,
403 // Fold [4,7] -> [3,0]
404 scol = icol < 4 ? icol : -icol + 7;
405 // For example, the safe rows for band 8 are 71 - 79. However row 70 is
406 // allowed if scol = [2,3] and row 80 is allowed if scol = [0,1].
407 if ( ! ( (srow == 70 && sband == 8 && scol >= 2) ||
408 (srow == 71 && sband == 7 && scol <= 2) ||
409 (srow == 79 && sband == 9 && scol >= 1) ||
410 (srow == 80 && sband == 8 && scol <= 1) ) )
411 irow = maxutmSrow_;
412 }
413 return irow;
414 }
415
416 void MGRS::Check() {
417 real lat, lon, x, y, t = tile_; int zone; bool northp;
418 UTMUPS::Reverse(31, true , 1*t, 0*t, lat, lon);
419 if (!( lon < 0 ))
420 throw GeographicErr("MGRS::Check: equator coverage failure");
421 UTMUPS::Reverse(31, true , 1*t, 95*t, lat, lon);
422 if (!( lat > 84 ))
423 throw GeographicErr("MGRS::Check: UTM doesn't reach latitude = 84");
424 UTMUPS::Reverse(31, false, 1*t, 10*t, lat, lon);
425 if (!( lat < -80 ))
426 throw GeographicErr("MGRS::Check: UTM doesn't reach latitude = -80");
427 UTMUPS::Forward(56, 3, zone, northp, x, y, 32);
428 if (!( x > 1*t ))
429 throw GeographicErr("MGRS::Check: Norway exception creates a gap");
430 UTMUPS::Forward(72, 21, zone, northp, x, y, 35);
431 if (!( x > 1*t ))
432 throw GeographicErr("MGRS::Check: Svalbard exception creates a gap");
433 UTMUPS::Reverse(0, true , 20*t, 13*t, lat, lon);
434 if (!( lat < 84 ))
435 throw
436 GeographicErr("MGRS::Check: North UPS doesn't reach latitude = 84");
437 UTMUPS::Reverse(0, false, 20*t, 8*t, lat, lon);
438 if (!( lat > -80 ))
439 throw
440 GeographicErr("MGRS::Check: South UPS doesn't reach latitude = -80");
441 // Entries are [band, x, y] either side of the band boundaries. Units for
442 // x, y are t = 100km.
443 const short tab[] = {
444 0, 5, 0, 0, 9, 0, // south edge of band 0
445 0, 5, 8, 0, 9, 8, // north edge of band 0
446 1, 5, 9, 1, 9, 9, // south edge of band 1
447 1, 5, 17, 1, 9, 17, // north edge of band 1
448 2, 5, 18, 2, 9, 18, // etc.
449 2, 5, 26, 2, 9, 26,
450 3, 5, 27, 3, 9, 27,
451 3, 5, 35, 3, 9, 35,
452 4, 5, 36, 4, 9, 36,
453 4, 5, 44, 4, 9, 44,
454 5, 5, 45, 5, 9, 45,
455 5, 5, 53, 5, 9, 53,
456 6, 5, 54, 6, 9, 54,
457 6, 5, 62, 6, 9, 62,
458 7, 5, 63, 7, 9, 63,
459 7, 5, 70, 7, 7, 70, 7, 7, 71, 7, 9, 71, // y = 71t crosses boundary
460 8, 5, 71, 8, 6, 71, 8, 6, 72, 8, 9, 72, // between bands 7 and 8.
461 8, 5, 79, 8, 8, 79, 8, 8, 80, 8, 9, 80, // y = 80t crosses boundary
462 9, 5, 80, 9, 7, 80, 9, 7, 81, 9, 9, 81, // between bands 8 and 9.
463 9, 5, 95, 9, 9, 95, // north edge of band 9
464 };
465 const int bandchecks = sizeof(tab) / (3 * sizeof(short));
466 for (int i = 0; i < bandchecks; ++i) {
467 UTMUPS::Reverse(38, true, tab[3*i+1]*t, tab[3*i+2]*t, lat, lon);
468 if (!( LatitudeBand(lat) == tab[3*i+0] ))
469 throw GeographicErr("MGRS::Check: Band error, b = " +
470 Utility::str(tab[3*i+0]) + ", x = " +
471 Utility::str(tab[3*i+1]) + "00km, y = " +
472 Utility::str(tab[3*i+2]) + "00km");
473 }
474 }
475
476} // namespace GeographicLib
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
Header for GeographicLib::MGRS class.
#define GEOGRAPHICLIB_VOLATILE
Definition: Math.hpp:58
Header for GeographicLib::Utility class.
Exception handling for GeographicLib.
Definition: Constants.hpp:316
static void Reverse(const std::string &mgrs, int &zone, bool &northp, real &x, real &y, int &prec, bool centerp=true)
Definition: MGRS.cpp:162
static void Check()
Definition: MGRS.cpp:416
static void Forward(int zone, bool northp, real x, real y, int prec, std::string &mgrs)
Definition: MGRS.cpp:127
static int digits()
Definition: Math.cpp:26
static T NaN()
Definition: Math.cpp:250
@ qd
degrees per quarter turn
Definition: Math.hpp:141
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:70
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:124
static int lookup(const std::string &s, char c)
Definition: Utility.cpp:160
static std::string str(T x, int p=-1)
Definition: Utility.hpp:161
Namespace for GeographicLib.
Definition: Accumulator.cpp:12