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, 2011) <charles@karney.com> and licensed
00006  * under the LGPL.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  **********************************************************************/
00009 
00010 #include "GeographicLib/OSGB.hpp"
00011 
00012 #define GEOGRAPHICLIB_OSGB_CPP "$Id: OSGB.cpp 6954 2011-02-16 14:36:56Z karney $"
00013 
00014 RCSID_DECL(GEOGRAPHICLIB_OSGB_CPP)
00015 RCSID_DECL(GEOGRAPHICLIB_OSGB_HPP)
00016 
00017 namespace GeographicLib {
00018 
00019   using namespace std;
00020 
00021   const string OSGB::letters = "ABCDEFGHJKLMNOPQRSTUVWXYZ";
00022   const string OSGB::digits = "0123456789";
00023 
00024   const TransverseMercator
00025   OSGB::OSGBTM(MajorRadius(), InverseFlattening(), CentralScale());
00026 
00027   Math::real OSGB::computenorthoffset() throw() {
00028     real x, y;
00029     OSGBTM.Forward(real(0), OriginLatitude(), real(0), x, y);
00030     return FalseNorthing() - y;
00031   }
00032 
00033   const Math::real OSGB::northoffset = computenorthoffset();
00034 
00035   void OSGB::GridReference(real x, real y, int prec, std::string& gridref) {
00036     CheckCoords(x, y);
00037     if (!(prec >= 0 || prec <= maxprec))
00038       throw GeographicErr("OSGB precision " + str(prec) + " not in [0, "
00039                           + str(int(maxprec)) + "]");
00040     char grid[2 + 2 * maxprec];
00041     int
00042       xh = int(floor(x)) / tile,
00043       yh = int(floor(y)) / tile;
00044     real
00045       xf = x - tile * xh,
00046       yf = y - tile * yh;
00047     xh += tileoffx;
00048     yh += tileoffy;
00049     int z = 0;
00050     grid[z++] = letters[(tilegrid - (yh / tilegrid) - 1)
00051                         * tilegrid + (xh / tilegrid)];
00052     grid[z++] = letters[(tilegrid - (yh % tilegrid) - 1)
00053                         * tilegrid + (xh % tilegrid)];
00054     real mult = pow(real(base), max(tilelevel - prec, 0));
00055     int
00056       ix = int(floor(xf / mult)),
00057       iy = int(floor(yf / mult));
00058     for (int c = min(prec, int(tilelevel)); c--;) {
00059       grid[z + c] = digits[ ix % base ];
00060       ix /= base;
00061       grid[z + c + prec] = digits[ iy % base ];
00062       iy /= base;
00063     }
00064     if (prec > tilelevel) {
00065       xf -= floor(xf / mult);
00066       yf -= floor(yf / mult);
00067       mult = pow(real(base), prec - tilelevel);
00068       ix = int(floor(xf * mult));
00069       iy = int(floor(yf * mult));
00070       for (int c = prec - tilelevel; c--;) {
00071         grid[z + c + tilelevel] = digits[ ix % base ];
00072         ix /= base;
00073         grid[z + c + tilelevel + prec] = digits[ iy % base ];
00074         iy /= base;
00075       }
00076     }
00077     int mlen = z + 2 * prec;
00078     gridref.resize(mlen);
00079     copy(grid, grid + mlen, gridref.begin());
00080   }
00081 
00082   void OSGB::GridReference(const std::string& gridref,
00083                            real& x, real& y, int& prec,
00084                            bool centerp) {
00085     int
00086       len = int(gridref.size()),
00087       p = 0;
00088     char grid[2 + 2 * maxprec];
00089     for (int i = 0; i < len; ++i) {
00090       if (!isspace(gridref[i])) {
00091         if (p >= 2 + 2 * maxprec)
00092           throw GeographicErr("OSGB string " + gridref + " too long");
00093         grid[p++] = gridref[i];
00094       }
00095     }
00096     len = p;
00097     p = 0;
00098     if (len < 2)
00099       throw GeographicErr("OSGB string " + gridref + " too short");
00100     if (len % 2)
00101       throw GeographicErr("OSGB string " + gridref +
00102                           " has odd number of characters");
00103     int
00104       xh = 0,
00105       yh = 0;
00106     while (p < 2) {
00107       int i = lookup(letters, grid[p++]);
00108       if (i < 0)
00109         throw GeographicErr("Illegal prefix character " + gridref);
00110       yh = yh * tilegrid + tilegrid - (i / tilegrid) - 1;
00111       xh = xh * tilegrid + (i % tilegrid);
00112     }
00113     xh -= tileoffx;
00114     yh -= tileoffy;
00115 
00116     int prec1 = (len - p)/2;
00117     real
00118       unit = tile,
00119       x1 = unit * xh,
00120       y1 = unit * yh;
00121     for (int i = 0; i < prec1; ++i) {
00122       unit /= base;
00123       int
00124         ix = lookup(digits, grid[p + i]),
00125         iy = lookup(digits, grid[p + i + prec1]);
00126       if (ix < 0 || iy < 0)
00127         throw GeographicErr("Encountered a non-digit in " + gridref);
00128       x1 += unit * ix;
00129       y1 += unit * iy;
00130     }
00131     if (centerp) {
00132       x1 += unit/2;
00133       y1 += unit/2;
00134     }
00135     x = x1;
00136     y = y1;
00137     prec = prec1;
00138   }
00139 
00140   void OSGB::CheckCoords(real x, real y) {
00141     // Limits are all multiples of 100km and are all closed on the lower end
00142     // and open on the upper end -- and this is reflected in the error
00143     // messages.
00144     if (! (x >= minx && x < maxx) )
00145       throw GeographicErr("Easting " + str(int(floor(x/1000)))
00146                           + "km not in OSGB range ["
00147                           + str(minx/1000) + "km, "
00148                           + str(maxx/1000) + "km)");
00149     if (! (y >= miny && y < maxy) )
00150       throw GeographicErr("Northing " + str(int(floor(y/1000)))
00151                           + "km not in OSGB range ["
00152                           + str(miny/1000) + "km, "
00153                           + str(maxy/1000) + "km)");
00154   }
00155 } // namespace GeographicLib