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 <iomanip>
00015
00016 #define GEOGRAPHICLIB_UTMUPS_CPP "$Id: UTMUPS.cpp 6918 2010-12-21 12:56:07Z karney $"
00017
00018 RCSID_DECL(GEOGRAPHICLIB_UTMUPS_CPP)
00019 RCSID_DECL(GEOGRAPHICLIB_UTMUPS_HPP)
00020
00021 namespace GeographicLib {
00022
00023 using namespace std;
00024
00025 const Math::real UTMUPS::falseeasting[4] =
00026 { MGRS::upseasting * MGRS::tile, MGRS::upseasting * MGRS::tile,
00027 MGRS::utmeasting * MGRS::tile, MGRS::utmeasting* MGRS::tile };
00028 const Math::real UTMUPS::falsenorthing[4] =
00029 { MGRS::upseasting * MGRS::tile, MGRS::upseasting * MGRS::tile,
00030 MGRS::maxutmSrow * MGRS::tile, MGRS::minutmNrow * MGRS::tile };
00031 const Math::real UTMUPS::mineasting[4] =
00032 { MGRS::minupsSind * MGRS::tile, MGRS::minupsNind * MGRS::tile,
00033 MGRS::minutmcol * MGRS::tile, MGRS::minutmcol * MGRS::tile };
00034 const Math::real UTMUPS::maxeasting[4] =
00035 { MGRS::maxupsSind * MGRS::tile, MGRS::maxupsNind * MGRS::tile,
00036 MGRS::maxutmcol * MGRS::tile, MGRS::maxutmcol * MGRS::tile };
00037 const Math::real UTMUPS::minnorthing[4] =
00038 { MGRS::minupsSind * MGRS::tile, MGRS::minupsNind * MGRS::tile,
00039 MGRS::minutmSrow * MGRS::tile,
00040 (MGRS::minutmNrow + MGRS::minutmSrow - MGRS::maxutmSrow) * MGRS::tile };
00041 const Math::real UTMUPS::maxnorthing[4] =
00042 { MGRS::maxupsSind * MGRS::tile, MGRS::maxupsNind * MGRS::tile,
00043 (MGRS::maxutmSrow + MGRS::maxutmNrow - MGRS::minutmNrow) * MGRS::tile,
00044 MGRS::maxutmNrow * MGRS::tile };
00045
00046 int UTMUPS::StandardZone(real lat, real lon, int setzone) {
00047 if (setzone < MINPSEUDOZONE || setzone > MAXZONE)
00048 throw GeographicErr("Illegal zone requested " + str(setzone));
00049 if (setzone >= MINZONE || setzone == INVALID)
00050 return setzone;
00051 if (Math::isnan(lat) || Math::isnan(lon))
00052 return INVALID;
00053
00054 if (setzone == UTM || (lat >= -80 && lat < 84)) {
00055
00056 int ilon = int(floor(lon));
00057 if (ilon >= 180)
00058 ilon -= 360;
00059 int zone = (ilon + 186)/6;
00060 int band = MGRS::LatitudeBand(lat);
00061 if (band == 7 && zone == 31 && ilon >= 3)
00062 zone = 32;
00063 else if (band == 9 && ilon >= 0 && ilon < 42)
00064 zone = 2 * ((ilon + 183)/12) + 1;
00065 return zone;
00066 } else
00067 return UPS;
00068 }
00069
00070 void UTMUPS::Forward(real lat, real lon,
00071 int& zone, bool& northp, real& x, real& y,
00072 real& gamma, real& k,
00073 int setzone, bool mgrslimits) {
00074 CheckLatLon(lat, lon);
00075 bool northp1 = lat >= 0;
00076 int zone1 = StandardZone(lat, lon, setzone);
00077 if (zone1 == INVALID) {
00078 zone = zone1;
00079 northp = northp1;
00080 x = y = gamma = k = Math::NaN();
00081 return;
00082 }
00083 real x1, y1, gamma1, k1;
00084 bool utmp = zone1 != UPS;
00085 if (utmp) {
00086 real
00087 lon0 = CentralMeridian(zone1),
00088 dlon = lon - lon0;
00089 dlon = abs(dlon - 360 * floor((dlon + 180)/360));
00090 if (dlon > 60)
00091
00092
00093 throw GeographicErr("Longitude " + str(lon)
00094 + "d more than 60d from center of UTM zone "
00095 + str(zone1));
00096 TransverseMercator::UTM.Forward(lon0, lat, lon, x1, y1, gamma1, k1);
00097 } else {
00098 if (abs(lat) < 70)
00099
00100 throw GeographicErr("Latitude " + str(lat) + "d more than 20d from "
00101 + (northp1 ? "N" : "S") + " pole");
00102 PolarStereographic::UPS.Forward(northp1, lat, lon, x1, y1, gamma1, k1);
00103 }
00104 int ind = (utmp ? 2 : 0) + (northp1 ? 1 : 0);
00105 x1 += falseeasting[ind];
00106 y1 += falsenorthing[ind];
00107 if (! CheckCoords(zone1 != UPS, northp1, x1, y1, mgrslimits, false) )
00108 throw GeographicErr("Latitude " + str(lat) + ", longitude " + str(lon)
00109 + " out of legal range for "
00110 + (utmp ? "UTM zone " + str(zone1) : "UPS"));
00111 zone = zone1;
00112 northp = northp1;
00113 x = x1;
00114 y = y1;
00115 gamma = gamma1;
00116 k = k1;
00117 }
00118
00119 void UTMUPS::Reverse(int zone, bool northp, real x, real y,
00120 real& lat, real& lon, real& gamma, real& k,
00121 bool mgrslimits) {
00122 if (zone == INVALID || Math::isnan(x) || Math::isnan(y)) {
00123 lat = lon = gamma = k = Math::NaN();
00124 return;
00125 }
00126 if (! (zone >= MINZONE && zone <= MAXZONE))
00127 throw GeographicErr("Zone " + str(zone) + " 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 (lat < -90 || lat > 90)
00142 throw GeographicErr("Latitude " + str(lat) + "d not in [-90d, 90d]");
00143 if (lon < -180 || lon > 360)
00144 throw GeographicErr("Latitude " + str(lon) + "d not in [-180d, 360d]");
00145 }
00146
00147 bool UTMUPS::CheckCoords(bool utmp, bool northp, real x, real y,
00148 bool mgrslimits, bool throwp) {
00149
00150
00151 real slop = mgrslimits ? 0 : MGRS::tile;
00152 int ind = (utmp ? 2 : 0) + (northp ? 1 : 0);
00153 if (x < mineasting[ind] - slop || x > maxeasting[ind] + slop) {
00154 if (!throwp) return false;
00155 throw GeographicErr("Easting " + str(x/1000) + "km not in "
00156 + (mgrslimits ? "MGRS/" : "")
00157 + (utmp ? "UTM" : "UPS") + " range for "
00158 + (northp ? "N" : "S" ) + " hemisphere ["
00159 + str((mineasting[ind] - slop)/1000) + "km, "
00160 + str((maxeasting[ind] + slop)/1000) + "km]");
00161 }
00162 if (y < minnorthing[ind] - slop || y > maxnorthing[ind] + slop) {
00163 if (!throwp) return false;
00164 throw GeographicErr("Northing " + str(y/1000) + "km not in "
00165 + (mgrslimits ? "MGRS/" : "")
00166 + (utmp ? "UTM" : "UPS") + " range for "
00167 + (northp ? "N" : "S" ) + " hemisphere ["
00168 + str((minnorthing[ind] - slop)/1000) + "km, "
00169 + str((maxnorthing[ind] + slop)/1000) + "km]");
00170 }
00171 return true;
00172 }
00173
00174 void UTMUPS::DecodeZone(const std::string& zonestr, int& zone, bool& northp) {
00175 unsigned zlen = unsigned(zonestr.size());
00176 if (zlen == 0)
00177 throw GeographicErr("Empty zone specification");
00178 if (zlen > 3)
00179 throw GeographicErr("More than 3 characters in zone specification "
00180 + zonestr);
00181 if (zlen == 3 &&
00182 toupper(zonestr[0]) == 'I' &&
00183 toupper(zonestr[1]) == 'N' &&
00184 toupper(zonestr[2]) == 'V') {
00185 zone = INVALID;
00186 northp = false;
00187 return;
00188 }
00189 char hemi = toupper(zonestr[zlen - 1]);
00190 bool northp1 = hemi == 'N';
00191 if (! (northp1 || hemi == 'S'))
00192 throw GeographicErr(string("Illegal hemisphere letter ") + hemi + " in "
00193 + zonestr + ", specify N or S");
00194 if (zlen == 1)
00195 zone = UPS;
00196 else {
00197 const char* c = zonestr.c_str();
00198 char* q;
00199 int zone1 = strtol(c, &q, 10);
00200 if (q == c)
00201 throw GeographicErr("No zone number found in " + zonestr);
00202 if (q - c != int(zlen) - 1)
00203 throw GeographicErr("Extra text " +
00204 zonestr.substr(q - c, int(zlen) - 1 - (q - c)) +
00205 " in UTM/UPS zone " + zonestr);
00206 if (zone1 == UPS)
00207
00208 throw GeographicErr("Illegal zone 0 in " + zonestr +
00209 ", use just " + hemi + " for UPS");
00210 if (!(zone1 >= MINUTMZONE && zone1 <= MAXUTMZONE))
00211 throw GeographicErr("Zone " + str(zone1) + " not in range [1, 60]");
00212 zone = zone1;
00213 }
00214 northp = northp1;
00215 }
00216
00217 std::string UTMUPS::EncodeZone(int zone, bool northp) {
00218 if (zone == INVALID)
00219 return string("INV");
00220 if (! (zone >= MINZONE && zone <= MAXZONE))
00221 throw GeographicErr("Zone " + str(zone) + " not in range [0, 60]");
00222 ostringstream os;
00223 if (zone != UPS)
00224 os << setfill('0') << setw(2) << zone;
00225 os << (northp ? 'N' : 'S');
00226 return os.str();
00227 }
00228
00229 Math::real UTMUPS::UTMShift() throw() { return real(MGRS::utmNshift); }
00230
00231 }