GeographicLib
1.21
|
00001 /** 00002 * \file MGRS.cpp 00003 * \brief Implementation for GeographicLib::MGRS class 00004 * 00005 * Copyright (c) Charles Karney (2008-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/MGRS.hpp> 00011 #include <GeographicLib/Utility.hpp> 00012 00013 #define GEOGRAPHICLIB_MGRS_CPP "$Id: e4e6b419c8cd8544b3edab85b3352add0d1dd7cb $" 00014 00015 RCSID_DECL(GEOGRAPHICLIB_MGRS_CPP) 00016 RCSID_DECL(GEOGRAPHICLIB_MGRS_HPP) 00017 00018 namespace GeographicLib { 00019 00020 using namespace std; 00021 00022 const Math::real MGRS::eps_ = 00023 // 25 = ceil(log_2(2e7)) -- use half circumference here because northing 00024 // 195e5 is a legal in the "southern" hemisphere. 00025 pow(real(0.5), numeric_limits<real>::digits - 25); 00026 const Math::real MGRS::angeps_ = 00027 // 7 = ceil(log_2(90)) 00028 pow(real(0.5), numeric_limits<real>::digits - 7); 00029 const string MGRS::hemispheres_ = "SN"; 00030 const string MGRS::utmcols_[3] = { "ABCDEFGH", "JKLMNPQR", "STUVWXYZ" }; 00031 const string MGRS::utmrow_ = "ABCDEFGHJKLMNPQRSTUV"; 00032 const string MGRS::upscols_[4] = 00033 { "JKLPQRSTUXYZ", "ABCFGHJKLPQR", "RSTUXYZ", "ABCFGHJ" }; 00034 const string MGRS::upsrows_[2] = 00035 { "ABCDEFGHJKLMNPQRSTUVWXYZ", "ABCDEFGHJKLMNP" }; 00036 const string MGRS::latband_ = "CDEFGHJKLMNPQRSTUVWX"; 00037 const string MGRS::upsband_ = "ABYZ"; 00038 const string MGRS::digits_ = "0123456789"; 00039 00040 const int MGRS::mineasting_[4] = 00041 { minupsSind_, minupsNind_, minutmcol_, minutmcol_ }; 00042 const int MGRS::maxeasting_[4] = 00043 { maxupsSind_, maxupsNind_, maxutmcol_, maxutmcol_ }; 00044 const int MGRS::minnorthing_[4] = 00045 { minupsSind_, minupsNind_, 00046 minutmSrow_, minutmSrow_ - (maxutmSrow_ - minutmNrow_) }; 00047 const int MGRS::maxnorthing_[4] = 00048 { maxupsSind_, maxupsNind_, 00049 maxutmNrow_ + (maxutmSrow_ - minutmNrow_), maxutmNrow_ }; 00050 00051 void MGRS::Forward(int zone, bool northp, real x, real y, real lat, 00052 int prec, std::string& mgrs) { 00053 if (zone == UTMUPS::INVALID || 00054 Math::isnan(x) || Math::isnan(y) || Math::isnan(lat)) { 00055 prec = -1; 00056 mgrs = "INVALID"; 00057 return; 00058 } 00059 bool utmp = zone != 0; 00060 CheckCoords(utmp, northp, x, y); 00061 if (!(zone >= UTMUPS::MINZONE && zone <= UTMUPS::MAXZONE)) 00062 throw GeographicErr("Zone " + Utility::str(zone) + " not in [0,60]"); 00063 if (!(prec >= 0 && prec <= maxprec_)) 00064 throw GeographicErr("MGRS precision " + Utility::str(prec) 00065 + " not in [0, " 00066 + Utility::str(int(maxprec_)) + "]"); 00067 // Fixed char array for accumulating string. Allow space for zone, 3 block 00068 // letters, easting + northing. Don't need to allow for terminating null. 00069 char mgrs1[2 + 3 + 2 * maxprec_]; 00070 int 00071 zone1 = zone - 1, 00072 z = utmp ? 2 : 0, 00073 mlen = z + 3 + 2 * prec; 00074 if (utmp) { 00075 mgrs1[0] = digits_[ zone / base_ ]; 00076 mgrs1[1] = digits_[ zone % base_ ]; 00077 // This isn't necessary...! Keep y non-neg 00078 // if (!northp) y -= maxutmSrow_ * tile_; 00079 } 00080 int 00081 xh = int(floor(x)) / tile_, 00082 yh = int(floor(y)) / tile_; 00083 real 00084 xf = x - tile_ * xh, 00085 yf = y - tile_ * yh; 00086 if (utmp) { 00087 int 00088 // Correct fuzziness in latitude near equator 00089 iband = abs(lat) > angeps_ ? LatitudeBand(lat) : (northp ? 0 : -1), 00090 icol = xh - minutmcol_, 00091 irow = UTMRow(iband, icol, yh % utmrowperiod_); 00092 if (irow != yh - (northp ? minutmNrow_ : maxutmSrow_)) 00093 throw GeographicErr("Latitude " + Utility::str(lat) 00094 + " is inconsistent with UTM coordinates"); 00095 mgrs1[z++] = latband_[10 + iband]; 00096 mgrs1[z++] = utmcols_[zone1 % 3][icol]; 00097 mgrs1[z++] = utmrow_[(yh + (zone1 & 1 ? utmevenrowshift_ : 0)) 00098 % utmrowperiod_]; 00099 } else { 00100 bool eastp = xh >= upseasting_; 00101 int iband = (northp ? 2 : 0) + (eastp ? 1 : 0); 00102 mgrs1[z++] = upsband_[iband]; 00103 mgrs1[z++] = upscols_[iband][xh - (eastp ? upseasting_ : 00104 (northp ? minupsNind_ : minupsSind_))]; 00105 mgrs1[z++] = upsrows_[northp][yh - (northp ? minupsNind_ : minupsSind_)]; 00106 } 00107 real mult = pow(real(base_), max(tilelevel_ - prec, 0)); 00108 int 00109 ix = int(floor(xf / mult)), 00110 iy = int(floor(yf / mult)); 00111 for (int c = min(prec, int(tilelevel_)); c--;) { 00112 mgrs1[z + c] = digits_[ ix % base_ ]; 00113 ix /= base_; 00114 mgrs1[z + c + prec] = digits_[ iy % base_ ]; 00115 iy /= base_; 00116 } 00117 if (prec > tilelevel_) { 00118 xf -= floor(xf / mult); 00119 yf -= floor(yf / mult); 00120 mult = pow(real(base_), prec - tilelevel_); 00121 ix = int(floor(xf * mult)); 00122 iy = int(floor(yf * mult)); 00123 for (int c = prec - tilelevel_; c--;) { 00124 mgrs1[z + c + tilelevel_] = digits_[ ix % base_ ]; 00125 ix /= base_; 00126 mgrs1[z + c + tilelevel_ + prec] = digits_[ iy % base_ ]; 00127 iy /= base_; 00128 } 00129 } 00130 mgrs.resize(mlen); 00131 copy(mgrs1, mgrs1 + mlen, mgrs.begin()); 00132 } 00133 00134 void MGRS::Forward(int zone, bool northp, real x, real y, 00135 int prec, std::string& mgrs) { 00136 real lat, lon; 00137 if (zone > 0) 00138 UTMUPS::Reverse(zone, northp, x, y, lat, lon); 00139 else 00140 // Latitude isn't needed for UPS specs or for INVALID 00141 lat = 0; 00142 Forward(zone, northp, x, y, lat, prec, mgrs); 00143 } 00144 00145 void MGRS::Reverse(const std::string& mgrs, 00146 int& zone, bool& northp, real& x, real& y, 00147 int& prec, bool centerp) { 00148 int 00149 p = 0, 00150 len = int(mgrs.size()); 00151 if (len >= 3 && 00152 toupper(mgrs[0]) == 'I' && 00153 toupper(mgrs[1]) == 'N' && 00154 toupper(mgrs[2]) == 'V') { 00155 zone = UTMUPS::INVALID; 00156 northp = false; 00157 x = y = Math::NaN<real>(); 00158 prec = -1; 00159 return; 00160 } 00161 int zone1 = 0; 00162 while (p < len) { 00163 int i = Utility::lookup(digits_, mgrs[p]); 00164 if (i < 0) 00165 break; 00166 zone1 = 10 * zone1 + i; 00167 ++p; 00168 } 00169 if (p > 0 && !(zone1 >= UTMUPS::MINUTMZONE && zone1 <= UTMUPS::MAXUTMZONE)) 00170 throw GeographicErr("Zone " + Utility::str(zone1) + " not in [1,60]"); 00171 if (p > 2) 00172 throw GeographicErr("More than 2 digits_ at start of MGRS " 00173 + mgrs.substr(0, p)); 00174 if (len - p < 3) 00175 throw GeographicErr("MGRS string too short " + mgrs); 00176 bool utmp = zone1 != UTMUPS::UPS; 00177 int zonem1 = zone1 - 1; 00178 const string& band = utmp ? latband_ : upsband_; 00179 int iband = Utility::lookup(band, mgrs[p++]); 00180 if (iband < 0) 00181 throw GeographicErr("Band letter " + Utility::str(mgrs[p-1]) + " not in " 00182 + (utmp ? "UTM" : "UPS") + " set " + band); 00183 bool northp1 = iband >= (utmp ? 10 : 2); 00184 const string& col = utmp ? utmcols_[zonem1 % 3] : upscols_[iband]; 00185 const string& row = utmp ? utmrow_ : upsrows_[northp1]; 00186 int icol = Utility::lookup(col, mgrs[p++]); 00187 if (icol < 0) 00188 throw GeographicErr("Column letter " + Utility::str(mgrs[p-1]) 00189 + " not in " 00190 + (utmp ? "zone " + mgrs.substr(0, p-2) : 00191 "UPS band " + Utility::str(mgrs[p-2])) 00192 + " set " + col ); 00193 int irow = Utility::lookup(row, mgrs[p++]); 00194 if (irow < 0) 00195 throw GeographicErr("Row letter " + Utility::str(mgrs[p-1]) + " not in " 00196 + (utmp ? "UTM" : 00197 "UPS " + Utility::str(hemispheres_[northp1])) 00198 + " set " + row); 00199 if (utmp) { 00200 if (zonem1 & 1) 00201 irow = (irow + utmrowperiod_ - utmevenrowshift_) % utmrowperiod_; 00202 iband -= 10; 00203 irow = UTMRow(iband, icol, irow); 00204 if (irow == maxutmSrow_) 00205 throw GeographicErr("Block " + mgrs.substr(p-2, 2) 00206 + " not in zone/band " + mgrs.substr(0, p-2)); 00207 00208 irow = northp1 ? irow : irow + 100; 00209 icol = icol + minutmcol_; 00210 } else { 00211 bool eastp = iband & 1; 00212 icol += eastp ? upseasting_ : (northp1 ? minupsNind_ : minupsSind_); 00213 irow += northp1 ? minupsNind_ : minupsSind_; 00214 } 00215 int prec1 = (len - p)/2; 00216 real 00217 unit = tile_, 00218 x1 = unit * icol, 00219 y1 = unit * irow; 00220 for (int i = 0; i < prec1; ++i) { 00221 unit /= base_; 00222 int 00223 ix = Utility::lookup(digits_, mgrs[p + i]), 00224 iy = Utility::lookup(digits_, mgrs[p + i + prec1]); 00225 if (ix < 0 || iy < 0) 00226 throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p)); 00227 x1 += unit * ix; 00228 y1 += unit * iy; 00229 } 00230 if ((len - p) % 2) { 00231 if (Utility::lookup(digits_, mgrs[len - 1]) < 0) 00232 throw GeographicErr("Encountered a non-digit in " + mgrs.substr(p)); 00233 else 00234 throw GeographicErr("Not an even number of digits_ in " 00235 + mgrs.substr(p)); 00236 } 00237 if (prec1 > maxprec_) 00238 throw GeographicErr("More than " + Utility::str(2*maxprec_) 00239 + " digits_ in " 00240 + mgrs.substr(p)); 00241 if (centerp) { 00242 x1 += unit/2; 00243 y1 += unit/2; 00244 } 00245 zone = zone1; 00246 northp = northp1; 00247 x = x1; 00248 y = y1; 00249 prec = prec1; 00250 } 00251 00252 void MGRS::CheckCoords(bool utmp, bool& northp, real& x, real& y) { 00253 // Limits are all multiples of 100km and are all closed on the lower end 00254 // and open on the upper end -- and this is reflected in the error 00255 // messages. However if a coordinate lies on the excluded upper end (e.g., 00256 // after rounding), it is shifted down by eps_. This also folds UTM 00257 // northings to the correct N/S hemisphere. 00258 int 00259 ix = int(floor(x / tile_)), 00260 iy = int(floor(y / tile_)), 00261 ind = (utmp ? 2 : 0) + (northp ? 1 : 0); 00262 if (! (ix >= mineasting_[ind] && ix < maxeasting_[ind]) ) { 00263 if (ix == maxeasting_[ind] && x == maxeasting_[ind] * tile_) 00264 x -= eps_; 00265 else 00266 throw GeographicErr("Easting " + Utility::str(int(floor(x/1000))) 00267 + "km not in MGRS/" 00268 + (utmp ? "UTM" : "UPS") + " range for " 00269 + (northp ? "N" : "S" ) + " hemisphere [" 00270 + Utility::str(mineasting_[ind]*tile_/1000) 00271 + "km, " 00272 + Utility::str(maxeasting_[ind]*tile_/1000) 00273 + "km)"); 00274 } 00275 if (! (iy >= minnorthing_[ind] && iy < maxnorthing_[ind]) ) { 00276 if (iy == maxnorthing_[ind] && y == maxnorthing_[ind] * tile_) 00277 y -= eps_; 00278 else 00279 throw GeographicErr("Northing " + Utility::str(int(floor(y/1000))) 00280 + "km not in MGRS/" 00281 + (utmp ? "UTM" : "UPS") + " range for " 00282 + (northp ? "N" : "S" ) + " hemisphere [" 00283 + Utility::str(minnorthing_[ind]*tile_/1000) 00284 + "km, " 00285 + Utility::str(maxnorthing_[ind]*tile_/1000) 00286 + "km)"); 00287 } 00288 00289 // Correct the UTM northing and hemisphere if necessary 00290 if (utmp) { 00291 if (northp && iy < minutmNrow_) { 00292 northp = false; 00293 y += utmNshift_; 00294 } else if (!northp && iy >= maxutmSrow_) { 00295 if (y == maxutmSrow_ * tile_) 00296 // If on equator retain S hemisphere 00297 y -= eps_; 00298 else { 00299 northp = true; 00300 y -= utmNshift_; 00301 } 00302 } 00303 } 00304 } 00305 00306 int MGRS::UTMRow(int iband, int icol, int irow) throw() { 00307 // Input is MGRS (periodic) row index and output is true row index. Band 00308 // index is in [-10, 10) (as returned by LatitudeBand). Column index 00309 // origin is easting = 100km. Returns maxutmSrow_ if irow and iband are 00310 // incompatible. Row index origin is equator. 00311 00312 // Estimate center row number for latitude band 00313 // 90 deg = 100 tiles; 1 band = 8 deg = 100*8/90 tiles 00314 real c = 100 * (8 * iband + 4)/real(90); 00315 bool northp = iband >= 0; 00316 int 00317 minrow = iband > -10 ? 00318 int(floor(c - real(4.3) - real(0.1) * northp)) : -90, 00319 maxrow = iband < 9 ? 00320 int(floor(c + real(4.4) - real(0.1) * northp)) : 94, 00321 baserow = (minrow + maxrow) / 2 - utmrowperiod_ / 2; 00322 // Add maxutmSrow_ = 5 * utmrowperiod_ to ensure operand is positive 00323 irow = (irow - baserow + maxutmSrow_) % utmrowperiod_ + baserow; 00324 if (irow < minrow || irow > maxrow) { 00325 // Northing = 71*100km and 80*100km intersect band boundaries 00326 // The following deals with these special cases. 00327 int 00328 // Fold [-10,-1] -> [9,0] 00329 sband = iband >= 0 ? iband : -iband - 1, 00330 // Fold [-90,-1] -> [89,0] 00331 srow = irow >= 0 ? irow : -irow - 1, 00332 // Fold [4,7] -> [3,0] 00333 scol = icol < 4 ? icol : -icol + 7; 00334 if ( ! ( (srow == 70 && sband == 8 && scol >= 2) || 00335 (srow == 71 && sband == 7 && scol <= 2) || 00336 (srow == 79 && sband == 9 && scol >= 1) || 00337 (srow == 80 && sband == 8 && scol <= 1) ) ) 00338 irow = maxutmSrow_; 00339 } 00340 return irow; 00341 } 00342 00343 } // namespace GeographicLib