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