Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
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
00142
00143
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 }