GeographicLib
1.21
|
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