GeographicLib  1.21
MGRS.cpp
Go to the documentation of this file.
00001 /**
00002  * \file MGRS.cpp
00003  * \brief Implementation for GeographicLib::MGRS class
00004  *
00005  * Copyright (c) Charles Karney (2008-2012) <charles@karney.com> and licensed
00006  * under the MIT/X11 License.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  **********************************************************************/
00009 
00010 #include <GeographicLib/MGRS.hpp>
00011 #include <GeographicLib/Utility.hpp>
00012 
00013 #define GEOGRAPHICLIB_MGRS_CPP "$Id: e4e6b419c8cd8544b3edab85b3352add0d1dd7cb $"
00014 
00015 RCSID_DECL(GEOGRAPHICLIB_MGRS_CPP)
00016 RCSID_DECL(GEOGRAPHICLIB_MGRS_HPP)
00017 
00018 namespace GeographicLib {
00019 
00020   using namespace std;
00021 
00022   const Math::real MGRS::eps_ =
00023     // 25 = ceil(log_2(2e7)) -- use half circumference here because northing
00024     // 195e5 is a legal in the "southern" hemisphere.
00025     pow(real(0.5), numeric_limits<real>::digits - 25);
00026   const Math::real MGRS::angeps_ =
00027     // 7 = ceil(log_2(90))
00028     pow(real(0.5), numeric_limits<real>::digits - 7);
00029   const string MGRS::hemispheres_ = "SN";
00030   const string MGRS::utmcols_[3] = { "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     if (zone == UTMUPS::INVALID ||
00054         Math::isnan(x) || Math::isnan(y) || Math::isnan(lat)) {
00055       prec = -1;
00056       mgrs = "INVALID";
00057       return;
00058     }
00059     bool utmp = zone != 0;
00060     CheckCoords(utmp, northp, x, y);
00061     if (!(zone >= UTMUPS::MINZONE && zone <= UTMUPS::MAXZONE))
00062       throw GeographicErr("Zone " + Utility::str(zone) + " not in [0,60]");
00063     if (!(prec >= 0 && prec <= maxprec_))
00064       throw GeographicErr("MGRS precision " + Utility::str(prec)
00065                           + " not in [0, "
00066                           + Utility::str(int(maxprec_)) + "]");
00067     // Fixed char array for accumulating string.  Allow space for zone, 3 block
00068     // letters, easting + northing.  Don't need to allow for terminating null.
00069     char mgrs1[2 + 3 + 2 * maxprec_];
00070     int
00071       zone1 = zone - 1,
00072       z = utmp ? 2 : 0,
00073       mlen = z + 3 + 2 * prec;
00074     if (utmp) {
00075       mgrs1[0] = digits_[ zone / base_ ];
00076       mgrs1[1] = digits_[ zone % base_ ];
00077       // This isn't necessary...!  Keep y non-neg
00078       // if (!northp) y -= maxutmSrow_ * tile_;
00079     }
00080     int
00081       xh = int(floor(x)) / tile_,
00082       yh = int(floor(y)) / tile_;
00083     real
00084       xf = x - tile_ * xh,
00085       yf = y - tile_ * yh;
00086     if (utmp) {
00087       int
00088         // Correct fuzziness in latitude near equator
00089         iband = abs(lat) > angeps_ ? LatitudeBand(lat) : (northp ? 0 : -1),
00090         icol = xh - minutmcol_,
00091         irow = UTMRow(iband, icol, yh % utmrowperiod_);
00092       if (irow != yh - (northp ? minutmNrow_ : maxutmSrow_))
00093         throw GeographicErr("Latitude " + Utility::str(lat)
00094                             + " is inconsistent with UTM coordinates");
00095       mgrs1[z++] = latband_[10 + iband];
00096       mgrs1[z++] = utmcols_[zone1 % 3][icol];
00097       mgrs1[z++] = utmrow_[(yh + (zone1 & 1 ? utmevenrowshift_ : 0))
00098                          % utmrowperiod_];
00099     } else {
00100       bool eastp = xh >= upseasting_;
00101       int iband = (northp ? 2 : 0) + (eastp ? 1 : 0);
00102       mgrs1[z++] = upsband_[iband];
00103       mgrs1[z++] = upscols_[iband][xh - (eastp ? upseasting_ :
00104                                          (northp ? minupsNind_ : minupsSind_))];
00105       mgrs1[z++] = upsrows_[northp][yh - (northp ? minupsNind_ : minupsSind_)];
00106     }
00107     real mult = pow(real(base_), max(tilelevel_ - prec, 0));
00108     int
00109       ix = int(floor(xf / mult)),
00110       iy = int(floor(yf / mult));
00111     for (int c = min(prec, int(tilelevel_)); c--;) {
00112       mgrs1[z + c] = digits_[ ix % base_ ];
00113       ix /= base_;
00114       mgrs1[z + c + prec] = digits_[ iy % base_ ];
00115       iy /= base_;
00116     }
00117     if (prec > tilelevel_) {
00118       xf -= floor(xf / mult);
00119       yf -= floor(yf / mult);
00120       mult = pow(real(base_), prec - tilelevel_);
00121       ix = int(floor(xf * mult));
00122       iy = int(floor(yf * mult));
00123       for (int c = prec - tilelevel_; c--;) {
00124         mgrs1[z + c + tilelevel_] = digits_[ ix % base_ ];
00125         ix /= base_;
00126         mgrs1[z + c + tilelevel_ + prec] = digits_[ iy % base_ ];
00127         iy /= base_;
00128       }
00129     }
00130     mgrs.resize(mlen);
00131     copy(mgrs1, mgrs1 + mlen, mgrs.begin());
00132   }
00133 
00134   void MGRS::Forward(int zone, bool northp, real x, real y,
00135                      int prec, std::string& mgrs) {
00136     real lat, lon;
00137     if (zone > 0)
00138       UTMUPS::Reverse(zone, northp, x, y, lat, lon);
00139     else
00140       // Latitude isn't needed for UPS specs or for INVALID
00141       lat = 0;
00142     Forward(zone, northp, x, y, lat, prec, mgrs);
00143   }
00144 
00145   void MGRS::Reverse(const std::string& mgrs,
00146                      int& zone, bool& northp, real& x, real& y,
00147                      int& prec, bool centerp) {
00148     int
00149       p = 0,
00150       len = int(mgrs.size());
00151     if (len >= 3 &&
00152         toupper(mgrs[0]) == 'I' &&
00153         toupper(mgrs[1]) == 'N' &&
00154         toupper(mgrs[2]) == 'V') {
00155       zone = UTMUPS::INVALID;
00156       northp = false;
00157       x = y = Math::NaN<real>();
00158       prec = -1;
00159       return;
00160     }
00161     int zone1 = 0;
00162     while (p < len) {
00163       int i = Utility::lookup(digits_, mgrs[p]);
00164       if (i < 0)
00165         break;
00166       zone1 = 10 * zone1 + i;
00167       ++p;
00168     }
00169     if (p > 0 && !(zone1 >= UTMUPS::MINUTMZONE && zone1 <= UTMUPS::MAXUTMZONE))
00170       throw GeographicErr("Zone " + Utility::str(zone1) + " not in [1,60]");
00171     if (p > 2)
00172       throw GeographicErr("More than 2 digits_ at start of MGRS "
00173                           + mgrs.substr(0, p));
00174     if (len - p < 3)
00175       throw GeographicErr("MGRS string too short " + mgrs);
00176     bool utmp = zone1 != UTMUPS::UPS;
00177     int zonem1 = zone1 - 1;
00178     const string& band = utmp ? latband_ : upsband_;
00179     int iband = Utility::lookup(band, mgrs[p++]);
00180     if (iband < 0)
00181       throw GeographicErr("Band letter " + Utility::str(mgrs[p-1]) + " not in "
00182                           + (utmp ? "UTM" : "UPS") + " set " + band);
00183     bool northp1 = iband >= (utmp ? 10 : 2);
00184     const string& col = utmp ? utmcols_[zonem1 % 3] : upscols_[iband];
00185     const string& row = utmp ? utmrow_ : upsrows_[northp1];
00186     int icol = Utility::lookup(col, mgrs[p++]);
00187     if (icol < 0)
00188       throw GeographicErr("Column letter " + Utility::str(mgrs[p-1])
00189                           + " not in "
00190                           + (utmp ? "zone " + mgrs.substr(0, p-2) :
00191                              "UPS band " + Utility::str(mgrs[p-2]))
00192                           + " set " + col );
00193     int irow = Utility::lookup(row, mgrs[p++]);
00194     if (irow < 0)
00195       throw GeographicErr("Row letter " + Utility::str(mgrs[p-1]) + " not in "
00196                           + (utmp ? "UTM" :
00197                              "UPS " + Utility::str(hemispheres_[northp1]))
00198                           + " set " + row);
00199     if (utmp) {
00200       if (zonem1 & 1)
00201         irow = (irow + utmrowperiod_ - utmevenrowshift_) % utmrowperiod_;
00202       iband -= 10;
00203       irow = UTMRow(iband, icol, irow);
00204       if (irow == maxutmSrow_)
00205         throw GeographicErr("Block " + mgrs.substr(p-2, 2)
00206                             + " not in zone/band " + mgrs.substr(0, p-2));
00207 
00208       irow = northp1 ? irow : irow + 100;
00209       icol = icol + minutmcol_;
00210     } else {
00211       bool eastp = iband & 1;
00212       icol += eastp ? upseasting_ : (northp1 ? minupsNind_ : minupsSind_);
00213       irow += northp1 ? minupsNind_ : minupsSind_;
00214     }
00215     int prec1 = (len - p)/2;
00216     real
00217       unit = tile_,
00218       x1 = unit * icol,
00219       y1 = unit * irow;
00220     for (int i = 0; i < prec1; ++i) {
00221       unit /= base_;
00222       int
00223         ix = Utility::lookup(digits_, mgrs[p + i]),
00224         iy = Utility::lookup(digits_, mgrs[p + i + prec1]);
00225       if (ix < 0 || iy < 0)
00226         throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p));
00227       x1 += unit * ix;
00228       y1 += unit * iy;
00229     }
00230     if ((len - p) % 2) {
00231       if (Utility::lookup(digits_, mgrs[len - 1]) < 0)
00232         throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p));
00233       else
00234         throw GeographicErr("Not an even number of digits_ in "
00235                             + mgrs.substr(p));
00236     }
00237     if (prec1 > maxprec_)
00238       throw GeographicErr("More than " + Utility::str(2*maxprec_)
00239                           + " digits_ in "
00240                           + mgrs.substr(p));
00241     if (centerp) {
00242       x1 += unit/2;
00243       y1 += unit/2;
00244     }
00245     zone = zone1;
00246     northp = northp1;
00247     x = x1;
00248     y = y1;
00249     prec = prec1;
00250   }
00251 
00252   void MGRS::CheckCoords(bool utmp, bool& northp, real& x, real& y) {
00253     // Limits are all multiples of 100km and are all closed on the lower end
00254     // and open on the upper end -- and this is reflected in the error
00255     // messages.  However if a coordinate lies on the excluded upper end (e.g.,
00256     // after rounding), it is shifted down by eps_.  This also folds UTM
00257     // northings to the correct N/S hemisphere.
00258     int
00259       ix = int(floor(x / tile_)),
00260       iy = int(floor(y / tile_)),
00261       ind = (utmp ? 2 : 0) + (northp ? 1 : 0);
00262     if (! (ix >= mineasting_[ind] && ix < maxeasting_[ind]) ) {
00263       if (ix == maxeasting_[ind] && x == maxeasting_[ind] * tile_)
00264         x -= eps_;
00265       else
00266         throw GeographicErr("Easting " + Utility::str(int(floor(x/1000)))
00267                             + "km not in MGRS/"
00268                             + (utmp ? "UTM" : "UPS") + " range for "
00269                             + (northp ? "N" : "S" ) + " hemisphere ["
00270                             + Utility::str(mineasting_[ind]*tile_/1000)
00271                             + "km, "
00272                             + Utility::str(maxeasting_[ind]*tile_/1000)
00273                             + "km)");
00274     }
00275     if (! (iy >= minnorthing_[ind] && iy < maxnorthing_[ind]) ) {
00276       if (iy == maxnorthing_[ind] && y == maxnorthing_[ind] * tile_)
00277         y -= eps_;
00278       else
00279         throw GeographicErr("Northing " + Utility::str(int(floor(y/1000)))
00280                             + "km not in MGRS/"
00281                             + (utmp ? "UTM" : "UPS") + " range for "
00282                             + (northp ? "N" : "S" ) + " hemisphere ["
00283                             + Utility::str(minnorthing_[ind]*tile_/1000)
00284                             + "km, "
00285                             + Utility::str(maxnorthing_[ind]*tile_/1000)
00286                             + "km)");
00287     }
00288 
00289     // Correct the UTM northing and hemisphere if necessary
00290     if (utmp) {
00291       if (northp && iy < minutmNrow_) {
00292         northp = false;
00293         y += utmNshift_;
00294       } else if (!northp && iy >= maxutmSrow_) {
00295         if (y == maxutmSrow_ * tile_)
00296           // If on equator retain S hemisphere
00297           y -= eps_;
00298         else {
00299           northp = true;
00300           y -= utmNshift_;
00301         }
00302       }
00303     }
00304   }
00305 
00306   int MGRS::UTMRow(int iband, int icol, int irow) throw() {
00307     // Input is MGRS (periodic) row index and output is true row index.  Band
00308     // index is in [-10, 10) (as returned by LatitudeBand).  Column index
00309     // origin is easting = 100km.  Returns maxutmSrow_ if irow and iband are
00310     // incompatible.  Row index origin is equator.
00311 
00312     // Estimate center row number for latitude band
00313     // 90 deg = 100 tiles; 1 band = 8 deg = 100*8/90 tiles
00314     real c = 100 * (8 * iband + 4)/real(90);
00315     bool northp = iband >= 0;
00316     int
00317       minrow = iband > -10 ?
00318       int(floor(c - real(4.3) - real(0.1) * northp)) : -90,
00319       maxrow = iband <   9 ?
00320       int(floor(c + real(4.4) - real(0.1) * northp)) :  94,
00321       baserow = (minrow + maxrow) / 2 - utmrowperiod_ / 2;
00322     // Add maxutmSrow_ = 5 * utmrowperiod_ to ensure operand is positive
00323     irow = (irow - baserow + maxutmSrow_) % utmrowperiod_ + baserow;
00324     if (irow < minrow || irow > maxrow) {
00325       // Northing = 71*100km and 80*100km intersect band boundaries
00326       // The following deals with these special cases.
00327       int
00328         // Fold [-10,-1] -> [9,0]
00329         sband = iband >= 0 ? iband : -iband - 1,
00330         // Fold [-90,-1] -> [89,0]
00331         srow = irow >= 0 ? irow : -irow - 1,
00332         // Fold [4,7] -> [3,0]
00333         scol = icol < 4 ? icol : -icol + 7;
00334       if ( ! ( (srow == 70 && sband == 8 && scol >= 2) ||
00335                (srow == 71 && sband == 7 && scol <= 2) ||
00336                (srow == 79 && sband == 9 && scol >= 1) ||
00337                (srow == 80 && sband == 8 && scol <= 1) ) )
00338         irow = maxutmSrow_;
00339     }
00340     return irow;
00341   }
00342 
00343 } // namespace GeographicLib