00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "GeographicLib/MGRS.hpp"
00011
00012 #define GEOGRAPHICLIB_MGRS_CPP "$Id: MGRS.cpp 6918 2010-12-21 12:56:07Z karney $"
00013
00014 RCSID_DECL(GEOGRAPHICLIB_MGRS_CPP)
00015 RCSID_DECL(GEOGRAPHICLIB_MGRS_HPP)
00016
00017 namespace GeographicLib {
00018
00019 using namespace std;
00020
00021 const Math::real MGRS::eps =
00022
00023
00024 pow(real(0.5), numeric_limits<real>::digits - 25);
00025 const Math::real MGRS::angeps =
00026
00027 pow(real(0.5), numeric_limits<real>::digits - 7);
00028 const string MGRS::hemispheres = "SN";
00029 const string MGRS::utmcols[3] = { "ABCDEFGH", "JKLMNPQR", "STUVWXYZ" };
00030 const string MGRS::utmrow = "ABCDEFGHJKLMNPQRSTUV";
00031 const string MGRS::upscols[4] =
00032 { "JKLPQRSTUXYZ", "ABCFGHJKLPQR", "RSTUXYZ", "ABCFGHJ" };
00033 const string MGRS::upsrows[2] =
00034 { "ABCDEFGHJKLMNPQRSTUVWXYZ", "ABCDEFGHJKLMNP" };
00035 const string MGRS::latband = "CDEFGHJKLMNPQRSTUVWX";
00036 const string MGRS::upsband = "ABYZ";
00037 const string MGRS::digits = "0123456789";
00038
00039 const int MGRS::mineasting[4] =
00040 { minupsSind, minupsNind, minutmcol, minutmcol };
00041 const int MGRS::maxeasting[4] =
00042 { maxupsSind, maxupsNind, maxutmcol, maxutmcol };
00043 const int MGRS::minnorthing[4] =
00044 { minupsSind, minupsNind,
00045 minutmSrow, minutmSrow - (maxutmSrow - minutmNrow) };
00046 const int MGRS::maxnorthing[4] =
00047 { maxupsSind, maxupsNind,
00048 maxutmNrow + (maxutmSrow - minutmNrow), maxutmNrow };
00049
00050 void MGRS::Forward(int zone, bool northp, real x, real y, real lat,
00051 int prec, std::string& mgrs) {
00052 if (zone == UTMUPS::INVALID ||
00053 Math::isnan(x) || Math::isnan(y) || Math::isnan(lat)) {
00054 prec = -1;
00055 mgrs = "INVALID";
00056 return;
00057 }
00058 bool utmp = zone != 0;
00059 CheckCoords(utmp, northp, x, y);
00060 if (!(zone >= 0 || zone <= 60))
00061 throw GeographicErr("Zone " + str(zone) + " not in [0,60]");
00062 if (!(prec >= 0 || prec <= maxprec))
00063 throw GeographicErr("MGRS precision " + str(prec) + " not in [0, "
00064 + str(int(maxprec)) + "]");
00065
00066
00067 char mgrs1[2 + 3 + 2 * maxprec];
00068 int
00069 zone1 = zone - 1,
00070 z = utmp ? 2 : 0,
00071 mlen = z + 3 + 2 * prec;
00072 if (utmp) {
00073 mgrs1[0] = digits[ zone / base ];
00074 mgrs1[1] = digits[ zone % base ];
00075
00076
00077 }
00078 int
00079 xh = int(floor(x)) / tile,
00080 yh = int(floor(y)) / tile;
00081 real
00082 xf = x - tile * xh,
00083 yf = y - tile * yh;
00084 if (utmp) {
00085 int
00086
00087 iband = abs(lat) > angeps ? LatitudeBand(lat) : (northp ? 0 : -1),
00088 icol = xh - minutmcol,
00089 irow = UTMRow(iband, icol, yh % utmrowperiod);
00090 if (irow != yh - (northp ? minutmNrow : maxutmSrow))
00091 throw GeographicErr("Latitude " + str(lat)
00092 + " is inconsistent with UTM coordinates");
00093 mgrs1[z++] = latband[10 + iband];
00094 mgrs1[z++] = utmcols[zone1 % 3][icol];
00095 mgrs1[z++] = utmrow[(yh + (zone1 & 1 ? utmevenrowshift : 0))
00096 % utmrowperiod];
00097 } else {
00098 bool eastp = xh >= upseasting;
00099 int iband = (northp ? 2 : 0) + (eastp ? 1 : 0);
00100 mgrs1[z++] = upsband[iband];
00101 mgrs1[z++] = upscols[iband][xh - (eastp ? upseasting :
00102 northp ? minupsNind : minupsSind)];
00103 mgrs1[z++] = upsrows[northp][yh - (northp ? minupsNind : minupsSind)];
00104 }
00105 real mult = pow(real(base), max(tilelevel - prec, 0));
00106 int
00107 ix = int(floor(xf / mult)),
00108 iy = int(floor(yf / mult));
00109 for (int c = min(prec, int(tilelevel)); c--;) {
00110 mgrs1[z + c] = digits[ ix % base ];
00111 ix /= base;
00112 mgrs1[z + c + prec] = digits[ iy % base ];
00113 iy /= base;
00114 }
00115 if (prec > tilelevel) {
00116 xf -= floor(xf / mult);
00117 yf -= floor(yf / mult);
00118 mult = pow(real(base), prec - tilelevel);
00119 ix = int(floor(xf * mult));
00120 iy = int(floor(yf * mult));
00121 for (int c = prec - tilelevel; c--;) {
00122 mgrs1[z + c + tilelevel] = digits[ ix % base ];
00123 ix /= base;
00124 mgrs1[z + c + tilelevel + prec] = digits[ iy % base ];
00125 iy /= base;
00126 }
00127 }
00128 mgrs.resize(mlen);
00129 copy(mgrs1, mgrs1 + mlen, mgrs.begin());
00130 }
00131
00132 void MGRS::Forward(int zone, bool northp, real x, real y,
00133 int prec, std::string& mgrs) {
00134 real lat, lon;
00135 if (zone > 0)
00136 UTMUPS::Reverse(zone, northp, x, y, lat, lon);
00137 else
00138
00139 lat = 0;
00140 Forward(zone, northp, x, y, lat, prec, mgrs);
00141 }
00142
00143 void MGRS::Reverse(const std::string& mgrs,
00144 int& zone, bool& northp, real& x, real& y,
00145 int& prec, bool centerp) {
00146 int
00147 p = 0,
00148 len = int(mgrs.size());
00149 if (len >= 3 &&
00150 toupper(mgrs[0]) == 'I' &&
00151 toupper(mgrs[1]) == 'N' &&
00152 toupper(mgrs[2]) == 'V') {
00153 zone = UTMUPS::INVALID;
00154 northp = false;
00155 x = y = Math::NaN();
00156 prec = -1;
00157 return;
00158 }
00159 int zone1 = 0;
00160 while (p < len) {
00161 int i = lookup(digits, mgrs[p]);
00162 if (i < 0)
00163 break;
00164 zone1 = 10 * zone1 + i;
00165 ++p;
00166 }
00167 if (p > 0 && (zone1 == 0 || zone1 > 60))
00168 throw GeographicErr("Zone " + str(zone1) + " not in [1,60]");
00169 if (p > 2)
00170 throw GeographicErr("More than 2 digits at start of MGRS "
00171 + mgrs.substr(0, p));
00172 if (len - p < 3)
00173 throw GeographicErr("MGRS string " + mgrs + " too short");
00174 bool utmp = zone1 != 0;
00175 int zonem1 = zone1 - 1;
00176 const string& band = utmp ? latband : upsband;
00177 int iband = lookup(band, mgrs[p++]);
00178 if (iband < 0)
00179 throw GeographicErr("Band letter " + str(mgrs[p-1]) + " not in "
00180 + (utmp ? "UTM" : "UPS") + " set " + band);
00181 bool northp1 = iband >= (utmp ? 10 : 2);
00182 const string& col = utmp ? utmcols[zonem1 % 3] : upscols[iband];
00183 const string& row = utmp ? utmrow : upsrows[northp1];
00184 int icol = lookup(col, mgrs[p++]);
00185 if (icol < 0)
00186 throw GeographicErr("Column letter " + str(mgrs[p-1]) + " not in "
00187 + (utmp ? "zone " + mgrs.substr(0, p-2) :
00188 "UPS band " + str(mgrs[p-2]))
00189 + " set " + col );
00190 int irow = lookup(row, mgrs[p++]);
00191 if (irow < 0)
00192 throw GeographicErr("Row letter " + str(mgrs[p-1]) + " not in "
00193 + (utmp ? "UTM" :
00194 "UPS " + str(hemispheres[northp1]))
00195 + " set " + row);
00196 if (utmp) {
00197 if (zonem1 & 1)
00198 irow = (irow + utmrowperiod - utmevenrowshift) % utmrowperiod;
00199 iband -= 10;
00200 irow = UTMRow(iband, icol, irow);
00201 if (irow == maxutmSrow)
00202 throw GeographicErr("Block " + mgrs.substr(p-2, 2)
00203 + " not in zone/band " + mgrs.substr(0, p-2));
00204
00205 irow = northp1 ? irow : irow + 100;
00206 icol = icol + minutmcol;
00207 } else {
00208 bool eastp = iband & 1;
00209 icol += eastp ? upseasting : northp1 ? minupsNind : minupsSind;
00210 irow += northp1 ? minupsNind : minupsSind;
00211 }
00212 int prec1 = (len - p)/2;
00213 real
00214 unit = tile,
00215 x1 = unit * icol,
00216 y1 = unit * irow;
00217 for (int i = 0; i < prec1; ++i) {
00218 unit /= base;
00219 int
00220 ix = lookup(digits, mgrs[p + i]),
00221 iy = lookup(digits, mgrs[p + i + prec1]);
00222 if (ix < 0 || iy < 0)
00223 throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p));
00224 x1 += unit * ix;
00225 y1 += unit * iy;
00226 }
00227 if ((len - p) % 2) {
00228 if (lookup(digits, mgrs[len - 1]) < 0)
00229 throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p));
00230 else
00231 throw GeographicErr("Not an even number of digits in "
00232 + mgrs.substr(p));
00233 }
00234 if (prec1 > maxprec)
00235 throw GeographicErr("More than " + str(2*maxprec) + " digits in "
00236 + mgrs.substr(p));
00237 if (centerp) {
00238 x1 += unit/2;
00239 y1 += unit/2;
00240 }
00241 zone = zone1;
00242 northp = northp1;
00243 x = x1;
00244 y = y1;
00245 prec = prec1;
00246 }
00247
00248 void MGRS::CheckCoords(bool utmp, bool& northp, real& x, real& y) {
00249
00250
00251
00252
00253
00254 int
00255 ix = int(floor(x / tile)),
00256 iy = int(floor(y / tile)),
00257 ind = (utmp ? 2 : 0) + (northp ? 1 : 0);
00258 if (! (ix >= mineasting[ind] && ix < maxeasting[ind]) ) {
00259 if (ix == maxeasting[ind] && x == maxeasting[ind] * tile)
00260 x -= eps;
00261 else
00262 throw GeographicErr("Easting " + str(int(floor(x/1000)))
00263 + "km not in MGRS/"
00264 + (utmp ? "UTM" : "UPS") + " range for "
00265 + (northp ? "N" : "S" ) + " hemisphere ["
00266 + str(mineasting[ind]*tile/1000) + "km, "
00267 + str(maxeasting[ind]*tile/1000) + "km)");
00268 }
00269 if (! (iy >= minnorthing[ind] && iy < maxnorthing[ind]) ) {
00270 if (iy == maxnorthing[ind] && y == maxnorthing[ind] * tile)
00271 y -= eps;
00272 else
00273 throw GeographicErr("Northing " + str(int(floor(y/1000)))
00274 + "km not in MGRS/"
00275 + (utmp ? "UTM" : "UPS") + " range for "
00276 + (northp ? "N" : "S" ) + " hemisphere ["
00277 + str(minnorthing[ind]*tile/1000) + "km, "
00278 + str(maxnorthing[ind]*tile/1000) + "km)");
00279 }
00280
00281
00282 if (utmp) {
00283 if (northp && iy < minutmNrow) {
00284 northp = false;
00285 y += utmNshift;
00286 } else if (!northp && iy >= maxutmSrow) {
00287 if (y == maxutmSrow * tile)
00288
00289 y -= eps;
00290 else {
00291 northp = true;
00292 y -= utmNshift;
00293 }
00294 }
00295 }
00296 }
00297
00298 int MGRS::UTMRow(int iband, int icol, int irow) throw() {
00299
00300
00301
00302
00303
00304
00305
00306 real c = 100 * (8 * iband + 4)/real(90);
00307 bool northp = iband >= 0;
00308 int
00309 minrow = iband > -10 ?
00310 int(floor(c - real(4.3) - real(0.1) * northp)) : -90,
00311 maxrow = iband < 9 ?
00312 int(floor(c + real(4.4) - real(0.1) * northp)) : 94,
00313 baserow = (minrow + maxrow) / 2 - utmrowperiod / 2;
00314
00315 irow = (irow - baserow + maxutmSrow) % utmrowperiod + baserow;
00316 if (irow < minrow || irow > maxrow) {
00317
00318
00319 int
00320
00321 sband = iband >= 0 ? iband : - iband - 1,
00322
00323 srow = irow >= 0 ? irow : -irow - 1,
00324
00325 scol = icol < 4 ? icol : -icol + 7;
00326 if ( ! ( (srow == 70 && sband == 8 && scol >= 2) ||
00327 (srow == 71 && sband == 7 && scol <= 2) ||
00328 (srow == 79 && sband == 9 && scol >= 1) ||
00329 (srow == 80 && sband == 8 && scol <= 1) ) )
00330 irow = maxutmSrow;
00331 }
00332 return irow;
00333 }
00334
00335 }