00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <GeographicLib/MGRS.hpp>
00011 #include <GeographicLib/Utility.hpp>
00012
00013 namespace GeographicLib {
00014
00015 using namespace std;
00016
00017 const string MGRS::hemispheres_ = "SN";
00018 const string MGRS::utmcols_[3] = { "ABCDEFGH", "JKLMNPQR", "STUVWXYZ" };
00019 const string MGRS::utmrow_ = "ABCDEFGHJKLMNPQRSTUV";
00020 const string MGRS::upscols_[4] =
00021 { "JKLPQRSTUXYZ", "ABCFGHJKLPQR", "RSTUXYZ", "ABCFGHJ" };
00022 const string MGRS::upsrows_[2] =
00023 { "ABCDEFGHJKLMNPQRSTUVWXYZ", "ABCDEFGHJKLMNP" };
00024 const string MGRS::latband_ = "CDEFGHJKLMNPQRSTUVWX";
00025 const string MGRS::upsband_ = "ABYZ";
00026 const string MGRS::digits_ = "0123456789";
00027
00028 const int MGRS::mineasting_[4] =
00029 { minupsSind_, minupsNind_, minutmcol_, minutmcol_ };
00030 const int MGRS::maxeasting_[4] =
00031 { maxupsSind_, maxupsNind_, maxutmcol_, maxutmcol_ };
00032 const int MGRS::minnorthing_[4] =
00033 { minupsSind_, minupsNind_,
00034 minutmSrow_, minutmSrow_ - (maxutmSrow_ - minutmNrow_) };
00035 const int MGRS::maxnorthing_[4] =
00036 { maxupsSind_, maxupsNind_,
00037 maxutmNrow_ + (maxutmSrow_ - minutmNrow_), maxutmNrow_ };
00038
00039 void MGRS::Forward(int zone, bool northp, real x, real y, real lat,
00040 int prec, std::string& mgrs) {
00041 if (zone == UTMUPS::INVALID ||
00042 Math::isnan(x) || Math::isnan(y) || Math::isnan(lat)) {
00043 mgrs = "INVALID";
00044 return;
00045 }
00046 bool utmp = zone != 0;
00047 CheckCoords(utmp, northp, x, y);
00048 if (!(zone >= UTMUPS::MINZONE && zone <= UTMUPS::MAXZONE))
00049 throw GeographicErr("Zone " + Utility::str(zone) + " not in [0,60]");
00050 if (!(prec >= -1 && prec <= maxprec_))
00051 throw GeographicErr("MGRS precision " + Utility::str(prec)
00052 + " not in [-1, "
00053 + Utility::str(int(maxprec_)) + "]");
00054
00055
00056 char mgrs1[2 + 3 + 2 * maxprec_];
00057 int
00058 zone1 = zone - 1,
00059 z = utmp ? 2 : 0,
00060 mlen = z + 3 + 2 * prec;
00061 if (utmp) {
00062 mgrs1[0] = digits_[ zone / base_ ];
00063 mgrs1[1] = digits_[ zone % base_ ];
00064
00065
00066 }
00067 int
00068 xh = int(floor(x)) / tile_,
00069 yh = int(floor(y)) / tile_;
00070 real
00071 xf = x - tile_ * xh,
00072 yf = y - tile_ * yh;
00073 if (utmp) {
00074 int
00075
00076 iband = abs(lat) > angeps() ? LatitudeBand(lat) : (northp ? 0 : -1),
00077 icol = xh - minutmcol_,
00078 irow = UTMRow(iband, icol, yh % utmrowperiod_);
00079 if (irow != yh - (northp ? minutmNrow_ : maxutmSrow_))
00080 throw GeographicErr("Latitude " + Utility::str(lat)
00081 + " is inconsistent with UTM coordinates");
00082 mgrs1[z++] = latband_[10 + iband];
00083 mgrs1[z++] = utmcols_[zone1 % 3][icol];
00084 mgrs1[z++] = utmrow_[(yh + (zone1 & 1 ? utmevenrowshift_ : 0))
00085 % utmrowperiod_];
00086 } else {
00087 bool eastp = xh >= upseasting_;
00088 int iband = (northp ? 2 : 0) + (eastp ? 1 : 0);
00089 mgrs1[z++] = upsband_[iband];
00090 mgrs1[z++] = upscols_[iband][xh - (eastp ? upseasting_ :
00091 (northp ? minupsNind_ : minupsSind_))];
00092 mgrs1[z++] = upsrows_[northp][yh - (northp ? minupsNind_ : minupsSind_)];
00093 }
00094 if (prec > 0) {
00095 real mult = pow(real(base_), max(tilelevel_ - prec, 0));
00096 int
00097 ix = int(floor(xf / mult)),
00098 iy = int(floor(yf / mult));
00099 for (int c = min(prec, int(tilelevel_)); c--;) {
00100 mgrs1[z + c] = digits_[ ix % base_ ];
00101 ix /= base_;
00102 mgrs1[z + c + prec] = digits_[ iy % base_ ];
00103 iy /= base_;
00104 }
00105 if (prec > tilelevel_) {
00106 xf -= floor(xf / mult);
00107 yf -= floor(yf / mult);
00108 mult = pow(real(base_), prec - tilelevel_);
00109 ix = int(floor(xf * mult));
00110 iy = int(floor(yf * mult));
00111 for (int c = prec - tilelevel_; c--;) {
00112 mgrs1[z + c + tilelevel_] = digits_[ ix % base_ ];
00113 ix /= base_;
00114 mgrs1[z + c + tilelevel_ + prec] = digits_[ iy % base_ ];
00115 iy /= base_;
00116 }
00117 }
00118 }
00119 mgrs.resize(mlen);
00120 copy(mgrs1, mgrs1 + mlen, mgrs.begin());
00121 }
00122
00123 void MGRS::Forward(int zone, bool northp, real x, real y,
00124 int prec, std::string& mgrs) {
00125 real lat, lon;
00126 if (zone > 0) {
00127
00128 real ys = northp ? y : y - utmNshift_;
00129
00130
00131
00132
00133
00134
00135 ys /= tile_;
00136 if (abs(ys) < 1)
00137 lat = 0.9 * ys;
00138 else {
00139 real
00140
00141
00142 latp = real(0.901) * ys + (ys > 0 ? 1 : -1) * real(0.135),
00143
00144
00145 late = real(0.902) * ys * (1 - real(1.85e-6) * ys * ys);
00146 if (LatitudeBand(latp) == LatitudeBand(late))
00147 lat = latp;
00148 else
00149
00150 UTMUPS::Reverse(zone, northp, x, y, lat, lon);
00151 }
00152 } else
00153
00154 lat = 0;
00155 Forward(zone, northp, x, y, lat, prec, mgrs);
00156 }
00157
00158 void MGRS::Reverse(const std::string& mgrs,
00159 int& zone, bool& northp, real& x, real& y,
00160 int& prec, bool centerp) {
00161 int
00162 p = 0,
00163 len = int(mgrs.size());
00164 if (len >= 3 &&
00165 toupper(mgrs[0]) == 'I' &&
00166 toupper(mgrs[1]) == 'N' &&
00167 toupper(mgrs[2]) == 'V') {
00168 zone = UTMUPS::INVALID;
00169 northp = false;
00170 x = y = Math::NaN();
00171 prec = -2;
00172 return;
00173 }
00174 int zone1 = 0;
00175 while (p < len) {
00176 int i = Utility::lookup(digits_, mgrs[p]);
00177 if (i < 0)
00178 break;
00179 zone1 = 10 * zone1 + i;
00180 ++p;
00181 }
00182 if (p > 0 && !(zone1 >= UTMUPS::MINUTMZONE && zone1 <= UTMUPS::MAXUTMZONE))
00183 throw GeographicErr("Zone " + Utility::str(zone1) + " not in [1,60]");
00184 if (p > 2)
00185 throw GeographicErr("More than 2 digits_ at start of MGRS "
00186 + mgrs.substr(0, p));
00187 if (len - p < 1)
00188 throw GeographicErr("MGRS string too short " + mgrs);
00189 bool utmp = zone1 != UTMUPS::UPS;
00190 int zonem1 = zone1 - 1;
00191 const string& band = utmp ? latband_ : upsband_;
00192 int iband = Utility::lookup(band, mgrs[p++]);
00193 if (iband < 0)
00194 throw GeographicErr("Band letter " + Utility::str(mgrs[p-1]) + " not in "
00195 + (utmp ? "UTM" : "UPS") + " set " + band);
00196 bool northp1 = iband >= (utmp ? 10 : 2);
00197 if (p == len) {
00198
00199 real deg = real(utmNshift_) / (90 * tile_);
00200 zone = zone1;
00201 northp = northp1;
00202 if (utmp) {
00203
00204 x = ((zone == 31 && iband == 17) ? 4 : 5) * tile_;
00205
00206 y = floor(8 * (iband - real(9.5)) * deg + real(0.5)) * tile_
00207 + (northp ? 0 : utmNshift_);
00208 } else {
00209
00210 x = ((iband & 1 ? 1 : -1) * floor(4 * deg + real(0.5))
00211 + upseasting_) * tile_;
00212
00213 y = upseasting_ * tile_;
00214 }
00215 prec = -1;
00216 return;
00217 } else if (len - p < 2)
00218 throw GeographicErr("Missing row letter in " + mgrs);
00219 const string& col = utmp ? utmcols_[zonem1 % 3] : upscols_[iband];
00220 const string& row = utmp ? utmrow_ : upsrows_[northp1];
00221 int icol = Utility::lookup(col, mgrs[p++]);
00222 if (icol < 0)
00223 throw GeographicErr("Column letter " + Utility::str(mgrs[p-1])
00224 + " not in "
00225 + (utmp ? "zone " + mgrs.substr(0, p-2) :
00226 "UPS band " + Utility::str(mgrs[p-2]))
00227 + " set " + col );
00228 int irow = Utility::lookup(row, mgrs[p++]);
00229 if (irow < 0)
00230 throw GeographicErr("Row letter " + Utility::str(mgrs[p-1]) + " not in "
00231 + (utmp ? "UTM" :
00232 "UPS " + Utility::str(hemispheres_[northp1]))
00233 + " set " + row);
00234 if (utmp) {
00235 if (zonem1 & 1)
00236 irow = (irow + utmrowperiod_ - utmevenrowshift_) % utmrowperiod_;
00237 iband -= 10;
00238 irow = UTMRow(iband, icol, irow);
00239 if (irow == maxutmSrow_)
00240 throw GeographicErr("Block " + mgrs.substr(p-2, 2)
00241 + " not in zone/band " + mgrs.substr(0, p-2));
00242
00243 irow = northp1 ? irow : irow + 100;
00244 icol = icol + minutmcol_;
00245 } else {
00246 bool eastp = iband & 1;
00247 icol += eastp ? upseasting_ : (northp1 ? minupsNind_ : minupsSind_);
00248 irow += northp1 ? minupsNind_ : minupsSind_;
00249 }
00250 int prec1 = (len - p)/2;
00251 real
00252 unit = tile_,
00253 x1 = unit * icol,
00254 y1 = unit * irow;
00255 for (int i = 0; i < prec1; ++i) {
00256 unit /= base_;
00257 int
00258 ix = Utility::lookup(digits_, mgrs[p + i]),
00259 iy = Utility::lookup(digits_, mgrs[p + i + prec1]);
00260 if (ix < 0 || iy < 0)
00261 throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p));
00262 x1 += unit * ix;
00263 y1 += unit * iy;
00264 }
00265 if ((len - p) % 2) {
00266 if (Utility::lookup(digits_, mgrs[len - 1]) < 0)
00267 throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p));
00268 else
00269 throw GeographicErr("Not an even number of digits_ in "
00270 + mgrs.substr(p));
00271 }
00272 if (prec1 > maxprec_)
00273 throw GeographicErr("More than " + Utility::str(2*maxprec_)
00274 + " digits_ in "
00275 + mgrs.substr(p));
00276 if (centerp) {
00277 x1 += unit/2;
00278 y1 += unit/2;
00279 }
00280 zone = zone1;
00281 northp = northp1;
00282 x = x1;
00283 y = y1;
00284 prec = prec1;
00285 }
00286
00287 void MGRS::CheckCoords(bool utmp, bool& northp, real& x, real& y) {
00288
00289
00290
00291
00292
00293 int
00294 ix = int(floor(x / tile_)),
00295 iy = int(floor(y / tile_)),
00296 ind = (utmp ? 2 : 0) + (northp ? 1 : 0);
00297 if (! (ix >= mineasting_[ind] && ix < maxeasting_[ind]) ) {
00298 if (ix == maxeasting_[ind] && x == maxeasting_[ind] * tile_)
00299 x -= eps();
00300 else
00301 throw GeographicErr("Easting " + Utility::str(int(floor(x/1000)))
00302 + "km not in MGRS/"
00303 + (utmp ? "UTM" : "UPS") + " range for "
00304 + (northp ? "N" : "S" ) + " hemisphere ["
00305 + Utility::str(mineasting_[ind]*tile_/1000)
00306 + "km, "
00307 + Utility::str(maxeasting_[ind]*tile_/1000)
00308 + "km)");
00309 }
00310 if (! (iy >= minnorthing_[ind] && iy < maxnorthing_[ind]) ) {
00311 if (iy == maxnorthing_[ind] && y == maxnorthing_[ind] * tile_)
00312 y -= eps();
00313 else
00314 throw GeographicErr("Northing " + Utility::str(int(floor(y/1000)))
00315 + "km not in MGRS/"
00316 + (utmp ? "UTM" : "UPS") + " range for "
00317 + (northp ? "N" : "S" ) + " hemisphere ["
00318 + Utility::str(minnorthing_[ind]*tile_/1000)
00319 + "km, "
00320 + Utility::str(maxnorthing_[ind]*tile_/1000)
00321 + "km)");
00322 }
00323
00324
00325 if (utmp) {
00326 if (northp && iy < minutmNrow_) {
00327 northp = false;
00328 y += utmNshift_;
00329 } else if (!northp && iy >= maxutmSrow_) {
00330 if (y == maxutmSrow_ * tile_)
00331
00332 y -= eps();
00333 else {
00334 northp = true;
00335 y -= utmNshift_;
00336 }
00337 }
00338 }
00339 }
00340
00341 int MGRS::UTMRow(int iband, int icol, int irow) {
00342
00343
00344
00345
00346
00347
00348
00349 real c = 100 * (8 * iband + 4)/real(90);
00350 bool northp = iband >= 0;
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372 int
00373 minrow = iband > -10 ?
00374 int(floor(c - real(4.3) - real(0.1) * northp)) : -90,
00375 maxrow = iband < 9 ?
00376 int(floor(c + real(4.4) - real(0.1) * northp)) : 94,
00377 baserow = (minrow + maxrow) / 2 - utmrowperiod_ / 2;
00378
00379 irow = (irow - baserow + maxutmSrow_) % utmrowperiod_ + baserow;
00380 if (irow < minrow || irow > maxrow) {
00381
00382
00383 int
00384
00385 sband = iband >= 0 ? iband : -iband - 1,
00386
00387 srow = irow >= 0 ? irow : -irow - 1,
00388
00389 scol = icol < 4 ? icol : -icol + 7;
00390 if ( ! ( (srow == 70 && sband == 8 && scol >= 2) ||
00391 (srow == 71 && sband == 7 && scol <= 2) ||
00392 (srow == 79 && sband == 9 && scol >= 1) ||
00393 (srow == 80 && sband == 8 && scol <= 1) ) )
00394 irow = maxutmSrow_;
00395 }
00396 return irow;
00397 }
00398
00399 }