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 6814 2010-02-04 22:33:41Z 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)
00050 return setzone;
00051
00052 if (setzone == UTM || (lat >= -80 && lat < 84)) {
00053
00054 int ilon = int(floor(lon));
00055 if (ilon >= 180)
00056 ilon -= 360;
00057 int zone = (ilon + 186)/6;
00058 int band = MGRS::LatitudeBand(lat);
00059 if (band == 7 && zone == 31 && ilon >= 3)
00060 zone = 32;
00061 else if (band == 9 && ilon >= 0 && ilon < 42)
00062 zone = 2 * ((ilon + 183)/12) + 1;
00063 return zone;
00064 } else
00065 return UPS;
00066 }
00067
00068 void UTMUPS::Forward(real lat, real lon,
00069 int& zone, bool& northp, real& x, real& y,
00070 real& gamma, real& k,
00071 int setzone, bool mgrslimits) {
00072 CheckLatLon(lat, lon);
00073 bool northp1 = lat >= 0;
00074 int zone1 = StandardZone(lat, lon, setzone);
00075 real x1, y1, gamma1, k1;
00076 bool utmp = zone1 != UPS;
00077 if (utmp) {
00078 real
00079 lon0 = CentralMeridian(zone1),
00080 dlon = lon - lon0;
00081 dlon = abs(dlon - 360 * floor((dlon + 180)/360));
00082 if (dlon > 60)
00083
00084
00085 throw GeographicErr("Longitude " + str(lon)
00086 + "d more than 60d from center of UTM zone "
00087 + str(zone1));
00088 TransverseMercator::UTM.Forward(lon0, lat, lon, x1, y1, gamma1, k1);
00089 } else {
00090 if (abs(lat) < 70)
00091
00092 throw GeographicErr("Latitude " + str(lat) + "d more than 20d from "
00093 + (northp1 ? "N" : "S") + " pole");
00094 PolarStereographic::UPS.Forward(northp1, lat, lon, x1, y1, gamma1, k1);
00095 }
00096 int ind = (utmp ? 2 : 0) + (northp1 ? 1 : 0);
00097 x1 += falseeasting[ind];
00098 y1 += falsenorthing[ind];
00099 if (! CheckCoords(zone1 != UPS, northp1, x1, y1, mgrslimits, false) )
00100 throw GeographicErr("Latitude " + str(lat) + ", longitude " + str(lon)
00101 + " out of legal range for "
00102 + (utmp ? "UTM zone " + str(zone1) : "UPS"));
00103 zone = zone1;
00104 northp = northp1;
00105 x = x1;
00106 y = y1;
00107 gamma = gamma1;
00108 k = k1;
00109 }
00110
00111 void UTMUPS::Reverse(int zone, bool northp, real x, real y,
00112 real& lat, real& lon, real& gamma, real& k,
00113 bool mgrslimits) {
00114 if (! (zone >= MINZONE && zone <= MAXZONE))
00115 throw GeographicErr("Zone " + str(zone) + " not in range [0, 60]");
00116 bool utmp = zone != UPS;
00117 CheckCoords(utmp, northp, x, y, mgrslimits);
00118 int ind = (utmp ? 2 : 0) + (northp ? 1 : 0);
00119 x -= falseeasting[ind];
00120 y -= falsenorthing[ind];
00121 if (utmp)
00122 TransverseMercator::UTM.Reverse(CentralMeridian(zone),
00123 x, y, lat, lon, gamma, k);
00124 else
00125 PolarStereographic::UPS.Reverse(northp, x, y, lat, lon, gamma, k);
00126 }
00127
00128 void UTMUPS::CheckLatLon(real lat, real lon) {
00129 if (! (lat >= -90 && lat <= 90))
00130 throw GeographicErr("Latitude " + str(lat) + "d not in [-90d, 90d]");
00131 if (! (lon >= -180 && lon <= 360))
00132 throw GeographicErr("Latitude " + str(lon) + "d not in [-180d, 360d]");
00133 }
00134
00135 bool UTMUPS::CheckCoords(bool utmp, bool northp, real x, real y,
00136 bool mgrslimits, bool throwp) {
00137
00138
00139 real slop = mgrslimits ? 0 : MGRS::tile;
00140 int ind = (utmp ? 2 : 0) + (northp ? 1 : 0);
00141 if (! (x >= mineasting[ind] - slop && x <= maxeasting[ind] + slop) ) {
00142 if (!throwp) return false;
00143 throw GeographicErr("Easting " + str(x/1000) + "km not in "
00144 + (mgrslimits ? "MGRS/" : "")
00145 + (utmp ? "UTM" : "UPS") + " range for "
00146 + (northp ? "N" : "S" ) + " hemisphere ["
00147 + str((mineasting[ind] - slop)/1000) + "km, "
00148 + str((maxeasting[ind] + slop)/1000) + "km]");
00149 }
00150 if (! (y >= minnorthing[ind] - slop && y <= maxnorthing[ind] + slop) ) {
00151 if (!throwp) return false;
00152 throw GeographicErr("Northing " + str(y/1000) + "km not in "
00153 + (mgrslimits ? "MGRS/" : "")
00154 + (utmp ? "UTM" : "UPS") + " range for "
00155 + (northp ? "N" : "S" ) + " hemisphere ["
00156 + str((minnorthing[ind] - slop)/1000) + "km, "
00157 + str((maxnorthing[ind] + slop)/1000) + "km]");
00158 }
00159 return true;
00160 }
00161
00162 void UTMUPS::DecodeZone(const std::string& zonestr, int& zone, bool& northp) {
00163 unsigned zlen = unsigned(zonestr.size());
00164 if (zlen == 0)
00165 throw GeographicErr("Empty zone specification");
00166 if (zlen > 3)
00167 throw GeographicErr("More than 3 characters in zone specification "
00168 + zonestr);
00169 char hemi = toupper(zonestr[zlen - 1]);
00170 bool northp1 = hemi == 'N';
00171 if (! (northp1 || hemi == 'S'))
00172 throw GeographicErr(string("Illegal hemisphere letter ") + hemi + " in "
00173 + zonestr + ", specify N or S");
00174 if (zlen == 1)
00175 zone = UPS;
00176 else {
00177 const char* c = zonestr.c_str();
00178 char* q;
00179 int zone1 = strtol(c, &q, 10);
00180 if (q == c)
00181 throw GeographicErr("No zone number found in " + zonestr);
00182 if (q - c != int(zlen) - 1)
00183 throw GeographicErr("Extra text " +
00184 zonestr.substr(q - c, int(zlen) - 1 - (q - c)) +
00185 " in UTM/UPS zone " + zonestr);
00186 if (zone1 == UPS)
00187
00188 throw GeographicErr("Illegal zone 0 in " + zonestr +
00189 ", use just " + hemi + " for UPS");
00190 if (!(zone1 >= MINUTMZONE && zone1 <= MAXUTMZONE))
00191 throw GeographicErr("Zone " + str(zone1) + " not in range [1, 60]");
00192 zone = zone1;
00193 }
00194 northp = northp1;
00195 }
00196
00197 std::string UTMUPS::EncodeZone(int zone, bool northp) {
00198 if (! (zone >= MINZONE && zone <= MAXZONE))
00199 throw GeographicErr("Zone " + str(zone) + " not in range [0, 60]");
00200 ostringstream os;
00201 if (zone != UPS)
00202 os << setfill('0') << setw(2) << zone;
00203 os << (northp ? 'N' : 'S');
00204 return os.str();
00205 }
00206
00207 Math::real UTMUPS::UTMShift() throw() { return real(MGRS::utmNshift); }
00208
00209 }