DMS.cpp

Go to the documentation of this file.
00001 /**
00002  * \file DMS.cpp
00003  * \brief Implementation for GeographicLib::DMS class
00004  *
00005  * Copyright (c) Charles Karney (2008, 2009, 2010, 2011) <charles@karney.com>
00006  * and licensed under the LGPL.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  **********************************************************************/
00009 
00010 #include "GeographicLib/DMS.hpp"
00011 #include <algorithm>
00012 
00013 #define GEOGRAPHICLIB_DMS_CPP "$Id: DMS.cpp 6956 2011-02-17 23:20:13Z karney $"
00014 
00015 RCSID_DECL(GEOGRAPHICLIB_DMS_CPP)
00016 RCSID_DECL(GEOGRAPHICLIB_DMS_HPP)
00017 
00018 namespace GeographicLib {
00019 
00020   using namespace std;
00021 
00022   const string DMS::hemispheres = "SNWE";
00023   const string DMS::signs = "-+";
00024   const string DMS::digits = "0123456789";
00025   const string DMS::dmsindicators = "D'\"";
00026   const string DMS::components[] = {"degrees", "minutes", "seconds"};
00027 
00028   Math::real DMS::Decode(const std::string& dms, flag& ind) {
00029     try {
00030       int sign = 1;
00031       unsigned
00032         beg = 0,
00033         end = unsigned(dms.size());
00034       while (beg < end && isspace(dms[beg]))
00035         ++beg;
00036       while (beg < end && isspace(dms[end - 1]))
00037         --end;
00038       flag ind1 = NONE;
00039       int k = -1;
00040       if (end > beg && (k = lookup(hemispheres, dms[beg])) >= 0) {
00041         ind1 = (k / 2) ? LONGITUDE : LATITUDE;
00042         sign = k % 2 ? 1 : -1;
00043         ++beg;
00044       }
00045       if (end > beg && (k = lookup(hemispheres, dms[end-1])) >= 0) {
00046         if (k >= 0) {
00047           if (ind1 != NONE) {
00048             if (toupper(dms[beg - 1]) == toupper(dms[end - 1]))
00049               throw GeographicErr("Repeated hemisphere indicators "
00050                                   + str(dms[beg - 1]) + " in "
00051                                   + dms.substr(beg - 1, end - beg + 1));
00052             else
00053               throw GeographicErr("Contradictory hemisphere indicators "
00054                                   + str(dms[beg - 1]) + " and "
00055                                   + str(dms[end - 1]) + " in "
00056                                   + dms.substr(beg - 1, end - beg + 1));
00057           }
00058           ind1 = (k / 2) ? LONGITUDE : LATITUDE;
00059           sign = k % 2 ? 1 : -1;
00060           --end;
00061         }
00062       }
00063       if (end > beg && (k = lookup(signs, dms[beg])) >= 0) {
00064         if (k >= 0) {
00065           sign *= k ? 1 : -1;
00066           ++beg;
00067         }
00068       }
00069       if (end == beg)
00070         throw GeographicErr("Empty or incomplete DMS string " + dms);
00071       real ipieces[] = {0, 0, 0};
00072       real fpieces[] = {0, 0, 0};
00073       unsigned npiece = 0;
00074       real icurrent = 0;
00075       real fcurrent = 0;
00076       unsigned ncurrent = 0, p = beg;
00077       bool pointseen = false;
00078       unsigned digcount = 0;
00079       while (p < end) {
00080         char x = dms[p++];
00081         if ((k = lookup(digits, x)) >= 0) {
00082           ++ncurrent;
00083           if (digcount > 0)
00084             ++digcount;           // Count of decimal digits
00085           else
00086             icurrent = 10 * icurrent + k;
00087         } else if (x == '.') {
00088           if (pointseen)
00089             throw GeographicErr("Multiple decimal points in "
00090                                 + dms.substr(beg, end - beg));
00091           pointseen = true;
00092           digcount = 1;
00093         } else if ((k = lookup(dmsindicators, x)) >= 0) {
00094           if (unsigned(k) == npiece - 1)
00095             throw GeographicErr("Repeated " + components[k] +
00096                                 " component in " + dms.substr(beg, end - beg));
00097           else if (unsigned(k) < npiece)
00098             throw GeographicErr(components[k] + " component follows "
00099                                 + components[npiece - 1] + " component in "
00100                                 + dms.substr(beg, end - beg));
00101           if (ncurrent == 0)
00102             throw GeographicErr("Missing numbers in " + components[k] +
00103                                 " component of " + dms.substr(beg, end - beg));
00104           if (digcount > 1) {
00105             istringstream s(dms.substr(p - digcount - 1, digcount));
00106             s >> fcurrent;
00107           }
00108           ipieces[k] = icurrent;
00109           fpieces[k] = icurrent + fcurrent;
00110           if (p < end) {
00111             npiece = k + 1;
00112             icurrent = fcurrent = 0;
00113             ncurrent = digcount = 0;
00114           }
00115         } else if (lookup(signs, x) >= 0)
00116           throw GeographicErr("Internal sign in DMS string "
00117                               + dms.substr(beg, end - beg));
00118         else
00119           throw GeographicErr("Illegal character " + str(x)
00120                               + " in DMS string "
00121                               + dms.substr(beg, end - beg));
00122       }
00123       if (lookup(dmsindicators, dms[p - 1]) < 0) {
00124         if (npiece >= 3)
00125           throw GeographicErr("Extra text following seconds in DMS string "
00126                               + dms.substr(beg, end - beg));
00127         if (ncurrent == 0)
00128           throw GeographicErr("Missing numbers in " + components[k]
00129                               + " component of " + dms.substr(beg, end - beg));
00130         if (digcount > 1) {
00131           istringstream s(dms.substr(p - digcount, digcount));
00132           s >> fcurrent;
00133         }
00134         ipieces[npiece] = icurrent;
00135         fpieces[npiece] = icurrent + fcurrent;
00136       }
00137       if (pointseen && digcount == 0)
00138         throw GeographicErr("Decimal point in non-terminal component of "
00139                             + dms.substr(beg, end - beg));
00140       // Note that we accept 59.999999... even though it rounds to 60.
00141       if (ipieces[1] >= 60)
00142         throw GeographicErr("Minutes " + str(fpieces[1])
00143                             + " not in range [0, 60)");
00144       if (ipieces[2] >= 60)
00145         throw GeographicErr("Seconds " + str(fpieces[2])
00146                             + " not in range [0, 60)");
00147       ind = ind1;
00148       // Assume check on range of result is made by calling routine (which
00149       // might be able to offer a better diagnostic).
00150       return real(sign) * (fpieces[0] + (fpieces[1] + fpieces[2] / 60) / 60);
00151     }
00152     catch (const GeographicErr&) {
00153       real res = NumMatch(dms);
00154       if (res == 0)
00155         throw;
00156       else
00157         ind = NONE;
00158       return res;
00159     }
00160   }
00161 
00162   Math::real DMS::NumMatch(const std::string& s) {
00163     if (s.length() < 3)
00164       return 0;
00165     string t;
00166     t.resize(s.length());
00167     for (size_t i = s.length(); i--;)
00168       t[i] = toupper(s[i]);
00169     int sign = t[0] == '-' ? -1 : 1;
00170     string::size_type p0 = t[0] == '-' || t[0] == '+' ? 1 : 0;
00171     string::size_type p1 = t.find_last_not_of('0');
00172     if (p1 == string::npos || p1 + 1 < p0 + 3)
00173       return 0;
00174     // Strip off sign and trailing 0s
00175     t = t.substr(p0, p1 + 1 - p0);  // Length at least 3
00176     if (t == "NAN" || t == "1.#QNAN" || t == "1.#SNAN" || t == "1.#IND")
00177       return sign * Math::NaN();
00178     else if (t == "INF" || t == "1.#INF")
00179       return sign * Math::infinity();
00180     return 0;
00181   }
00182 
00183   Math::real DMS::Decode(const std::string& str) {
00184     istringstream is(str);
00185     real num;
00186     try {
00187       if (!(is >> num))
00188         throw GeographicErr("Could not read number: " + str);
00189       // On some platforms, is >> num gobbles final E in 1234E, so look for
00190       // last character which is legal as the final character in a number
00191       // (digit or period).
00192       int pos = min(int(is.tellg()), int(str.find_last_of("0123456789.")) + 1);
00193       if (pos != int(str.size()))
00194         throw GeographicErr("Extra text " + str.substr(pos) +
00195                             " in number " + str);
00196     }
00197     catch (const GeographicErr&) {
00198       num = NumMatch(str);
00199       if (num == 0)
00200         throw;
00201     }
00202     return num;
00203   }
00204 
00205   void DMS::DecodeLatLon(const std::string& stra, const std::string& strb,
00206                          real& lat, real& lon) {
00207       real a, b;
00208       flag ia, ib;
00209       a = Decode(stra, ia);
00210       b = Decode(strb, ib);
00211       if (ia == NONE && ib == NONE) {
00212         // Default to lat, long
00213         ia = LATITUDE;
00214         ib = LONGITUDE;
00215       } else if (ia == NONE)
00216         ia = flag(LATITUDE + LONGITUDE - ib);
00217       else if (ib == NONE)
00218         ib = flag(LATITUDE + LONGITUDE - ia);
00219       if (ia == ib)
00220         throw GeographicErr("Both " + stra + " and "
00221                             + strb + " interpreted as "
00222                             + (ia == LATITUDE ? "latitudes" : "longitudes"));
00223       real
00224         lat1 = ia == LATITUDE ? a : b,
00225         lon1 = ia == LATITUDE ? b : a;
00226       if (lat1 < -90 || lat1 > 90)
00227         throw GeographicErr("Latitude " + str(lat1) + "d not in [-90d, 90d]");
00228       if (lon1 < -180 || lon1 > 360)
00229         throw GeographicErr("Latitude " + str(lon1)
00230                             + "d not in [-180d, 360d]");
00231       if (lon1 >= 180)
00232         lon1 -= 360;
00233       lat = lat1;
00234       lon = lon1;
00235   }
00236 
00237   Math::real DMS::DecodeAngle(const std::string& angstr) {
00238     DMS::flag ind;
00239     real ang = Decode(angstr, ind);
00240     if (ind != DMS::NONE)
00241       throw GeographicErr("Arc angle " + angstr
00242                           + " includes a hemisphere, N/E/W/S");
00243     return ang;
00244   }
00245 
00246   Math::real DMS::DecodeAzimuth(const std::string& azistr) {
00247     DMS::flag ind;
00248     real azi = Decode(azistr, ind);
00249     if (ind == DMS::LATITUDE)
00250       throw GeographicErr("Azimuth " + azistr
00251                           + " has a latitude hemisphere, N/S");
00252     if (azi < -180 || azi > 360)
00253       throw GeographicErr("Azimuth " + azistr + " not in range [-180d,360d]");
00254     if (azi >= 180) azi -= 360;
00255     return azi;
00256   }
00257 
00258   string DMS::Encode(real angle, component trailing, unsigned prec, flag ind) {
00259     // Assume check on range of input angle has been made by calling
00260     // routine (which might be able to offer a better diagnostic).
00261     if (!Math::isfinite(angle))
00262       return angle < 0 ? string("-inf") :
00263         (angle > 0 ? string("inf") : string("nan"));
00264 
00265     // 15 - 2 * trailing = ceiling(log10(2^53/90/60^trailing)).
00266     // This suffices to give full real precision for numbers in [-90,90]
00267     prec = min(15 - 2 * unsigned(trailing), prec);
00268     real scale = 1;
00269     for (unsigned i = 0; i < unsigned(trailing); ++i)
00270       scale *= 60;
00271     for (unsigned i = 0; i < prec; ++i)
00272       scale *= 10;
00273     if (ind == AZIMUTH)
00274       angle -= floor(angle/360) * 360;
00275     int sign = angle < 0 ? -1 : 1;
00276     angle *= sign;
00277 
00278     // Break off integer part to preserve precision in manipulation of
00279     // fractional part.
00280     real
00281       idegree = floor(angle),
00282       fdegree = floor((angle - idegree) * scale + real(0.5)) / scale;
00283     if (fdegree >= 1) {
00284       idegree += 1;
00285       fdegree -= 1;
00286     }
00287     real pieces[3] = {fdegree, 0, 0};
00288     for (unsigned i = 1; i <= unsigned(trailing); ++i) {
00289       real
00290         ip = floor(pieces[i - 1]),
00291         fp = pieces[i - 1] - ip;
00292       pieces[i] = fp * 60;
00293       pieces[i - 1] = ip;
00294     }
00295     pieces[0] += idegree;
00296     ostringstream s;
00297     s << fixed  << setfill('0');
00298     if (ind == NONE && sign < 0)
00299       s << '-';
00300     switch (trailing) {
00301     case DEGREE:
00302       if (ind != NONE)
00303         s << setw(1 + min(int(ind), 2) + prec + (prec ? 1 : 0));
00304       s << setprecision(prec) << pieces[0];
00305       // Don't include degree designator (d) if it is the trailing component.
00306       break;
00307     default:
00308       if (ind != NONE)
00309         s << setw(1 + min(int(ind), 2));
00310       s << setprecision(0) << pieces[0] << char(tolower(dmsindicators[0]));
00311       switch (trailing) {
00312       case MINUTE:
00313         s << setw(2 + prec + (prec ? 1 : 0)) << setprecision(prec)
00314           << pieces[1] <<  char(tolower(dmsindicators[1]));
00315         break;
00316       case SECOND:
00317         s << setw(2) << pieces[1] <<  char(tolower(dmsindicators[1]))
00318           << setw(2 + prec + (prec ? 1 : 0)) << setprecision(prec)
00319           << pieces[2] <<  char(tolower(dmsindicators[2]));
00320         break;
00321       default:
00322         break;
00323       }
00324     }
00325     if (ind != NONE && ind != AZIMUTH)
00326       s << hemispheres[(ind == LATITUDE ? 0 : 2) + (sign < 0 ? 0 : 1)];
00327     return s.str();
00328   }
00329 
00330 } // namespace GeographicLib