00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
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
00198
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
00262 std::cout << "ERROR: " << e.what() << "\n";
00263 retval = 1;
00264 }
00265 }
00266 return retval;
00267 }