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 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     // 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] =
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     // Fixed char array for accumulating string.  Allow space for zone, 3 block
00061     // letters, easting + northing.  Don't need to allow for terminating null.
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       // This isn't necessary...!  Keep y non-neg
00071       // if (!northp) y -= maxutmSrow * tile;
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         // Correct fuzziness in latitude near equator
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       // Latitude isn't needed for UPS specs.
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     // Limits are all multiples of 100km and are all closed on the lower end
00235     // and open on the upper end -- and this is reflected in the error
00236     // messages.  However if a coordinate lies on the excluded upper end (e.g.,
00237     // after rounding), it is shifted down by eps.  This also folds UTM
00238     // northings to the correct N/S hemisphere.
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     // Correct the UTM northing and hemisphere if necessary
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           // If on equator retain S hemisphere
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     // Input is MGRS (periodic) row index and output is true row index.  Band
00285     // index is in [-10, 10) (as returned by LatitudeBand).  Column index
00286     // origin is easting = 100km.  Returns maxutmSrow if irow and iband are
00287     // incompatible.  Row index origin is equator.
00288 
00289     // Estimate center row number for latitude band
00290     // 90 deg = 100 tiles; 1 band = 8 deg = 100*8/90 tiles
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     // Add maxutmSrow = 5 * utmrowperiod to ensure operand is positive
00300     irow = (irow - baserow + maxutmSrow) % utmrowperiod + baserow;
00301     if (irow < minrow || irow > maxrow) {
00302       // Northing = 71*100km and 80*100km intersect band boundaries
00303       // The following deals with these special cases.
00304       int
00305         // Fold [-10,-1] -> [9,0]
00306         sband = iband >= 0 ? iband : - iband  - 1,
00307         // Fold [-90,-1] -> [89,0]
00308         srow = irow >= 0 ? irow : -irow - 1,
00309         // Fold [4,7] -> [3,0]
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 } // namespace GeographicLib

Generated on 21 May 2010 for GeographicLib by  doxygen 1.6.1