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, 2009, 2010) <charles@karney.com>
00006  * and licensed under the LGPL.  For more information, see
00007  * http://geographiclib.sourceforge.net/
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)) // Check if lat or lon is a NaN
00052       return INVALID;
00053     // Assume lon is in [-180, 360].
00054     if (setzone == UTM || (lat >= -80 && lat < 84)) {
00055       // Assume lon is in [-180, 360].
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         // Check isn't really necessary because CheckCoords catches this case.
00092         // But this allows a more meaningful error message to be given.
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         // Check isn't really necessary ... (see above).
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     // Limits are all multiples of 100km and are all closed on the both ends.
00150     // Failure tests are such that NaNs succeed.
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         // Don't allow 0N as an alternative to N for UPS coordinates
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 } // namespace GeographicLib