GeographicLib
1.21
|
00001 /** 00002 * \file UTMUPS.cpp 00003 * \brief Implementation for GeographicLib::UTMUPS 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/UTMUPS.hpp> 00011 #include <iomanip> 00012 #include <GeographicLib/MGRS.hpp> 00013 #include <GeographicLib/PolarStereographic.hpp> 00014 #include <GeographicLib/TransverseMercator.hpp> 00015 #include <GeographicLib/Utility.hpp> 00016 00017 #define GEOGRAPHICLIB_UTMUPS_CPP \ 00018 "$Id: 5672b003ee47cd660377c111e3fca2b81da86323 $" 00019 00020 RCSID_DECL(GEOGRAPHICLIB_UTMUPS_CPP) 00021 RCSID_DECL(GEOGRAPHICLIB_UTMUPS_HPP) 00022 00023 namespace GeographicLib { 00024 00025 using namespace std; 00026 00027 const Math::real UTMUPS::falseeasting_[4] = 00028 { MGRS::upseasting_ * MGRS::tile_, MGRS::upseasting_ * MGRS::tile_, 00029 MGRS::utmeasting_ * MGRS::tile_, MGRS::utmeasting_ * MGRS::tile_ }; 00030 const Math::real UTMUPS::falsenorthing_[4] = 00031 { MGRS::upseasting_ * MGRS::tile_, MGRS::upseasting_ * MGRS::tile_, 00032 MGRS::maxutmSrow_ * MGRS::tile_, MGRS::minutmNrow_ * MGRS::tile_ }; 00033 const Math::real UTMUPS::mineasting_[4] = 00034 { MGRS::minupsSind_ * MGRS::tile_, MGRS::minupsNind_ * MGRS::tile_, 00035 MGRS::minutmcol_ * MGRS::tile_, MGRS::minutmcol_ * MGRS::tile_ }; 00036 const Math::real UTMUPS::maxeasting_[4] = 00037 { MGRS::maxupsSind_ * MGRS::tile_, MGRS::maxupsNind_ * MGRS::tile_, 00038 MGRS::maxutmcol_ * MGRS::tile_, MGRS::maxutmcol_ * MGRS::tile_ }; 00039 const Math::real UTMUPS::minnorthing_[4] = 00040 { MGRS::minupsSind_ * MGRS::tile_, MGRS::minupsNind_ * MGRS::tile_, 00041 MGRS::minutmSrow_ * MGRS::tile_, 00042 (MGRS::minutmNrow_ + MGRS::minutmSrow_ - MGRS::maxutmSrow_) 00043 * MGRS::tile_ }; 00044 const Math::real UTMUPS::maxnorthing_[4] = 00045 { MGRS::maxupsSind_ * MGRS::tile_, MGRS::maxupsNind_ * MGRS::tile_, 00046 (MGRS::maxutmSrow_ + MGRS::maxutmNrow_ - MGRS::minutmNrow_) * MGRS::tile_, 00047 MGRS::maxutmNrow_ * MGRS::tile_ }; 00048 00049 int UTMUPS::StandardZone(real lat, real lon, int setzone) { 00050 if (!(setzone >= MINPSEUDOZONE && setzone <= MAXZONE)) 00051 throw GeographicErr("Illegal zone requested " + Utility::str(setzone)); 00052 if (setzone >= MINZONE || setzone == INVALID) 00053 return setzone; 00054 if (Math::isnan(lat) || Math::isnan(lon)) // Check if lat or lon is a NaN 00055 return INVALID; 00056 // Assume lon is in [-180, 360]. 00057 if (setzone == UTM || (lat >= -80 && lat < 84)) { 00058 // Assume lon is in [-180, 360]. 00059 int ilon = int(floor(lon)); 00060 if (ilon >= 180) 00061 ilon -= 360; 00062 int zone = (ilon + 186)/6; 00063 int band = MGRS::LatitudeBand(lat); 00064 if (band == 7 && zone == 31 && ilon >= 3) 00065 zone = 32; 00066 else if (band == 9 && ilon >= 0 && ilon < 42) 00067 zone = 2 * ((ilon + 183)/12) + 1; 00068 return zone; 00069 } else 00070 return UPS; 00071 } 00072 00073 void UTMUPS::Forward(real lat, real lon, 00074 int& zone, bool& northp, real& x, real& y, 00075 real& gamma, real& k, 00076 int setzone, bool mgrslimits) { 00077 CheckLatLon(lat, lon); 00078 bool northp1 = lat >= 0; 00079 int zone1 = StandardZone(lat, lon, setzone); 00080 if (zone1 == INVALID) { 00081 zone = zone1; 00082 northp = northp1; 00083 x = y = gamma = k = Math::NaN<real>(); 00084 return; 00085 } 00086 real x1, y1, gamma1, k1; 00087 bool utmp = zone1 != UPS; 00088 if (utmp) { 00089 real 00090 lon0 = CentralMeridian(zone1), 00091 dlon = lon - lon0; 00092 dlon = abs(dlon - 360 * floor((dlon + 180)/360)); 00093 if (dlon > 60) 00094 // Check isn't really necessary because CheckCoords catches this case. 00095 // But this allows a more meaningful error message to be given. 00096 throw GeographicErr("Longitude " + Utility::str(lon) 00097 + "d more than 60d from center of UTM zone " 00098 + Utility::str(zone1)); 00099 TransverseMercator::UTM.Forward(lon0, lat, lon, x1, y1, gamma1, k1); 00100 } else { 00101 if (abs(lat) < 70) 00102 // Check isn't really necessary ... (see above). 00103 throw GeographicErr("Latitude " + Utility::str(lat) 00104 + "d more than 20d from " 00105 + (northp1 ? "N" : "S") + " pole"); 00106 PolarStereographic::UPS.Forward(northp1, lat, lon, x1, y1, gamma1, k1); 00107 } 00108 int ind = (utmp ? 2 : 0) + (northp1 ? 1 : 0); 00109 x1 += falseeasting_[ind]; 00110 y1 += falsenorthing_[ind]; 00111 if (! CheckCoords(zone1 != UPS, northp1, x1, y1, mgrslimits, false) ) 00112 throw GeographicErr("Latitude " + Utility::str(lat) 00113 + ", longitude " + Utility::str(lon) 00114 + " out of legal range for " 00115 + (utmp ? "UTM zone " + Utility::str(zone1) : "UPS")); 00116 zone = zone1; 00117 northp = northp1; 00118 x = x1; 00119 y = y1; 00120 gamma = gamma1; 00121 k = k1; 00122 } 00123 00124 void UTMUPS::Reverse(int zone, bool northp, real x, real y, 00125 real& lat, real& lon, real& gamma, real& k, 00126 bool mgrslimits) { 00127 if (zone == INVALID || Math::isnan(x) || Math::isnan(y)) { 00128 lat = lon = gamma = k = Math::NaN<real>(); 00129 return; 00130 } 00131 if (!(zone >= MINZONE && zone <= MAXZONE)) 00132 throw GeographicErr("Zone " + Utility::str(zone) 00133 + " not in range [0, 60]"); 00134 bool utmp = zone != UPS; 00135 CheckCoords(utmp, northp, x, y, mgrslimits); 00136 int ind = (utmp ? 2 : 0) + (northp ? 1 : 0); 00137 x -= falseeasting_[ind]; 00138 y -= falsenorthing_[ind]; 00139 if (utmp) 00140 TransverseMercator::UTM.Reverse(CentralMeridian(zone), 00141 x, y, lat, lon, gamma, k); 00142 else 00143 PolarStereographic::UPS.Reverse(northp, x, y, lat, lon, gamma, k); 00144 } 00145 00146 void UTMUPS::CheckLatLon(real lat, real lon) { 00147 if (lat < -90 || lat > 90) 00148 throw GeographicErr("Latitude " + Utility::str(lat) 00149 + "d not in [-90d, 90d]"); 00150 if (lon < -180 || lon > 360) 00151 throw GeographicErr("Latitude " + Utility::str(lon) 00152 + "d not in [-180d, 360d]"); 00153 } 00154 00155 bool UTMUPS::CheckCoords(bool utmp, bool northp, real x, real y, 00156 bool mgrslimits, bool throwp) { 00157 // Limits are all multiples of 100km and are all closed on the both ends. 00158 // Failure tests are such that NaNs succeed. 00159 real slop = mgrslimits ? 0 : MGRS::tile_; 00160 int ind = (utmp ? 2 : 0) + (northp ? 1 : 0); 00161 if (x < mineasting_[ind] - slop || x > maxeasting_[ind] + slop) { 00162 if (!throwp) return false; 00163 throw GeographicErr("Easting " + Utility::str(x/1000) + "km not in " 00164 + (mgrslimits ? "MGRS/" : "") 00165 + (utmp ? "UTM" : "UPS") + " range for " 00166 + (northp ? "N" : "S" ) + " hemisphere [" 00167 + Utility::str((mineasting_[ind] - slop)/1000) 00168 + "km, " 00169 + Utility::str((maxeasting_[ind] + slop)/1000) 00170 + "km]"); 00171 } 00172 if (y < minnorthing_[ind] - slop || y > maxnorthing_[ind] + slop) { 00173 if (!throwp) return false; 00174 throw GeographicErr("Northing " + Utility::str(y/1000) + "km not in " 00175 + (mgrslimits ? "MGRS/" : "") 00176 + (utmp ? "UTM" : "UPS") + " range for " 00177 + (northp ? "N" : "S" ) + " hemisphere [" 00178 + Utility::str((minnorthing_[ind] - slop)/1000) 00179 + "km, " 00180 + Utility::str((maxnorthing_[ind] + slop)/1000) 00181 + "km]"); 00182 } 00183 return true; 00184 } 00185 00186 void UTMUPS::DecodeZone(const std::string& zonestr, int& zone, bool& northp) { 00187 unsigned zlen = unsigned(zonestr.size()); 00188 if (zlen == 0) 00189 throw GeographicErr("Empty zone specification"); 00190 if (zlen > 3) 00191 throw GeographicErr("More than 3 characters in zone specification " 00192 + zonestr); 00193 if (zlen == 3 && 00194 toupper(zonestr[0]) == 'I' && 00195 toupper(zonestr[1]) == 'N' && 00196 toupper(zonestr[2]) == 'V') { 00197 zone = INVALID; 00198 northp = false; 00199 return; 00200 } 00201 char hemi = toupper(zonestr[zlen - 1]); 00202 bool northp1 = hemi == 'N'; 00203 if (! (northp1 || hemi == 'S')) 00204 throw GeographicErr(string("Illegal hemisphere letter ") + hemi + " in " 00205 + zonestr + ", specify N or S"); 00206 if (zlen == 1) 00207 zone = UPS; 00208 else { 00209 const char* c = zonestr.c_str(); 00210 char* q; 00211 int zone1 = strtol(c, &q, 10); 00212 if (q == c) 00213 throw GeographicErr("No zone number found in " + zonestr); 00214 if (q - c != int(zlen) - 1) 00215 throw GeographicErr("Extra text " + 00216 zonestr.substr(q - c, int(zlen) - 1 - (q - c)) + 00217 " in UTM/UPS zone " + zonestr); 00218 if (zone1 == UPS) 00219 // Don't allow 0N as an alternative to N for UPS coordinates 00220 throw GeographicErr("Illegal zone 0 in " + zonestr + 00221 ", use just " + hemi + " for UPS"); 00222 if (!(zone1 >= MINUTMZONE && zone1 <= MAXUTMZONE)) 00223 throw GeographicErr("Zone " + Utility::str(zone1) 00224 + " not in range [1, 60]"); 00225 zone = zone1; 00226 } 00227 northp = northp1; 00228 } 00229 00230 std::string UTMUPS::EncodeZone(int zone, bool northp) { 00231 if (zone == INVALID) 00232 return string("INV"); 00233 if (!(zone >= MINZONE && zone <= MAXZONE)) 00234 throw GeographicErr("Zone " + Utility::str(zone) 00235 + " not in range [0, 60]"); 00236 ostringstream os; 00237 if (zone != UPS) 00238 os << setfill('0') << setw(2) << zone; 00239 os << (northp ? 'N' : 'S'); 00240 return os.str(); 00241 } 00242 00243 Math::real UTMUPS::UTMShift() throw() { return real(MGRS::utmNshift_); } 00244 00245 } // namespace GeographicLib