GeographicLib  1.21
UTMUPS.cpp
Go to the documentation of this file.
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