Planimeter.cpp

Go to the documentation of this file.
00001 /**
00002  * \file Planimeter.cpp
00003  * \brief Command line utility for measuring the area of geodesic polygons
00004  *
00005  * Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licensed
00006  * under the LGPL.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  *
00009  * Compile with -I../include and link with Geodesic.o GeodesicLine.o DMS.o
00010  *
00011  * See the <a href="Planimeter.1.html">man page</a> for usage
00012  * information.
00013  **********************************************************************/
00014 
00015 #include "GeographicLib/Geodesic.hpp"
00016 #include "GeographicLib/DMS.hpp"
00017 #include "GeographicLib/GeoCoords.hpp"
00018 #include <iostream>
00019 
00020 #include "Planimeter.usage"
00021 
00022 int main(int argc, char* argv[]) {
00023   using namespace GeographicLib;
00024   typedef Math::real real;
00025 
00026   class Accumulator {
00027     // Compute a sum at double precision.  This is Algorithm 4.1, of T. Ogita,
00028     // S. M. Rump, S. Oishi, Accurate sum and dot product, SIAM J. Sci. Comp.,
00029     // 26(6) 1955-1988 (2005).
00030   private:
00031     // _s accumulates for the straight sum
00032     // _t accumulates the errors.
00033     real _s, _t;
00034     // Error free transformation of a sum;
00035     // see Knuth, TAOCP, Vol 2, 4.2.2, Theorem B.
00036     static inline real sum(real u, real v, real& t) {
00037       volatile real s = u + v;
00038       volatile real up = s - v;
00039       volatile real vpp = s - up;
00040       up -= u;
00041       vpp -= v;
00042       t = -(up + vpp);
00043       // u + v =       s      + t
00044       //       = round(u + v) + t
00045       return s;
00046     }
00047   public:
00048     Accumulator() throw() : _s(0), _t(0) {};
00049     void Clear() throw() { _s = 0; _t = 0; }
00050     // Accumulate y
00051     void Add(real y) throw() {
00052       _s = sum(_s, y, y);
00053       _t += y;
00054     }
00055     void Negate() throw() { _s *= -1; _t *= -1; }
00056     // Return sum + y (don't accumulate)
00057     real Sum(real y) const throw() {
00058       real s = sum(_s, y, y);
00059       return s + (_t + y);
00060     }
00061     // Return sum
00062     real Sum() const throw() { return _s + _t; }
00063   };
00064 
00065   class GeodesicPolygon {
00066   private:
00067     const Geodesic& _g;
00068     const real _area0;          // Full ellipsoid area
00069     unsigned _num;
00070     int _crossings;
00071     Accumulator _area, _perimeter;
00072     real _lat0, _lon0, _lat1, _lon1;
00073     // Copied from Geodesic class
00074     static inline real AngNormalize(real x) throw() {
00075       // Place angle in [-180, 180).  Assumes x is in [-540, 540).
00076       return x >= 180 ? x - 360 : x < -180 ? x + 360 : x;
00077     }
00078     static inline int transit(real lon1, real lon2) {
00079       // Return 1 or -1 if crossing prime meridian in east or west direction.
00080       // Otherwise return zero.
00081       lon1 = AngNormalize(lon1);
00082       lon2 = AngNormalize(lon2);
00083       // treat lon12 = -180 as an eastward geodesic, so convert to 180.
00084       real lon12 = -AngNormalize(lon1 - lon2); // In (-180, 180]
00085       int cross =
00086         lon1 < 0 && lon2 >= 0 && lon12 > 0 ? 1 :
00087         lon2 < 0 && lon1 >= 0 && lon12 < 0 ? -1 : 0;
00088       return cross;
00089     }
00090   public:
00091     GeodesicPolygon(const Geodesic& g) throw()
00092       : _g(g)
00093       , _area0(_g.EllipsoidArea())
00094     {
00095       Clear();
00096     }
00097     void Clear() throw() {
00098       _num = 0;
00099       _crossings = 0;
00100       _area.Clear();
00101       _perimeter.Clear();
00102       _lat0 = _lon0 = _lat1 = _lon1 = 0;
00103     }
00104     void AddPoint(real lat, real lon) throw() {
00105       if (_num == 0) {
00106         _lat0 = _lat1 = lat;
00107         _lon0 = _lon1 = lon;
00108       } else {
00109         real s12, S12, t;
00110         _g.GenInverse(_lat1, _lon1, lat, lon,
00111                       Geodesic::DISTANCE | Geodesic::AREA,
00112                       s12, t, t, t, t, t, S12);
00113         _perimeter.Add(s12);
00114         _area.Add(S12);
00115         if (_area.Sum() > _area0/2)
00116           _area.Add(-_area0);
00117         else if (_area.Sum() <= -_area0/2)
00118           _area.Add(_area0);
00119         _crossings += transit(_lon1, lon);
00120         _lat1 = lat;
00121         _lon1 = lon;
00122       }
00123       ++_num;
00124     }
00125     unsigned Compute(bool reverse, bool sign,
00126                      real& perimeter, real& area) const throw() {
00127       real s12, S12, t;
00128       if (_num < 2) {
00129         perimeter = area = 0;
00130         return _num;
00131       }
00132       _g.GenInverse(_lat1, _lon1, _lat0, _lon0,
00133                     Geodesic::DISTANCE | Geodesic::AREA,
00134                     s12, t, t, t, t, t, S12);
00135       perimeter = _perimeter.Sum(s12);
00136       Accumulator area1(_area);
00137       area1.Add(S12);
00138       if (area1.Sum() > _area0/2)
00139         area1.Add(-_area0);
00140       else if (area1.Sum() <= -_area0/2)
00141         area1.Add(_area0);
00142       int crossings = _crossings + transit(_lon1, _lon0);
00143       if (crossings & 1) {
00144         if (area1.Sum() < 0)
00145           area1.Sum(_area0/2);
00146         else
00147           area1.Sum(-_area0/2);
00148       }
00149       // area is with the clockwise sense.  If !reverse convert to
00150       // counter-clockwise convention.
00151       if (!reverse)
00152         area1.Negate();
00153       // If sign put area in (-area0/2, area0/2], else put area in [0, area0)
00154       if (sign) {
00155         if (area1.Sum() > _area0/2)
00156           area1.Add(-_area0);
00157         else if (area1.Sum() <= -_area0/2)
00158           area1.Add(_area0);
00159       } else {
00160         if (area1.Sum() >= _area0)
00161           area1.Add(-_area0);
00162         else if (area1.Sum() < 0)
00163           area1.Add(_area0);
00164       }
00165       area = area1.Sum();
00166       return _num;
00167     }
00168   };
00169 
00170   real
00171     a = Constants::WGS84_a<real>(),
00172     r = Constants::WGS84_r<real>();
00173   bool reverse = false, sign = false;
00174   for (int m = 1; m < argc; ++m) {
00175     std::string arg(argv[m]);
00176     if (arg == "-r")
00177       reverse = !reverse;
00178     else if (arg == "-s")
00179       sign = !sign;
00180     else if (arg == "-e") {
00181       if (m + 2 >= argc) return usage(1, true);
00182       try {
00183         a = DMS::Decode(std::string(argv[m + 1]));
00184         r = DMS::Decode(std::string(argv[m + 2]));
00185       }
00186       catch (const std::exception& e) {
00187         std::cerr << "Error decoding arguments of -e: " << e.what() << "\n";
00188         return 1;
00189       }
00190       m += 2;
00191     } else if (arg == "--version") {
00192       std::cout
00193         << PROGRAM_NAME
00194         << ": $Id: Planimeter.cpp 6978 2011-02-21 22:42:11Z karney $\n"
00195         << "GeographicLib version " << GEOGRAPHICLIB_VERSION << "\n";
00196       return 0;
00197     } else
00198       return usage(!(arg == "-h" || arg == "--help"), arg != "--help");
00199   }
00200 
00201   const Geodesic geod(a, r);
00202   GeodesicPolygon poly(geod);
00203   GeoCoords p;
00204 
00205   std::string s;
00206   real perimeter, area;
00207   unsigned num;
00208   while (std::getline(std::cin, s)) {
00209     try {
00210       p.Reset(s);
00211       if (p.Latitude() != p.Latitude() || p.Longitude() != p.Longitude())
00212         throw GeographicErr("NAN");
00213     }
00214     catch (const GeographicErr&) {
00215       num = poly.Compute(reverse, sign, perimeter, area);
00216       if (num > 0)
00217         std::cout << num << " "
00218                   << DMS::Encode(perimeter, 8, DMS::NUMBER) << " "
00219                   << DMS::Encode(area, 4, DMS::NUMBER) << "\n";
00220       poly.Clear();
00221       continue;
00222     }
00223     poly.AddPoint(p.Latitude(), p.Longitude());
00224   }
00225   num = poly.Compute(reverse, sign, perimeter, area);
00226   if (num > 0)
00227         std::cout << num << " "
00228                   << DMS::Encode(perimeter, 8, DMS::NUMBER) << " "
00229                   << DMS::Encode(area, 4, DMS::NUMBER) << "\n";
00230   poly.Clear();
00231   return 0;
00232 }