Geod.cpp

Go to the documentation of this file.
00001 /**
00002  * \file Geod.cpp
00003  * \brief Command line utility for geodesic calculations
00004  *
00005  * Copyright (c) Charles Karney (2009, 2010) <charles@karney.com>
00006  * and licensed under the LGPL.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  *
00009  * Compile with -I../include and link with Geodesic.o DMS.o
00010  *
00011  * See \ref geod for usage information.
00012  **********************************************************************/
00013 
00014 #include "GeographicLib/Geodesic.hpp"
00015 #include "GeographicLib/DMS.hpp"
00016 #include <iostream>
00017 #include <iomanip>
00018 #include <sstream>
00019 
00020 int usage(int retval) {
00021   ( retval ? std::cerr : std::cout ) <<
00022 "Usage: Geod [-l lat1 lon1 azi1 | -i] [-a] [-n | -e a r]\n\
00023             [-d] [-b] [-f] [-p prec] [-h]\n\
00024 $Id: Geod.cpp 6827 2010-05-20 19:56:18Z karney $\n\
00025 \n\
00026 Perform geodesic calculations.\n\
00027 \n\
00028 The shortest path between two points on the ellipsoid at (lat1, lon1) and\n\
00029 (lat2, lon2) is called the geodesic.  Its length is s12 and the geodesic\n\
00030 from point 1 to point 2 has azimuths azi1 and azi2 at the two end\n\
00031 points.  The reduced length of the geodesic, m12, is defined such that\n\
00032 if the initial azimuth is perturbed by dazi1 (radians) then the second\n\
00033 point is displaced by m12*dazi1 in the direction perpendicular to the\n\
00034 geodesic.  On a flat surface, we have m12 = s12.\n\
00035 \n\
00036 Geod operates in one of three modes:\n\
00037 \n\
00038 (1) It accepts lines on the standard input containing \"lat1 lon1 azi1\n\
00039     s12\" and prints \"lat2 lon2 azi2 m12\" on standard output.  This is\n\
00040     the direct geodesic calculation.\n\
00041 \n\
00042 (2) Command line arguments \"-l lat1 lon1 azi1\" specify a geodesic line.\n\
00043     Geod then accepts a sequence of s12 values (one per line) on\n\
00044     standard input and prints \"lat2 lon2 azi2 m12\" for each.  This\n\
00045     generates a sequence of points on a single geodesic.\n\
00046 \n\
00047 (3) With the -i command line argument, Geod performs the inverse\n\
00048     geodesic calculation.  It reads lines containing \"lat1 lon1 lat2\n\
00049     lon2\" and prints the corresponding values of \"azi1 azi2 s12 m12\".\n\
00050 \n\
00051 By default, the WGS84 ellipsoid is used.  Specifying \"-e a r\" sets the\n\
00052 equatorial radius of the ellipsoid to \"a\" and the reciprocal flattening\n\
00053 to r.  Setting r = 0 results in a sphere.  Specify r < 0 for a prolate\n\
00054 ellipsoid.  The -n option uses the international ellipsoid (equivalent to\n\
00055 \"-e 6378388 297\").\n\
00056 \n\
00057 Output of angles is as decimal degrees.  If -d is specified the output\n\
00058 is as degrees, minutes, seconds.  Input can be in either style.  d, ',\n\
00059 and \" are used to denote degrees, minutes, and seconds, with the least\n\
00060 significant designator optional.  By default, latitude precedes\n\
00061 longitude for each point; however on input either may be given first by\n\
00062 appending N or S to the latitude and E or W to the longitude.  Azimuths\n\
00063 (measured clockwise from north) give the heading of the geodesic.  The\n\
00064 azimuth azi2 is the forward azimuth (the heading beyond point 2).  If\n\
00065 the -b flag is given, azi2 is converted to a back azimuth (the direction\n\
00066 back to point 1) for output.\n\
00067 \n\
00068 s12 is given in meters, unless the -a flag is given.  In that case, s12\n\
00069 (on both input and output) are given as the arc length on the auxiliary\n\
00070 sphere a12 (measured in degrees).  In these terms, 180 degrees is the\n\
00071 distance from one equator crossing to the next or from the minimum\n\
00072 latitude to the maximum latitude.  Distances greater than 180 degrees do\n\
00073 not correspond to shortest paths.  m12 is always given in meters.\n\
00074 \n\
00075 The output lines consist of the four quantities needed to complete the\n\
00076 specification of the geodesic.  With the -f option, each line of output\n\
00077 is a complete geodesic specification consisting of nine quantities\n\
00078 \n\
00079     lat1 lon1 azi1 lat2 lon2 azi2 s12 a12 m12\n\
00080 \n\
00081 where here s12 is the distance and a12 the arc length.\n\
00082 \n\
00083 -p prec (default 3) gives the precision of the output relative to 1m.\n\
00084 The minimum value of prec is 0 (1 m accuracy) and the maximum value is\n\
00085 10 (0.1 nm accuracy, but then the last digits are unreliable).\n\
00086 \n\
00087 -h prints this help.\n";
00088   return retval;
00089 }
00090 
00091 typedef GeographicLib::Math::real real;
00092 
00093 std::string LatLonString(real lat, real lon, int prec, bool dms) {
00094   using namespace GeographicLib;
00095   return dms ?
00096     DMS::Encode(lat, prec + 5, DMS::LATITUDE) + " " +
00097     DMS::Encode(lon, prec + 5, DMS::LONGITUDE) :
00098     DMS::Encode(lat, prec + 5, DMS::NUMBER) + " " +
00099     DMS::Encode(lon, prec + 5, DMS::NUMBER);
00100 }
00101 
00102 std::string AzimuthString(real azi, int prec, bool dms) {
00103   using namespace GeographicLib;
00104   return dms ? DMS::Encode(azi, prec + 5, DMS::AZIMUTH) :
00105     DMS::Encode(azi >= 180 ? azi - 360 : azi, prec + 5, DMS::NUMBER);
00106 }
00107 
00108 std::string DistanceStrings(real s12, real a12,
00109                             bool full, bool arcmode, int prec, bool dms) {
00110   using namespace GeographicLib;
00111   std::string s;
00112   if (full || !arcmode)
00113     s += DMS::Encode(s12, prec, DMS::NUMBER);
00114   if (full)
00115     s += " ";
00116   if (full || arcmode)
00117     s += dms ? DMS::Encode(a12, prec + 5, DMS::NONE) :
00118       DMS::Encode(a12, prec + 5, DMS::NUMBER);
00119   return s;
00120 }
00121 
00122 real ReadDistance(const std::string& s, bool arcmode) {
00123   using namespace GeographicLib;
00124   return arcmode ? DMS::DecodeAngle(s) : DMS::Decode(s);
00125 }
00126 
00127 int main(int argc, char* argv[]) {
00128   using namespace GeographicLib;
00129   bool linecalc = false, inverse = false, arcmode = false,
00130     dms = false, full = false;
00131   real
00132     a = Constants::WGS84_a(),
00133     r = Constants::WGS84_r();
00134   real lat1, lon1, azi1, lat2, lon2, azi2, s12, m12, a12;
00135   real azi2sense = 0;
00136   int prec = 3;
00137 
00138   for (int m = 1; m < argc; ++m) {
00139     std::string arg(argv[m]);
00140     if (arg == "-i") {
00141       inverse = true;
00142       linecalc = false;
00143     } else if (arg == "-a")
00144       arcmode = true;
00145     else if (arg == "-l") {
00146       inverse = false;
00147       linecalc = true;
00148       if (m + 3 >= argc) return usage(1);
00149       try {
00150         DMS::DecodeLatLon(std::string(argv[m + 1]), std::string(argv[m + 2]),
00151                           lat1, lon1);
00152         azi1 = DMS::DecodeAzimuth(std::string(argv[m + 3]));
00153       }
00154       catch (const std::exception& e) {
00155         std::cerr << "Error decoding arguments of -l: " << e.what() << "\n";
00156         return 1;
00157       }
00158       m += 3;
00159     } else if (arg == "-n") {
00160       a = 6378388;
00161       r = 297;
00162     } else if (arg == "-e") {
00163       if (m + 2 >= argc) return usage(1);
00164       try {
00165         a = DMS::Decode(std::string(argv[m + 1]));
00166         r = DMS::Decode(std::string(argv[m + 2]));
00167       }
00168       catch (const std::exception& e) {
00169         std::cerr << "Error decoding arguments of -e: " << e.what() << "\n";
00170         return 1;
00171       }
00172       m += 2;
00173     }
00174     else if (arg == "-d")
00175       dms = true;
00176     else if (arg == "-b")
00177       azi2sense = 180;
00178     else if (arg == "-f")
00179       full = true;
00180     else if (arg == "-p") {
00181       if (++m == argc) return usage(1);
00182       std::istringstream str(argv[m]);
00183       char c;
00184       if (!(str >> prec) || (str >> c)) {
00185           std::cerr << "Precision " << argv[m] << " is not a number\n";
00186           return 1;
00187       }
00188     } else
00189       return usage(arg != "-h");
00190   }
00191 
00192   const Geodesic geod(a, r);
00193   GeodesicLine l;
00194   if (linecalc)
00195     l = geod.Line(lat1, lon1, azi1);
00196 
00197   // Max precision = 10: 0.1 nm in distance, 10^-15 deg (= 0.11 nm),
00198   // 10^-11 sec (= 0.3 nm).
00199   prec = std::min(10, std::max(0, prec));
00200   std::cout << std::fixed << std::setprecision(prec);
00201   std::string s;
00202   int retval = 0;
00203   while (std::getline(std::cin, s)) {
00204     try {
00205       std::istringstream str(s);
00206       if (inverse) {
00207         std::string slat1, slon1, slat2, slon2;
00208         if (!(str >> slat1 >> slon1 >> slat2 >> slon2))
00209           throw GeographicErr("Incomplete input: " + s);
00210         std::string strc;
00211         if (str >> strc)
00212           throw GeographicErr("Extraneous input: " + strc);
00213         DMS::DecodeLatLon(slat1, slon1, lat1, lon1);
00214         DMS::DecodeLatLon(slat2, slon2, lat2, lon2);
00215         a12 = geod.Inverse(lat1, lon1, lat2, lon2, s12, azi1, azi2, m12);
00216         if (full)
00217           std::cout << LatLonString(lat1, lon1, prec, dms) << " ";
00218         std::cout << AzimuthString(azi1, prec, dms) << " ";
00219         if (full)
00220           std::cout << LatLonString(lat2, lon2, prec, dms) << " ";
00221         std::cout << AzimuthString(azi2 + azi2sense, prec, dms) << " "
00222                   << DistanceStrings(s12, a12, full, arcmode, prec, dms) << " "
00223                   << DMS::Encode(m12, prec, DMS::NUMBER) << "\n";
00224       } else {
00225         if (linecalc) {
00226           std::string ss12;
00227           if (!(str >> ss12))
00228             throw GeographicErr("Incomplete input: " + s);
00229           std::string strc;
00230           if (str >> strc)
00231             throw GeographicErr("Extraneous input: " + strc);
00232           s12 = ReadDistance(ss12, arcmode);
00233           a12 = l.Position(s12, lat2, lon2, azi2, m12, arcmode);
00234         } else {
00235           std::string slat1, slon1, sazi1, ss12;
00236           if (!(str >> slat1 >> slon1 >> sazi1 >> ss12))
00237             throw GeographicErr("Incomplete input: " + s);
00238           std::string strc;
00239           if (str >> strc)
00240             throw GeographicErr("Extraneous input: " + strc);
00241           DMS::DecodeLatLon(slat1, slon1, lat1, lon1);
00242           azi1 = DMS::DecodeAzimuth(sazi1);
00243           s12 = ReadDistance(ss12, arcmode);
00244           a12 =
00245             geod.Direct(lat1, lon1, azi1, s12, lat2, lon2, azi2, m12, arcmode);
00246         }
00247         if (arcmode)
00248           std::swap(s12, a12);
00249         if (full)
00250           std::cout << LatLonString(lat1, lon1, prec, dms) << " "
00251                     << AzimuthString(azi1, prec, dms) << " ";
00252         std::cout << LatLonString(lat2, lon2, prec, dms) << " "
00253                   << AzimuthString(azi2 + azi2sense, prec, dms);
00254         if (full)
00255           std::cout << " "
00256                     << DistanceStrings(s12, a12, full, arcmode, prec, dms);
00257         std::cout << " " << DMS::Encode(m12, prec, DMS::NUMBER) << "\n";
00258       }
00259     }
00260     catch (const std::exception& e) {
00261       // Write error message cout so output lines match input lines
00262       std::cout << "ERROR: " << e.what() << "\n";
00263       retval = 1;
00264     }
00265   }
00266   return retval;
00267 }

Generated on 21 May 2010 for GeographicLib by  doxygen 1.6.1