00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <GeographicLib/UTMUPS.hpp>
00011 #include <GeographicLib/MGRS.hpp>
00012 #include <GeographicLib/PolarStereographic.hpp>
00013 #include <GeographicLib/TransverseMercator.hpp>
00014 #include <GeographicLib/Utility.hpp>
00015
00016 namespace GeographicLib {
00017
00018 using namespace std;
00019
00020 const int UTMUPS::falseeasting_[4] =
00021 { MGRS::upseasting_ * MGRS::tile_, MGRS::upseasting_ * MGRS::tile_,
00022 MGRS::utmeasting_ * MGRS::tile_, MGRS::utmeasting_ * MGRS::tile_ };
00023 const int UTMUPS::falsenorthing_[4] =
00024 { MGRS::upseasting_ * MGRS::tile_, MGRS::upseasting_ * MGRS::tile_,
00025 MGRS::maxutmSrow_ * MGRS::tile_, MGRS::minutmNrow_ * MGRS::tile_ };
00026 const int UTMUPS::mineasting_[4] =
00027 { MGRS::minupsSind_ * MGRS::tile_, MGRS::minupsNind_ * MGRS::tile_,
00028 MGRS::minutmcol_ * MGRS::tile_, MGRS::minutmcol_ * MGRS::tile_ };
00029 const int UTMUPS::maxeasting_[4] =
00030 { MGRS::maxupsSind_ * MGRS::tile_, MGRS::maxupsNind_ * MGRS::tile_,
00031 MGRS::maxutmcol_ * MGRS::tile_, MGRS::maxutmcol_ * MGRS::tile_ };
00032 const int UTMUPS::minnorthing_[4] =
00033 { MGRS::minupsSind_ * MGRS::tile_, MGRS::minupsNind_ * MGRS::tile_,
00034 MGRS::minutmSrow_ * MGRS::tile_,
00035 (MGRS::minutmNrow_ + MGRS::minutmSrow_ - MGRS::maxutmSrow_)
00036 * MGRS::tile_ };
00037 const int UTMUPS::maxnorthing_[4] =
00038 { MGRS::maxupsSind_ * MGRS::tile_, MGRS::maxupsNind_ * MGRS::tile_,
00039 (MGRS::maxutmSrow_ + MGRS::maxutmNrow_ - MGRS::minutmNrow_) * MGRS::tile_,
00040 MGRS::maxutmNrow_ * MGRS::tile_ };
00041
00042 int UTMUPS::StandardZone(real lat, real lon, int setzone) {
00043 if (!(setzone >= MINPSEUDOZONE && setzone <= MAXZONE))
00044 throw GeographicErr("Illegal zone requested " + Utility::str(setzone));
00045 if (setzone >= MINZONE || setzone == INVALID)
00046 return setzone;
00047 if (Math::isnan(lat) || Math::isnan(lon))
00048 return INVALID;
00049 if (setzone == UTM || (lat >= -80 && lat < 84)) {
00050
00051 int ilon = int(floor(lon));
00052 if (ilon >= 180)
00053 ilon -= 360;
00054 else if (ilon < -180)
00055 ilon += 360;
00056 int zone = (ilon + 186)/6;
00057 int band = MGRS::LatitudeBand(lat);
00058 if (band == 7 && zone == 31 && ilon >= 3)
00059 zone = 32;
00060 else if (band == 9 && ilon >= 0 && ilon < 42)
00061 zone = 2 * ((ilon + 183)/12) + 1;
00062 return zone;
00063 } else
00064 return UPS;
00065 }
00066
00067 void UTMUPS::Forward(real lat, real lon,
00068 int& zone, bool& northp, real& x, real& y,
00069 real& gamma, real& k,
00070 int setzone, bool mgrslimits) {
00071 CheckLatLon(lat, lon);
00072 bool northp1 = lat >= 0;
00073 int zone1 = StandardZone(lat, lon, setzone);
00074 if (zone1 == INVALID) {
00075 zone = zone1;
00076 northp = northp1;
00077 x = y = gamma = k = Math::NaN();
00078 return;
00079 }
00080 real x1, y1, gamma1, k1;
00081 bool utmp = zone1 != UPS;
00082 if (utmp) {
00083 real
00084 lon0 = CentralMeridian(zone1),
00085 dlon = lon - lon0;
00086 dlon = abs(dlon - 360 * floor((dlon + 180)/360));
00087 if (dlon > 60)
00088
00089
00090 throw GeographicErr("Longitude " + Utility::str(lon)
00091 + "d more than 60d from center of UTM zone "
00092 + Utility::str(zone1));
00093 TransverseMercator::UTM().Forward(lon0, lat, lon, x1, y1, gamma1, k1);
00094 } else {
00095 if (abs(lat) < 70)
00096
00097 throw GeographicErr("Latitude " + Utility::str(lat)
00098 + "d more than 20d from "
00099 + (northp1 ? "N" : "S") + " pole");
00100 PolarStereographic::UPS().Forward(northp1, lat, lon, x1, y1, gamma1, k1);
00101 }
00102 int ind = (utmp ? 2 : 0) + (northp1 ? 1 : 0);
00103 x1 += falseeasting_[ind];
00104 y1 += falsenorthing_[ind];
00105 if (! CheckCoords(zone1 != UPS, northp1, x1, y1, mgrslimits, false) )
00106 throw GeographicErr("Latitude " + Utility::str(lat)
00107 + ", longitude " + Utility::str(lon)
00108 + " out of legal range for "
00109 + (utmp ? "UTM zone " + Utility::str(zone1) : "UPS"));
00110 zone = zone1;
00111 northp = northp1;
00112 x = x1;
00113 y = y1;
00114 gamma = gamma1;
00115 k = k1;
00116 }
00117
00118 void UTMUPS::Reverse(int zone, bool northp, real x, real y,
00119 real& lat, real& lon, real& gamma, real& k,
00120 bool mgrslimits) {
00121 if (zone == INVALID || Math::isnan(x) || Math::isnan(y)) {
00122 lat = lon = gamma = k = Math::NaN();
00123 return;
00124 }
00125 if (!(zone >= MINZONE && zone <= MAXZONE))
00126 throw GeographicErr("Zone " + Utility::str(zone)
00127 + " not in range [0, 60]");
00128 bool utmp = zone != UPS;
00129 CheckCoords(utmp, northp, x, y, mgrslimits);
00130 int ind = (utmp ? 2 : 0) + (northp ? 1 : 0);
00131 x -= falseeasting_[ind];
00132 y -= falsenorthing_[ind];
00133 if (utmp)
00134 TransverseMercator::UTM().Reverse(CentralMeridian(zone),
00135 x, y, lat, lon, gamma, k);
00136 else
00137 PolarStereographic::UPS().Reverse(northp, x, y, lat, lon, gamma, k);
00138 }
00139
00140 void UTMUPS::CheckLatLon(real lat, real lon) {
00141 if (abs(lat) > 90)
00142 throw GeographicErr("Latitude " + Utility::str(lat)
00143 + "d not in [-90d, 90d]");
00144 if (lon < -540 || lon >= 540)
00145 throw GeographicErr("Longitude " + Utility::str(lon)
00146 + "d not in [-540d, 540d)");
00147 }
00148
00149 bool UTMUPS::CheckCoords(bool utmp, bool northp, real x, real y,
00150 bool mgrslimits, bool throwp) {
00151
00152
00153 real slop = mgrslimits ? 0 : MGRS::tile_;
00154 int ind = (utmp ? 2 : 0) + (northp ? 1 : 0);
00155 if (x < mineasting_[ind] - slop || x > maxeasting_[ind] + slop) {
00156 if (!throwp) return false;
00157 throw GeographicErr("Easting " + Utility::str(x/1000) + "km not in "
00158 + (mgrslimits ? "MGRS/" : "")
00159 + (utmp ? "UTM" : "UPS") + " range for "
00160 + (northp ? "N" : "S" ) + " hemisphere ["
00161 + Utility::str((mineasting_[ind] - slop)/1000)
00162 + "km, "
00163 + Utility::str((maxeasting_[ind] + slop)/1000)
00164 + "km]");
00165 }
00166 if (y < minnorthing_[ind] - slop || y > maxnorthing_[ind] + slop) {
00167 if (!throwp) return false;
00168 throw GeographicErr("Northing " + Utility::str(y/1000) + "km not in "
00169 + (mgrslimits ? "MGRS/" : "")
00170 + (utmp ? "UTM" : "UPS") + " range for "
00171 + (northp ? "N" : "S" ) + " hemisphere ["
00172 + Utility::str((minnorthing_[ind] - slop)/1000)
00173 + "km, "
00174 + Utility::str((maxnorthing_[ind] + slop)/1000)
00175 + "km]");
00176 }
00177 return true;
00178 }
00179
00180 void UTMUPS::Transfer(int zonein, bool northpin, real xin, real yin,
00181 int zoneout, bool northpout, real& xout, real& yout,
00182 int& zone) {
00183 bool northp = northpin;
00184 if (zonein != zoneout) {
00185
00186 real lat, lon;
00187 GeographicLib::UTMUPS::Reverse(zonein, northpin, xin, yin, lat, lon);
00188
00189 real x, y;
00190 int zone1;
00191 GeographicLib::UTMUPS::Forward(lat, lon, zone1, northp, x, y,
00192 zoneout == UTMUPS::MATCH
00193 ? zonein : zoneout);
00194 if (zone1 == 0 && northp != northpout)
00195 throw GeographicErr
00196 ("Attempt to transfer UPS coordinates between hemispheres");
00197 zone = zone1;
00198 xout = x;
00199 yout = y;
00200 } else {
00201 if (zoneout == 0 && northp != northpout)
00202 throw GeographicErr
00203 ("Attempt to transfer UPS coordinates between hemispheres");
00204 zone = zoneout;
00205 xout = xin;
00206 yout = yin;
00207 }
00208 if (northp != northpout)
00209
00210 yout += (northpout ? -1 : 1) * MGRS::utmNshift_;
00211 return;
00212 }
00213
00214 void UTMUPS::DecodeZone(const std::string& zonestr, int& zone, bool& northp) {
00215 unsigned zlen = unsigned(zonestr.size());
00216 if (zlen == 0)
00217 throw GeographicErr("Empty zone specification");
00218
00219 if (zlen > 7)
00220 throw GeographicErr("More than 7 characters in zone specification "
00221 + zonestr);
00222
00223 const char* c = zonestr.c_str();
00224 char* q;
00225 int zone1 = strtol(c, &q, 10);
00226
00227
00228 if (zone1 == UPS) {
00229 if (!(q == c))
00230
00231 throw GeographicErr("Illegal zone 0 in " + zonestr +
00232 ", use just the hemisphere for UPS");
00233 } else if (!(zone1 >= MINUTMZONE && zone1 <= MAXUTMZONE))
00234 throw GeographicErr("Zone " + Utility::str(zone1)
00235 + " not in range [1, 60]");
00236 else if (!isdigit(zonestr[0]))
00237 throw GeographicErr("Must use unsigned number for zone "
00238 + Utility::str(zone1));
00239 else if (q - c > 2)
00240 throw GeographicErr("More than 2 digits use to specify zone "
00241 + Utility::str(zone1));
00242
00243 string hemi = zonestr.substr(q - c);
00244 transform(hemi.begin(), hemi.end(), hemi.begin(), (int(*)(int))tolower);
00245 if (q == c && (hemi == "inv" || hemi == "invalid")) {
00246 zone = INVALID;
00247 northp = false;
00248 return;
00249 }
00250 bool northp1 = hemi == "north" || hemi == "n";
00251 if (!(northp1 || hemi == "south" || hemi == "s"))
00252 throw GeographicErr(string("Illegal hemisphere ") + hemi + " in "
00253 + zonestr + ", specify north or south");
00254 zone = zone1;
00255 northp = northp1;
00256 }
00257
00258 std::string UTMUPS::EncodeZone(int zone, bool northp, bool abbrev) {
00259 if (zone == INVALID)
00260 return string(abbrev ? "inv" : "invalid");
00261 if (!(zone >= MINZONE && zone <= MAXZONE))
00262 throw GeographicErr("Zone " + Utility::str(zone)
00263 + " not in range [0, 60]");
00264 ostringstream os;
00265 if (zone != UPS)
00266 os << setfill('0') << setw(2) << zone;
00267 if (abbrev)
00268 os << (northp ? 'n' : 's');
00269 else
00270 os << (northp ? "north" : "south");
00271 return os.str();
00272 }
00273
00274 void UTMUPS::DecodeEPSG(int epsg, int& zone, bool& northp) {
00275 northp = false;
00276 if (epsg >= epsg01N && epsg <= epsg60N) {
00277 zone = (epsg - epsg01N) + MINUTMZONE;
00278 northp = true;
00279 } else if (epsg == epsgN) {
00280 zone = UPS;
00281 northp = true;
00282 } else if (epsg >= epsg01S && epsg <= epsg60S) {
00283 zone = (epsg - epsg01S) + MINUTMZONE;
00284 } else if (epsg == epsgS) {
00285 zone = UPS;
00286 } else {
00287 zone = INVALID;
00288 }
00289 }
00290
00291 int UTMUPS::EncodeEPSG(int zone, bool northp) {
00292 int epsg = -1;
00293 if (zone == UPS)
00294 epsg = epsgS;
00295 else if (zone >= MINUTMZONE && zone <= MAXUTMZONE)
00296 epsg = (zone - MINUTMZONE) + epsg01S;
00297 if (epsg >= 0 && northp)
00298 epsg += epsgN - epsgS;
00299 return epsg;
00300 }
00301
00302 Math::real UTMUPS::UTMShift() { return real(MGRS::utmNshift_); }
00303
00304 }