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, 2009, 2010) <charles@karney.com>
00006  * and licensed under the LGPL.  For more information, see
00007  * http://geographiclib.sourceforge.net/
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     // 25 = ceil(log_2(2e7)) -- use half circumference here because northing
00023     // 195e5 is a legal in the "southern" hemisphere.
00024     pow(real(0.5), numeric_limits<real>::digits - 25);
00025   const Math::real MGRS::angeps =
00026     // 7 = ceil(log_2(90))
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     // Fixed char array for accumulating string.  Allow space for zone, 3 block
00066     // letters, easting + northing.  Don't need to allow for terminating null.
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       // This isn't necessary...!  Keep y non-neg
00076       // if (!northp) y -= maxutmSrow * tile;
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         // Correct fuzziness in latitude near equator
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       // Latitude isn't needed for UPS specs or for INVALID
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     // Limits are all multiples of 100km and are all closed on the lower end
00250     // and open on the upper end -- and this is reflected in the error
00251     // messages.  However if a coordinate lies on the excluded upper end (e.g.,
00252     // after rounding), it is shifted down by eps.  This also folds UTM
00253     // northings to the correct N/S hemisphere.
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     // Correct the UTM northing and hemisphere if necessary
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           // If on equator retain S hemisphere
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     // Input is MGRS (periodic) row index and output is true row index.  Band
00300     // index is in [-10, 10) (as returned by LatitudeBand).  Column index
00301     // origin is easting = 100km.  Returns maxutmSrow if irow and iband are
00302     // incompatible.  Row index origin is equator.
00303 
00304     // Estimate center row number for latitude band
00305     // 90 deg = 100 tiles; 1 band = 8 deg = 100*8/90 tiles
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     // Add maxutmSrow = 5 * utmrowperiod to ensure operand is positive
00315     irow = (irow - baserow + maxutmSrow) % utmrowperiod + baserow;
00316     if (irow < minrow || irow > maxrow) {
00317       // Northing = 71*100km and 80*100km intersect band boundaries
00318       // The following deals with these special cases.
00319       int
00320         // Fold [-10,-1] -> [9,0]
00321         sband = iband >= 0 ? iband : - iband  - 1,
00322         // Fold [-90,-1] -> [89,0]
00323         srow = irow >= 0 ? irow : -irow - 1,
00324         // Fold [4,7] -> [3,0]
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 } // namespace GeographicLib