GeographicLib: MGRS.cpp Source File (original) (raw)

15 using namespace std;

17 const char* const MGRS::hemispheres_ = "SN";

18 const char* const MGRS::utmcols_[] = { "ABCDEFGH", "JKLMNPQR", "STUVWXYZ" };

19 const char* const MGRS::utmrow_ = "ABCDEFGHJKLMNPQRSTUV";

20 const char* const MGRS::upscols_[] =

21 { "JKLPQRSTUXYZ", "ABCFGHJKLPQR", "RSTUXYZ", "ABCFGHJ" };

22 const char* const MGRS::upsrows_[] =

23 { "ABCDEFGHJKLMNPQRSTUVWXYZ", "ABCDEFGHJKLMNP" };

24 const char* const MGRS::latband_ = "CDEFGHJKLMNPQRSTUVWX";

25 const char* const MGRS::upsband_ = "ABYZ";

26 const char* const MGRS::digits_ = "0123456789";

27 const char* const MGRS::alpha_ =

28 "ABCDEFGHJKLMNPQRSTUVWXYZabcdefghjklmnpqrstuvwxyz";

30 const int MGRS::mineasting_[] =

31 { minupsSind_, minupsNind_, minutmcol_, minutmcol_ };

32 const int MGRS::maxeasting_[] =

33 { maxupsSind_, maxupsNind_, maxutmcol_, maxutmcol_ };

34 const int MGRS::minnorthing_[] =

35 { minupsSind_, minupsNind_,

36 minutmSrow_, minutmSrow_ - (maxutmSrow_ - minutmNrow_) };

37 const int MGRS::maxnorthing_[] =

38 { maxupsSind_, maxupsNind_,

39 maxutmNrow_ + (maxutmSrow_ - minutmNrow_), maxutmNrow_ };

41 void MGRS::Forward(int zone, bool northp, real x, real y, real lat,

42 int prec, std::string& mgrs) {

43 using std::isnan;

46 static const real angeps = ldexp(real(1), -(Math::digits() - 7));

48 isnan(x) || isnan(y) || isnan(lat)) {

49 mgrs = "INVALID";

50 return;

51 }

52 bool utmp = zone != 0;

53 CheckCoords(utmp, northp, x, y);

56 if (!(prec >= -1 && prec <= maxprec_))

58 + " not in [-1, "

62 char mgrs1[2 + 3 + 2 * maxprec_];

63 int

64 zone1 = zone - 1,

65 z = utmp ? 2 : 0,

66 mlen = z + 3 + 2 * prec;

67 if (utmp) {

68 mgrs1[0] = digits_[ zone / base_ ];

69 mgrs1[1] = digits_[ zone % base_ ];

72 }

75 static_assert(numeric_limits::digits >= 44,

76 "long long not wide enough to store 10e12");

84 long long

85 ix = (long long)(floor(xx)),

86 iy = (long long)(floor(yy)),

87 m = (long long)(mult_) * (long long)(tile_);

88 int xh = int(ix / m), yh = int(iy / m);

89 if (utmp) {

90 int

92 iband = fabs(lat) < angeps ? (northp ? 0 : -1) : LatitudeBand(lat),

93 icol = xh - minutmcol_,

94 irow = UTMRow(iband, icol, yh % utmrowperiod_);

95 if (irow != yh - (northp ? minutmNrow_ : maxutmSrow_))

97 + " is inconsistent with UTM coordinates");

98 mgrs1[z++] = latband_[10 + iband];

99 mgrs1[z++] = utmcols_[zone1 % 3][icol];

100 mgrs1[z++] = utmrow_[(yh + (zone1 & 1 ? utmevenrowshift_ : 0))

101 % utmrowperiod_];

102 } else {

103 bool eastp = xh >= upseasting_;

104 int iband = (northp ? 2 : 0) + (eastp ? 1 : 0);

105 mgrs1[z++] = upsband_[iband];

106 mgrs1[z++] = upscols_[iband][xh - (eastp ? upseasting_ :

107 (northp ? minupsNind_ :

108 minupsSind_))];

109 mgrs1[z++] = upsrows_[northp][yh - (northp ? minupsNind_ : minupsSind_)];

110 }

111 if (prec > 0) {

112 ix -= m * xh; iy -= m * yh;

113 long long d = (long long)(pow(real(base_), maxprec_ - prec));

114 ix /= d; iy /= d;

115 for (int c = prec; c--;) {

116 mgrs1[z + c ] = digits_[ix % base_]; ix /= base_;

117 mgrs1[z + c + prec] = digits_[iy % base_]; iy /= base_;

118 }

119 }

120 mgrs.resize(mlen);

121 copy(mgrs1, mgrs1 + mlen, mgrs.begin());

122 }

123

125 int prec, std::string& mgrs) {

126 real lat, lon;

127 if (zone > 0) {

128

129 real ys = northp ? y : y - utmNshift_;

130

131

132

133

134

135

136 ys /= tile_;

137 if (fabs(ys) < 1)

138 lat = real(0.9) * ys;

139 else {

140 real

141

142

143 latp = real(0.901) * ys + (ys > 0 ? 1 : -1) * real(0.135),

144

145

146 late = real(0.902) * ys * (1 - real(1.85e-6) * ys * ys);

147 if (LatitudeBand(latp) == LatitudeBand(late))

148 lat = latp;

149 else

150

152 }

153 } else

154

155 lat = 0;

156 Forward(zone, northp, x, y, lat, prec, mgrs);

157 }

158

160 int& zone, bool& northp, real& x, real& y,

161 int& prec, bool centerp) {

162 int

163 p = 0,

164 len = int(mgrs.length());

165 if (len >= 3 &&

166 toupper(mgrs[0]) == 'I' &&

167 toupper(mgrs[1]) == 'N' &&

168 toupper(mgrs[2]) == 'V') {

170 northp = false;

172 prec = -2;

173 return;

174 }

175 int zone1 = 0;

176 while (p < len) {

178 if (i < 0)

179 break;

180 zone1 = 10 * zone1 + i;

181 ++p;

182 }

185 if (p > 2)

186 throw GeographicErr("More than 2 digits at start of MGRS "

187 + mgrs.substr(0, p));

188 if (len - p < 1)

189 throw GeographicErr("MGRS string too short " + mgrs);

191 int zonem1 = zone1 - 1;

192 const char* band = utmp ? latband_ : upsband_;

194 if (iband < 0)

196 + (utmp ? "UTM" : "UPS") + " set " + band);

197 bool northp1 = iband >= (utmp ? 10 : 2);

198 if (p == len) {

199

200 real deg = real(utmNshift_) / (Math::qd * tile_);

201 zone = zone1;

202 northp = northp1;

203 if (utmp) {

204

205 x = ((zone == 31 && iband == 17) ? 4 : 5) * tile_;

206

207 y = floor(8 * (iband - real(9.5)) * deg + real(0.5)) * tile_

208 + (northp ? 0 : utmNshift_);

209 } else {

210

211 x = ((iband & 1 ? 1 : -1) * floor(4 * deg + real(0.5))

212 + upseasting_) * tile_;

213

214 y = upseasting_ * tile_;

215 }

216 prec = -1;

217 return;

218 } else if (len - p < 2)

219 throw GeographicErr("Missing row letter in " + mgrs);

220 const char* col = utmp ? utmcols_[zonem1 % 3] : upscols_[iband];

221 const char* row = utmp ? utmrow_ : upsrows_[northp1];

223 if (icol < 0)

225 + " not in "

226 + (utmp ? "zone " + mgrs.substr(0, p-2) :

228 + " set " + col );

230 if (irow < 0)

232 + (utmp ? "UTM" :

233 "UPS " + Utility::str(hemispheres_[northp1]))

234 + " set " + row);

235 if (utmp) {

236 if (zonem1 & 1)

237 irow = (irow + utmrowperiod_ - utmevenrowshift_) % utmrowperiod_;

238 iband -= 10;

239 irow = UTMRow(iband, icol, irow);

240 if (irow == maxutmSrow_)

241 throw GeographicErr("Block " + mgrs.substr(p-2, 2)

242 + " not in zone/band " + mgrs.substr(0, p-2));

243

244 irow = northp1 ? irow : irow + 100;

245 icol = icol + minutmcol_;

246 } else {

247 bool eastp = iband & 1;

248 icol += eastp ? upseasting_ : (northp1 ? minupsNind_ : minupsSind_);

249 irow += northp1 ? minupsNind_ : minupsSind_;

250 }

251 int prec1 = (len - p)/2;

252 real

253 unit = 1,

254 x1 = icol,

255 y1 = irow;

256 for (int i = 0; i < prec1; ++i) {

257 unit *= base_;

258 int

261 if (ix < 0 || iy < 0)

262 throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p));

263 x1 = base_ * x1 + ix;

264 y1 = base_ * y1 + iy;

265 }

266 if ((len - p) % 2) {

268 throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p));

269 else

270 throw GeographicErr("Not an even number of digits in "

271 + mgrs.substr(p));

272 }

273 if (prec1 > maxprec_)

275 + " digits in " + mgrs.substr(p));

276 if (centerp) {

277 unit *= 2; x1 = 2 * x1 + 1; y1 = 2 * y1 + 1;

278 }

279 zone = zone1;

280 northp = northp1;

281 x = (tile_ * x1) / unit;

282 y = (tile_ * y1) / unit;

283 prec = prec1;

284 }

285

286 void MGRS::CheckCoords(bool utmp, bool& northp, real& x, real& y) {

287

288

289

290

291

292

293

294

295

297 int

298 ix = int(floor(x / tile_)),

299 iy = int(floor(y / tile_)),

300 ind = (utmp ? 2 : 0) + (northp ? 1 : 0);

301 if (! (ix >= mineasting_[ind] && ix < maxeasting_[ind]) ) {

302 if (ix == maxeasting_[ind] && x == maxeasting_[ind] * tile_)

303 x -= eps;

304 else

306 + "km not in MGRS/"

307 + (utmp ? "UTM" : "UPS") + " range for "

308 + (northp ? "N" : "S" ) + " hemisphere ["

309 + Utility::str(mineasting_[ind]*tile_/1000)

310 + "km, "

311 + Utility::str(maxeasting_[ind]*tile_/1000)

312 + "km)");

313 }

314 if (! (iy >= minnorthing_[ind] && iy < maxnorthing_[ind]) ) {

315 if (iy == maxnorthing_[ind] && y == maxnorthing_[ind] * tile_)

316 y -= eps;

317 else

318 throw GeographicErr("Northing " + Utility::str(int(floor(y/1000)))

319 + "km not in MGRS/"

320 + (utmp ? "UTM" : "UPS") + " range for "

321 + (northp ? "N" : "S" ) + " hemisphere ["

322 + Utility::str(minnorthing_[ind]*tile_/1000)

323 + "km, "

324 + Utility::str(maxnorthing_[ind]*tile_/1000)

325 + "km)");

326 }

327

328

329 if (utmp) {

330 if (northp && iy < minutmNrow_) {

331 northp = false;

332 y += utmNshift_;

333 } else if (!northp && iy >= maxutmSrow_) {

334 if (y == maxutmSrow_ * tile_)

335

336 y -= eps;

337 else {

338 northp = true;

339 y -= utmNshift_;

340 }

341 }

342 }

343 }

344

345 int MGRS::UTMRow(int iband, int icol, int irow) {

346

347

348

349

350

351

352

353

355 bool northp = iband >= 0;

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378 int

379 minrow = iband > -10 ?

380 int(floor(c - real(4.3) - real(0.1) * northp)) : -90,

381 maxrow = iband < 9 ?

382 int(floor(c + real(4.4) - real(0.1) * northp)) : 94,

383 baserow = (minrow + maxrow) / 2 - utmrowperiod_ / 2;

384

385

386

387 irow = (irow - baserow + maxutmSrow_) % utmrowperiod_ + baserow;

388 if (!( irow >= minrow && irow <= maxrow )) {

389

390

391

392

393

394

395 int

396

397 sband = iband >= 0 ? iband : -iband - 1,

398

399 srow = irow >= 0 ? irow : -irow - 1,

400

401 scol = icol < 4 ? icol : -icol + 7;

402

403

404 if ( ! ( (srow == 70 && sband == 8 && scol >= 2) ||

405 (srow == 71 && sband == 7 && scol <= 2) ||

406 (srow == 79 && sband == 9 && scol >= 1) ||

407 (srow == 80 && sband == 8 && scol <= 1) ) )

408 irow = maxutmSrow_;

409 }

410 return irow;

411 }

412

414 string& gridzone, string& block,

415 string& easting, string& northing) {

416 string::size_type n = mgrs.length();

417 if (n >= 3 &&

418 toupper(mgrs[0]) == 'I' &&

419 toupper(mgrs[1]) == 'N' &&

420 toupper(mgrs[2]) == 'V') {

421 gridzone = mgrs.substr(0, 3);

422 block = easting = northing = "";

423 return;

424 }

425 string::size_type p0 = mgrs.find_first_not_of(digits_);

426 if (p0 == string::npos)

427 throw GeographicErr("MGRS::Decode: ref does not contain alpha chars");

428 if (!(p0 <= 2))

429 throw GeographicErr("MGRS::Decode: ref does not start with 0-2 digits");

430 string::size_type p1 = mgrs.find_first_of(alpha_, p0);

431 if (p1 != p0)

432 throw GeographicErr("MGRS::Decode: ref contains non alphanumeric chars");

433 p1 = min(mgrs.find_first_not_of(alpha_, p0), n);

434 if (!(p1 == p0 + 1 || p1 == p0 + 3))

435 throw GeographicErr("MGRS::Decode: ref must contain 1 or 3 alpha chars");

436 if (p1 == p0 + 1 && p1 < n)

437 throw GeographicErr("MGRS::Decode: ref contains junk after 1 alpha char");

438 if (p1 < n && (mgrs.find_first_of(digits_, p1) != p1 ||

439 mgrs.find_first_not_of(digits_, p1) != string::npos))

440 throw GeographicErr("MGRS::Decode: ref contains junk at end");

441 if ((n - p1) & 1u)

442 throw GeographicErr("MGRS::Decode: ref must end with even no of digits");

443

444 gridzone = mgrs.substr(0, p0+1);

445 block = mgrs.substr(p0+1, p1 - (p0 + 1));

446 easting = mgrs.substr(p1, (n - p1) / 2);

447 northing = mgrs.substr(p1 + (n - p1) / 2);

448 }

449

451 real lat, lon, x, y, t = tile_; int zone; bool northp;

453 if (!( lon < 0 ))

454 throw GeographicErr("MGRS::Check: equator coverage failure");

456 if (!( lat > 84 ))

457 throw GeographicErr("MGRS::Check: UTM doesn't reach latitude = 84");

459 if (!( lat < -80 ))

460 throw GeographicErr("MGRS::Check: UTM doesn't reach latitude = -80");

462 if (!( x > 1*t ))

463 throw GeographicErr("MGRS::Check: Norway exception creates a gap");

465 if (!( x > 1*t ))

466 throw GeographicErr("MGRS::Check: Svalbard exception creates a gap");

468 if (!( lat < 84 ))

469 throw

470 GeographicErr("MGRS::Check: North UPS doesn't reach latitude = 84");

472 if (!( lat > -80 ))

473 throw

474 GeographicErr("MGRS::Check: South UPS doesn't reach latitude = -80");

475

476

477 const short tab[] = {

478 0, 5, 0, 0, 9, 0,

479 0, 5, 8, 0, 9, 8,

480 1, 5, 9, 1, 9, 9,

481 1, 5, 17, 1, 9, 17,

482 2, 5, 18, 2, 9, 18,

483 2, 5, 26, 2, 9, 26,

484 3, 5, 27, 3, 9, 27,

485 3, 5, 35, 3, 9, 35,

486 4, 5, 36, 4, 9, 36,

487 4, 5, 44, 4, 9, 44,

488 5, 5, 45, 5, 9, 45,

489 5, 5, 53, 5, 9, 53,

490 6, 5, 54, 6, 9, 54,

491 6, 5, 62, 6, 9, 62,

492 7, 5, 63, 7, 9, 63,

493 7, 5, 70, 7, 7, 70, 7, 7, 71, 7, 9, 71,

494 8, 5, 71, 8, 6, 71, 8, 6, 72, 8, 9, 72,

495 8, 5, 79, 8, 8, 79, 8, 8, 80, 8, 9, 80,

496 9, 5, 80, 9, 7, 80, 9, 7, 81, 9, 9, 81,

497 9, 5, 95, 9, 9, 95,

498 };

499 const int bandchecks = sizeof(tab) / (3 * sizeof(short));

500 for (int i = 0; i < bandchecks; ++i) {

501 UTMUPS::Reverse(38, true, tab[3*i+1]*t, tab[3*i+2]*t, lat, lon);

502 if (!( LatitudeBand(lat) == tab[3*i+0] ))

503 throw GeographicErr("MGRS::Check: Band error, b = " +

507 }

508 }

509

510}

GeographicLib::Math::real real

Header for GeographicLib::MGRS class.

#define GEOGRAPHICLIB_VOLATILE

Header for GeographicLib::Utility class.

Exception handling for GeographicLib.

static void Reverse(const std::string &mgrs, int &zone, bool &northp, real &x, real &y, int &prec, bool centerp=true)

Definition MGRS.cpp:159

static void Check()

Definition MGRS.cpp:450

static void Decode(const std::string &mgrs, std::string &gridzone, std::string &block, std::string &easting, std::string &northing)

Definition MGRS.cpp:413

static void Forward(int zone, bool northp, real x, real y, int prec, std::string &mgrs)

Definition MGRS.cpp:124

static constexpr int qd

degrees per quarter turn

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)

static void Reverse(int zone, bool northp, real x, real y, real &lat, real &lon, real &gamma, real &k, bool mgrslimits=false)

Some utility routines for GeographicLib.

static int lookup(const std::string &s, char c)

static std::string str(T x, int p=-1)

Namespace for GeographicLib.