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