GeographicLib  1.21
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-2011) <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/DMS.hpp>
00011 #include <algorithm>
00012 #include <GeographicLib/Utility.hpp>
00013 
00014 #define GEOGRAPHICLIB_DMS_CPP "$Id: db38ddc05f7c27732da3aa820191a51200ce92ac $"
00015 
00016 RCSID_DECL(GEOGRAPHICLIB_DMS_CPP)
00017 RCSID_DECL(GEOGRAPHICLIB_DMS_HPP)
00018 RCSID_DECL(GEOGRAPHICLIB_CONSTANTS_HPP)
00019 RCSID_DECL(GEOGRAPHICLIB_MATH_HPP)
00020 
00021 namespace GeographicLib {
00022 
00023   using namespace std;
00024 
00025   const string DMS::hemispheres_ = "SNWE";
00026   const string DMS::signs_ = "-+";
00027   const string DMS::digits_ = "0123456789";
00028   const string DMS::dmsindicators_ = "D'\":";
00029   const string DMS::components_[] = {"degrees", "minutes", "seconds"};
00030 
00031   Math::real DMS::Decode(const std::string& dms, flag& ind) {
00032     string errormsg;
00033     string dmsa = dms;
00034     replace(dmsa, "\xc2\xb0", 'd'); // degree symbol (U+00b0 = UTF-8 c2 b0)
00035     replace(dmsa, "\xc2\xba", 'd'); // alt symbol (U+00ba = UTF-8 c2 ba)
00036     replace(dmsa, "\xe2\x81\xb0", 'd');  // sup zero (U+2070 = UTF-8 e2 81 b0)
00037     replace(dmsa, "\xe2\x80\xb2", '\''); // prime (U+2032 = UTF-8 e2 80 b2)
00038     replace(dmsa, "\xc2\xb4", '\'');     // acute accent (U+00b4 = UTF-8 c2 b4)
00039     replace(dmsa, "\xe2\x80\xb3", '"');  // dbl prime (U+2033 = UTF-8 e2 80 b3)
00040     replace(dmsa, "\xb0", 'd');          // bare degree symbol (b0)
00041     replace(dmsa, "\xba", 'd');          // bare alt symbol (ba)
00042     replace(dmsa, "\xb4", 'd');          // bare acute accent (b4)
00043     replace(dmsa, "''", '"');            // '' -> "
00044     do {                       // Executed once (provides the ability to break)
00045       int sign = 1;
00046       unsigned
00047         beg = 0,
00048         end = unsigned(dmsa.size());
00049       while (beg < end && isspace(dmsa[beg]))
00050         ++beg;
00051       while (beg < end && isspace(dmsa[end - 1]))
00052         --end;
00053       flag ind1 = NONE;
00054       int k = -1;
00055       if (end > beg && (k = Utility::lookup(hemispheres_, dmsa[beg])) >= 0) {
00056         ind1 = (k / 2) ? LONGITUDE : LATITUDE;
00057         sign = k % 2 ? 1 : -1;
00058         ++beg;
00059       }
00060       if (end > beg && (k = Utility::lookup(hemispheres_, dmsa[end-1])) >= 0) {
00061         if (k >= 0) {
00062           if (ind1 != NONE) {
00063             if (toupper(dmsa[beg - 1]) == toupper(dmsa[end - 1]))
00064               errormsg = "Repeated hemisphere indicators "
00065                 + Utility::str(dmsa[beg - 1])
00066                 + " in " + dmsa.substr(beg - 1, end - beg + 1);
00067             else
00068               errormsg = "Contradictory hemisphere indicators "
00069                 + Utility::str(dmsa[beg - 1]) + " and "
00070                 + Utility::str(dmsa[end - 1]) + " in "
00071                 + dmsa.substr(beg - 1, end - beg + 1);
00072             break;
00073           }
00074           ind1 = (k / 2) ? LONGITUDE : LATITUDE;
00075           sign = k % 2 ? 1 : -1;
00076           --end;
00077         }
00078       }
00079       if (end > beg && (k = Utility::lookup(signs_, dmsa[beg])) >= 0) {
00080         if (k >= 0) {
00081           sign *= k ? 1 : -1;
00082           ++beg;
00083         }
00084       }
00085       if (end == beg) {
00086         errormsg = "Empty or incomplete DMS string " + dmsa;
00087         break;
00088       }
00089       real ipieces[] = {0, 0, 0};
00090       real fpieces[] = {0, 0, 0};
00091       unsigned npiece = 0;
00092       real icurrent = 0;
00093       real fcurrent = 0;
00094       unsigned ncurrent = 0, p = beg;
00095       bool pointseen = false;
00096       unsigned digcount = 0;
00097       while (p < end) {
00098         char x = dmsa[p++];
00099         if ((k = Utility::lookup(digits_, x)) >= 0) {
00100           ++ncurrent;
00101           if (digcount > 0)
00102             ++digcount;         // Count of decimal digits
00103           else
00104             icurrent = 10 * icurrent + k;
00105         } else if (x == '.') {
00106           if (pointseen) {
00107             errormsg = "Multiple decimal points in "
00108               + dmsa.substr(beg, end - beg);
00109             break;
00110           }
00111           pointseen = true;
00112           digcount = 1;
00113         } else if ((k = Utility::lookup(dmsindicators_, x)) >= 0) {
00114           if (k >= 3) {
00115             if (p == end) {
00116               errormsg = "Illegal for : to appear at the end of " +
00117                 dmsa.substr(beg, end - beg);
00118               break;
00119             }
00120             k = npiece;
00121           }
00122           if (unsigned(k) == npiece - 1) {
00123             errormsg = "Repeated " + components_[k] +
00124               " component in " + dmsa.substr(beg, end - beg);
00125             break;
00126           } else if (unsigned(k) < npiece) {
00127             errormsg = components_[k] + " component follows "
00128               + components_[npiece - 1] + " component in "
00129               + dmsa.substr(beg, end - beg);
00130             break;
00131           }
00132           if (ncurrent == 0) {
00133             errormsg = "Missing numbers in " + components_[k] +
00134               " component of " + dmsa.substr(beg, end - beg);
00135             break;
00136           }
00137           if (digcount > 1) {
00138             istringstream s(dmsa.substr(p - digcount - 1, digcount));
00139             s >> fcurrent;
00140           }
00141           ipieces[k] = icurrent;
00142           fpieces[k] = icurrent + fcurrent;
00143           if (p < end) {
00144             npiece = k + 1;
00145             icurrent = fcurrent = 0;
00146             ncurrent = digcount = 0;
00147           }
00148         } else if (Utility::lookup(signs_, x) >= 0) {
00149           errormsg = "Internal sign in DMS string "
00150             + dmsa.substr(beg, end - beg);
00151           break;
00152         } else {
00153           errormsg = "Illegal character " + Utility::str(x) + " in DMS string "
00154             + dmsa.substr(beg, end - beg);
00155           break;
00156         }
00157       }
00158       if (!errormsg.empty())
00159         break;
00160       if (Utility::lookup(dmsindicators_, dmsa[p - 1]) < 0) {
00161         if (npiece >= 3) {
00162           errormsg = "Extra text following seconds in DMS string "
00163             + dmsa.substr(beg, end - beg);
00164           break;
00165         }
00166         if (ncurrent == 0) {
00167           errormsg = "Missing numbers in trailing component of "
00168             + dmsa.substr(beg, end - beg);
00169           break;
00170         }
00171         if (digcount > 1) {
00172           istringstream s(dmsa.substr(p - digcount, digcount));
00173           s >> fcurrent;
00174         }
00175         ipieces[npiece] = icurrent;
00176         fpieces[npiece] = icurrent + fcurrent;
00177       }
00178       if (pointseen && digcount == 0) {
00179         errormsg = "Decimal point in non-terminal component of "
00180           + dmsa.substr(beg, end - beg);
00181         break;
00182       }
00183       // Note that we accept 59.999999... even though it rounds to 60.
00184       if (ipieces[1] >= 60) {
00185         errormsg = "Minutes " + Utility::str(fpieces[1])
00186           + " not in range [0, 60)";
00187         break;
00188       }
00189       if (ipieces[2] >= 60) {
00190         errormsg = "Seconds " + Utility::str(fpieces[2])
00191           + " not in range [0, 60)";
00192         break;
00193       }
00194       ind = ind1;
00195       // Assume check on range of result is made by calling routine (which
00196       // might be able to offer a better diagnostic).
00197       return real(sign) * (fpieces[0] + (fpieces[1] + fpieces[2] / 60) / 60);
00198     } while (false);
00199     real val = Utility::nummatch<real>(dmsa);
00200     if (val == 0)
00201       throw GeographicErr(errormsg);
00202     else
00203       ind = NONE;
00204     return val;
00205   }
00206 
00207   void DMS::DecodeLatLon(const std::string& stra, const std::string& strb,
00208                          real& lat, real& lon, bool swaplatlong) {
00209     real a, b;
00210     flag ia, ib;
00211     a = Decode(stra, ia);
00212     b = Decode(strb, ib);
00213     if (ia == NONE && ib == NONE) {
00214       // Default to lat, long unless swaplatlong
00215       ia = swaplatlong ? LONGITUDE : LATITUDE;
00216       ib = swaplatlong ? LATITUDE : LONGITUDE;
00217     } else if (ia == NONE)
00218       ia = flag(LATITUDE + LONGITUDE - ib);
00219     else if (ib == NONE)
00220       ib = flag(LATITUDE + LONGITUDE - ia);
00221     if (ia == ib)
00222       throw GeographicErr("Both " + stra + " and "
00223                           + strb + " interpreted as "
00224                           + (ia == LATITUDE ? "latitudes" : "longitudes"));
00225     real
00226       lat1 = ia == LATITUDE ? a : b,
00227       lon1 = ia == LATITUDE ? b : a;
00228     if (lat1 < -90 || lat1 > 90)
00229       throw GeographicErr("Latitude " + Utility::str(lat1)
00230                           + "d not in [-90d, 90d]");
00231     if (lon1 < -180 || lon1 > 360)
00232       throw GeographicErr("Longitude " + Utility::str(lon1)
00233                           + "d not in [-180d, 360d]");
00234     if (lon1 >= 180)
00235       lon1 -= 360;
00236     lat = lat1;
00237     lon = lon1;
00238   }
00239 
00240   Math::real DMS::DecodeAngle(const std::string& angstr) {
00241     flag ind;
00242     real ang = Decode(angstr, ind);
00243     if (ind != NONE)
00244       throw GeographicErr("Arc angle " + angstr
00245                           + " includes a hemisphere, N/E/W/S");
00246     return ang;
00247   }
00248 
00249   Math::real DMS::DecodeAzimuth(const std::string& azistr) {
00250     flag ind;
00251     real azi = Decode(azistr, ind);
00252     if (ind == LATITUDE)
00253       throw GeographicErr("Azimuth " + azistr
00254                           + " has a latitude hemisphere, N/S");
00255     if (azi < -180 || azi > 360)
00256       throw GeographicErr("Azimuth " + azistr + " not in range [-180d, 360d]");
00257     if (azi >= 180) azi -= 360;
00258     return azi;
00259   }
00260 
00261   string DMS::Encode(real angle, component trailing, unsigned prec, flag ind,
00262                      char dmssep) {
00263     // Assume check on range of input angle has been made by calling
00264     // routine (which might be able to offer a better diagnostic).
00265     if (!Math::isfinite(angle))
00266       return angle < 0 ? string("-inf") :
00267         (angle > 0 ? string("inf") : string("nan"));
00268 
00269     // 15 - 2 * trailing = ceiling(log10(2^53/90/60^trailing)).
00270     // This suffices to give full real precision for numbers in [-90,90]
00271     prec = min(15 - 2 * unsigned(trailing), prec);
00272     real scale = 1;
00273     for (unsigned i = 0; i < unsigned(trailing); ++i)
00274       scale *= 60;
00275     for (unsigned i = 0; i < prec; ++i)
00276       scale *= 10;
00277     if (ind == AZIMUTH)
00278       angle -= floor(angle/360) * 360;
00279     int sign = angle < 0 ? -1 : 1;
00280     angle *= sign;
00281 
00282     // Break off integer part to preserve precision in manipulation of
00283     // fractional part.
00284     real
00285       idegree = floor(angle),
00286       fdegree = floor((angle - idegree) * scale + real(0.5)) / scale;
00287     if (fdegree >= 1) {
00288       idegree += 1;
00289       fdegree -= 1;
00290     }
00291     real pieces[3] = {fdegree, 0, 0};
00292     for (unsigned i = 1; i <= unsigned(trailing); ++i) {
00293       real
00294         ip = floor(pieces[i - 1]),
00295         fp = pieces[i - 1] - ip;
00296       pieces[i] = fp * 60;
00297       pieces[i - 1] = ip;
00298     }
00299     pieces[0] += idegree;
00300     ostringstream s;
00301     s << fixed << setfill('0');
00302     if (ind == NONE && sign < 0)
00303       s << '-';
00304     switch (trailing) {
00305     case DEGREE:
00306       if (ind != NONE)
00307         s << setw(1 + min(int(ind), 2) + prec + (prec ? 1 : 0));
00308       s << setprecision(prec) << pieces[0];
00309       // Don't include degree designator (d) if it is the trailing component.
00310       break;
00311     default:
00312       if (ind != NONE)
00313         s << setw(1 + min(int(ind), 2));
00314       s << setprecision(0) << pieces[0]
00315         << (dmssep ? dmssep : char(tolower(dmsindicators_[0])));
00316       switch (trailing) {
00317       case MINUTE:
00318         s << setw(2 + prec + (prec ? 1 : 0)) << setprecision(prec) << pieces[1];
00319         if (!dmssep)
00320           s << char(tolower(dmsindicators_[1]));
00321         break;
00322       case SECOND:
00323         s << setw(2)
00324           << pieces[1] << (dmssep ? dmssep : char(tolower(dmsindicators_[1])))
00325           << setw(2 + prec + (prec ? 1 : 0)) << setprecision(prec) << pieces[2];
00326         if (!dmssep)
00327           s << char(tolower(dmsindicators_[2]));
00328         break;
00329       default:
00330         break;
00331       }
00332     }
00333     if (ind != NONE && ind != AZIMUTH)
00334       s << hemispheres_[(ind == LATITUDE ? 0 : 2) + (sign < 0 ? 0 : 1)];
00335     return s.str();
00336   }
00337 
00338   string DMS::Encode(real angle, component trailing, unsigned prec, flag ind)
00339   { return Encode(angle, trailing, prec, ind, char(0)); }
00340 
00341 } // namespace GeographicLib