GeographicLib  1.21
OSGB.cpp
Go to the documentation of this file.
00001 /**
00002  * \file OSGB.cpp
00003  * \brief Implementation for GeographicLib::OSGB class
00004  *
00005  * Copyright (c) Charles Karney (2010-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/OSGB.hpp>
00011 #include <GeographicLib/Utility.hpp>
00012 
00013 #define GEOGRAPHICLIB_OSGB_CPP "$Id: 4bfb35c88866ed936faad797f3cef6f4ece36196 $"
00014 
00015 RCSID_DECL(GEOGRAPHICLIB_OSGB_CPP)
00016 RCSID_DECL(GEOGRAPHICLIB_OSGB_HPP)
00017 
00018 namespace GeographicLib {
00019 
00020   using namespace std;
00021 
00022   const string OSGB::letters_ = "ABCDEFGHJKLMNOPQRSTUVWXYZ";
00023   const string OSGB::digits_ = "0123456789";
00024 
00025   const TransverseMercator
00026   OSGB::OSGBTM_(MajorRadius(), Flattening(), CentralScale());
00027 
00028   Math::real OSGB::computenorthoffset() throw() {
00029     real x, y;
00030     OSGBTM_.Forward(real(0), OriginLatitude(), real(0), x, y);
00031     return FalseNorthing() - y;
00032   }
00033 
00034   const Math::real OSGB::northoffset_ = computenorthoffset();
00035 
00036   void OSGB::GridReference(real x, real y, int prec, std::string& gridref) {
00037     CheckCoords(x, y);
00038     if (!(prec >= 0 && prec <= maxprec_))
00039       throw GeographicErr("OSGB precision " + Utility::str(prec)
00040                           + " not in [0, "
00041                           + Utility::str(int(maxprec_)) + "]");
00042     char grid[2 + 2 * maxprec_];
00043     int
00044       xh = int(floor(x)) / tile_,
00045       yh = int(floor(y)) / tile_;
00046     real
00047       xf = x - tile_ * xh,
00048       yf = y - tile_ * yh;
00049     xh += tileoffx_;
00050     yh += tileoffy_;
00051     int z = 0;
00052     grid[z++] = letters_[(tilegrid_ - (yh / tilegrid_) - 1)
00053                         * tilegrid_ + (xh / tilegrid_)];
00054     grid[z++] = letters_[(tilegrid_ - (yh % tilegrid_) - 1)
00055                         * tilegrid_ + (xh % tilegrid_)];
00056     real mult = pow(real(base_), max(tilelevel_ - prec, 0));
00057     int
00058       ix = int(floor(xf / mult)),
00059       iy = int(floor(yf / mult));
00060     for (int c = min(prec, int(tilelevel_)); c--;) {
00061       grid[z + c] = digits_[ ix % base_ ];
00062       ix /= base_;
00063       grid[z + c + prec] = digits_[ iy % base_ ];
00064       iy /= base_;
00065     }
00066     if (prec > tilelevel_) {
00067       xf -= floor(xf / mult);
00068       yf -= floor(yf / mult);
00069       mult = pow(real(base_), prec - tilelevel_);
00070       ix = int(floor(xf * mult));
00071       iy = int(floor(yf * mult));
00072       for (int c = prec - tilelevel_; c--;) {
00073         grid[z + c + tilelevel_] = digits_[ ix % base_ ];
00074         ix /= base_;
00075         grid[z + c + tilelevel_ + prec] = digits_[ iy % base_ ];
00076         iy /= base_;
00077       }
00078     }
00079     int mlen = z + 2 * prec;
00080     gridref.resize(mlen);
00081     copy(grid, grid + mlen, gridref.begin());
00082   }
00083 
00084   void OSGB::GridReference(const std::string& gridref,
00085                            real& x, real& y, int& prec,
00086                            bool centerp) {
00087     int
00088       len = int(gridref.size()),
00089       p = 0;
00090     char grid[2 + 2 * maxprec_];
00091     for (int i = 0; i < len; ++i) {
00092       if (!isspace(gridref[i])) {
00093         if (p >= 2 + 2 * maxprec_)
00094           throw GeographicErr("OSGB string " + gridref + " too long");
00095         grid[p++] = gridref[i];
00096       }
00097     }
00098     len = p;
00099     p = 0;
00100     if (len < 2)
00101       throw GeographicErr("OSGB string " + gridref + " too short");
00102     if (len % 2)
00103       throw GeographicErr("OSGB string " + gridref +
00104                           " has odd number of characters");
00105     int
00106       xh = 0,
00107       yh = 0;
00108     while (p < 2) {
00109       int i = Utility::lookup(letters_, grid[p++]);
00110       if (i < 0)
00111         throw GeographicErr("Illegal prefix character " + gridref);
00112       yh = yh * tilegrid_ + tilegrid_ - (i / tilegrid_) - 1;
00113       xh = xh * tilegrid_ + (i % tilegrid_);
00114     }
00115     xh -= tileoffx_;
00116     yh -= tileoffy_;
00117 
00118     int prec1 = (len - p)/2;
00119     real
00120       unit = tile_,
00121       x1 = unit * xh,
00122       y1 = unit * yh;
00123     for (int i = 0; i < prec1; ++i) {
00124       unit /= base_;
00125       int
00126         ix = Utility::lookup(digits_, grid[p + i]),
00127         iy = Utility::lookup(digits_, grid[p + i + prec1]);
00128       if (ix < 0 || iy < 0)
00129         throw GeographicErr("Encountered a non-digit in " + gridref);
00130       x1 += unit * ix;
00131       y1 += unit * iy;
00132     }
00133     if (centerp) {
00134       x1 += unit/2;
00135       y1 += unit/2;
00136     }
00137     x = x1;
00138     y = y1;
00139     prec = prec1;
00140   }
00141 
00142   void OSGB::CheckCoords(real x, real y) {
00143     // Limits are all multiples of 100km and are all closed on the lower end
00144     // and open on the upper end -- and this is reflected in the error
00145     // messages.
00146     if (! (x >= minx_ && x < maxx_) )
00147       throw GeographicErr("Easting " + Utility::str(int(floor(x/1000)))
00148                           + "km not in OSGB range ["
00149                           + Utility::str(minx_/1000) + "km, "
00150                           + Utility::str(maxx_/1000) + "km)");
00151     if (! (y >= miny_ && y < maxy_) )
00152       throw GeographicErr("Northing " + Utility::str(int(floor(y/1000)))
00153                           + "km not in OSGB range ["
00154                           + Utility::str(miny_/1000) + "km, "
00155                           + Utility::str(maxy_/1000) + "km)");
00156   }
00157 
00158 } // namespace GeographicLib