GeodSolve.cpp
Go to the documentation of this file.
1 
12 #include <iostream>
13 #include <string>
14 #include <sstream>
15 #include <fstream>
20 #include <GeographicLib/DMS.hpp>
22 
23 #if defined(_MSC_VER)
24 // Squelch warnings about constant conditional expressions and potentially
25 // uninitialized local variables
26 # pragma warning (disable: 4127 4701)
27 #endif
28 
29 #include "GeodSolve.usage"
30 
32 
33 std::string LatLonString(real lat, real lon, int prec, bool dms, char dmssep,
34  bool longfirst) {
35  using namespace GeographicLib;
36  std::string
37  latstr = dms ? DMS::Encode(lat, prec + 5, DMS::LATITUDE, dmssep) :
38  DMS::Encode(lat, prec + 5, DMS::NUMBER),
39  lonstr = dms ? DMS::Encode(lon, prec + 5, DMS::LONGITUDE, dmssep) :
40  DMS::Encode(lon, prec + 5, DMS::NUMBER);
41  return
42  (longfirst ? lonstr : latstr) + " " + (longfirst ? latstr : lonstr);
43 }
44 
45 std::string AzimuthString(real azi, int prec, bool dms, char dmssep) {
46  using namespace GeographicLib;
47  return dms ? DMS::Encode(azi, prec + 5, DMS::AZIMUTH, dmssep) :
48  DMS::Encode(azi, prec + 5, DMS::NUMBER);
49 }
50 
51 std::string DistanceStrings(real s12, real a12,
52  bool full, bool arcmode, int prec, bool dms) {
53  using namespace GeographicLib;
54  std::string s;
55  if (full || !arcmode)
56  s += Utility::str(s12, prec);
57  if (full)
58  s += " ";
59  if (full || arcmode)
60  s += DMS::Encode(a12, prec + 5, dms ? DMS::NONE : DMS::NUMBER);
61  return s;
62 }
63 
64 real ReadDistance(const std::string& s, bool arcmode) {
65  using namespace GeographicLib;
66  return arcmode ? DMS::DecodeAngle(s) : Utility::val<real>(s);
67 }
68 
69 int main(int argc, const char* const argv[]) {
70  try {
71  using namespace GeographicLib;
72  enum { NONE = 0, LINE, DIRECT, INVERSE };
74  bool inverse = false, arcmode = false,
75  dms = false, full = false, exact = false, unroll = false,
76  longfirst = false, azi2back = false, fraction = false,
77  arcmodeline = false;
78  real
81  real lat1, lon1, azi1, lat2, lon2, azi2, s12, m12, a12, M12, M21, S12,
82  mult = 1;
83  int linecalc = NONE, prec = 3;
84  std::string istring, ifile, ofile, cdelim;
85  char lsep = ';', dmssep = char(0);
86 
87  for (int m = 1; m < argc; ++m) {
88  std::string arg(argv[m]);
89  if (arg == "-i") {
90  inverse = true;
91  linecalc = NONE;
92  } else if (arg == "-a")
93  arcmode = !arcmode;
94  else if (arg == "-F")
95  fraction = true;
96  else if (arg == "-L" || arg == "-l") { // -l is DEPRECATED
97  inverse = false;
98  linecalc = LINE;
99  if (m + 3 >= argc) return usage(1, true);
100  try {
101  DMS::DecodeLatLon(std::string(argv[m + 1]), std::string(argv[m + 2]),
102  lat1, lon1, longfirst);
103  azi1 = DMS::DecodeAzimuth(std::string(argv[m + 3]));
104  }
105  catch (const std::exception& e) {
106  std::cerr << "Error decoding arguments of -L: " << e.what() << "\n";
107  return 1;
108  }
109  m += 3;
110  } else if (arg == "-D") {
111  inverse = false;
112  linecalc = DIRECT;
113  if (m + 4 >= argc) return usage(1, true);
114  try {
115  DMS::DecodeLatLon(std::string(argv[m + 1]), std::string(argv[m + 2]),
116  lat1, lon1, longfirst);
117  azi1 = DMS::DecodeAzimuth(std::string(argv[m + 3]));
118  s12 = ReadDistance(std::string(argv[m + 4]), arcmode);
119  arcmodeline = arcmode;
120  }
121  catch (const std::exception& e) {
122  std::cerr << "Error decoding arguments of -D: " << e.what() << "\n";
123  return 1;
124  }
125  m += 4;
126  } else if (arg == "-I") {
127  inverse = false;
128  linecalc = INVERSE;
129  if (m + 4 >= argc) return usage(1, true);
130  try {
131  DMS::DecodeLatLon(std::string(argv[m + 1]), std::string(argv[m + 2]),
132  lat1, lon1, longfirst);
133  DMS::DecodeLatLon(std::string(argv[m + 3]), std::string(argv[m + 4]),
134  lat2, lon2, longfirst);
135  }
136  catch (const std::exception& e) {
137  std::cerr << "Error decoding arguments of -I: " << e.what() << "\n";
138  return 1;
139  }
140  m += 4;
141  } else if (arg == "-e") {
142  if (m + 2 >= argc) return usage(1, true);
143  try {
144  a = Utility::val<real>(std::string(argv[m + 1]));
145  f = Utility::fract<real>(std::string(argv[m + 2]));
146  }
147  catch (const std::exception& e) {
148  std::cerr << "Error decoding arguments of -e: " << e.what() << "\n";
149  return 1;
150  }
151  m += 2;
152  } else if (arg == "-u")
153  unroll = true;
154  else if (arg == "-d") {
155  dms = true;
156  dmssep = '\0';
157  } else if (arg == "-:") {
158  dms = true;
159  dmssep = ':';
160  } else if (arg == "-w")
161  longfirst = !longfirst;
162  else if (arg == "-b")
163  azi2back = true;
164  else if (arg == "-f")
165  full = true;
166  else if (arg == "-p") {
167  if (++m == argc) return usage(1, true);
168  try {
169  prec = Utility::val<int>(std::string(argv[m]));
170  }
171  catch (const std::exception&) {
172  std::cerr << "Precision " << argv[m] << " is not a number\n";
173  return 1;
174  }
175  } else if (arg == "-E")
176  exact = true;
177  else if (arg == "--input-string") {
178  if (++m == argc) return usage(1, true);
179  istring = argv[m];
180  } else if (arg == "--input-file") {
181  if (++m == argc) return usage(1, true);
182  ifile = argv[m];
183  } else if (arg == "--output-file") {
184  if (++m == argc) return usage(1, true);
185  ofile = argv[m];
186  } else if (arg == "--line-separator") {
187  if (++m == argc) return usage(1, true);
188  if (std::string(argv[m]).size() != 1) {
189  std::cerr << "Line separator must be a single character\n";
190  return 1;
191  }
192  lsep = argv[m][0];
193  } else if (arg == "--comment-delimiter") {
194  if (++m == argc) return usage(1, true);
195  cdelim = argv[m];
196  } else if (arg == "--version") {
197  std::cout << argv[0] << ": GeographicLib version "
198  << GEOGRAPHICLIB_VERSION_STRING << "\n";
199  return 0;
200  } else
201  return usage(!(arg == "-h" || arg == "--help"), arg != "--help");
202  }
203 
204  if (!ifile.empty() && !istring.empty()) {
205  std::cerr << "Cannot specify --input-string and --input-file together\n";
206  return 1;
207  }
208  if (ifile == "-") ifile.clear();
209  std::ifstream infile;
210  std::istringstream instring;
211  if (!ifile.empty()) {
212  infile.open(ifile.c_str());
213  if (!infile.is_open()) {
214  std::cerr << "Cannot open " << ifile << " for reading\n";
215  return 1;
216  }
217  } else if (!istring.empty()) {
218  std::string::size_type m = 0;
219  while (true) {
220  m = istring.find(lsep, m);
221  if (m == std::string::npos)
222  break;
223  istring[m] = '\n';
224  }
225  instring.str(istring);
226  }
227  std::istream* input = !ifile.empty() ? &infile :
228  (!istring.empty() ? &instring : &std::cin);
229 
230  std::ofstream outfile;
231  if (ofile == "-") ofile.clear();
232  if (!ofile.empty()) {
233  outfile.open(ofile.c_str());
234  if (!outfile.is_open()) {
235  std::cerr << "Cannot open " << ofile << " for writing\n";
236  return 1;
237  }
238  }
239  std::ostream* output = !ofile.empty() ? &outfile : &std::cout;
240 
241  // GeodesicExact mask values are the same as Geodesic
242  unsigned outmask = Geodesic::LATITUDE | Geodesic::LONGITUDE |
243  Geodesic::AZIMUTH; // basic output quantities
244  outmask |= inverse ? Geodesic::DISTANCE : // distance-related flags
246  // longitude unrolling
247  outmask |= unroll ? Geodesic::LONG_UNROLL : Geodesic::NONE;
248  // full output -- don't use Geodesic::ALL since this includes DISTANCE_IN
249  outmask |= full ? (Geodesic::DISTANCE | Geodesic::REDUCEDLENGTH |
252 
253  const Geodesic geods(a, f);
254  const GeodesicExact geode(a, f);
255  GeodesicLine ls;
257  if (linecalc) {
258  if (linecalc == LINE) fraction = false;
259  if (exact) {
260  le = linecalc == DIRECT ?
261  geode.GenDirectLine(lat1, lon1, azi1, arcmodeline, s12, outmask) :
262  linecalc == INVERSE ?
263  geode.InverseLine(lat1, lon1, lat2, lon2, outmask) :
264  // linecalc == LINE
265  geode.Line(lat1, lon1, azi1, outmask);
266  mult = fraction ? le.GenDistance(arcmode) : 1;
267  if (linecalc == INVERSE) azi1 = le.Azimuth();
268  } else {
269  ls = linecalc == DIRECT ?
270  geods.GenDirectLine(lat1, lon1, azi1, arcmodeline, s12, outmask) :
271  linecalc == INVERSE ?
272  geods.InverseLine(lat1, lon1, lat2, lon2, outmask) :
273  // linecalc == LINE
274  geods.Line(lat1, lon1, azi1, outmask);
275  mult = fraction ? ls.GenDistance(arcmode) : 1;
276  if (linecalc == INVERSE) azi1 = ls.Azimuth();
277  }
278  }
279 
280  // Max precision = 10: 0.1 nm in distance, 10^-15 deg (= 0.11 nm),
281  // 10^-11 sec (= 0.3 nm).
282  prec = std::min(10 + Math::extra_digits(), std::max(0, prec));
283  std::string s, eol, slat1, slon1, slat2, slon2, sazi1, ss12, strc;
284  std::istringstream str;
285  int retval = 0;
286  while (std::getline(*input, s)) {
287  try {
288  eol = "\n";
289  if (!cdelim.empty()) {
290  std::string::size_type m = s.find(cdelim);
291  if (m != std::string::npos) {
292  eol = " " + s.substr(m) + "\n";
293  s = s.substr(0, m);
294  }
295  }
296  str.clear(); str.str(s);
297  if (inverse) {
298  if (!(str >> slat1 >> slon1 >> slat2 >> slon2))
299  throw GeographicErr("Incomplete input: " + s);
300  if (str >> strc)
301  throw GeographicErr("Extraneous input: " + strc);
302  DMS::DecodeLatLon(slat1, slon1, lat1, lon1, longfirst);
303  DMS::DecodeLatLon(slat2, slon2, lat2, lon2, longfirst);
304  a12 = exact ?
305  geode.GenInverse(lat1, lon1, lat2, lon2, outmask,
306  s12, azi1, azi2, m12, M12, M21, S12) :
307  geods.GenInverse(lat1, lon1, lat2, lon2, outmask,
308  s12, azi1, azi2, m12, M12, M21, S12);
309  if (full) {
310  if (unroll) {
311  real e;
312  lon2 = lon1 + Math::AngDiff(lon1, lon2, e);
313  lon2 += e;
314  } else {
315  lon1 = Math::AngNormalize(lon1);
316  lon2 = Math::AngNormalize(lon2);
317  }
318  *output << LatLonString(lat1, lon1, prec, dms, dmssep, longfirst)
319  << " ";
320  }
321  *output << AzimuthString(azi1, prec, dms, dmssep) << " ";
322  if (full)
323  *output << LatLonString(lat2, lon2, prec, dms, dmssep, longfirst)
324  << " ";
325  if (azi2back)
326  azi2 += azi2 >= 0 ? -180 : 180;
327  *output << AzimuthString(azi2, prec, dms, dmssep) << " "
328  << DistanceStrings(s12, a12, full, arcmode, prec, dms);
329  if (full)
330  *output << " " << Utility::str(m12, prec)
331  << " " << Utility::str(M12, prec+7)
332  << " " << Utility::str(M21, prec+7)
333  << " " << Utility::str(S12, std::max(prec-7, 0));
334  *output << eol;
335  } else {
336  if (linecalc) {
337  if (!(str >> ss12))
338  throw GeographicErr("Incomplete input: " + s);
339  if (str >> strc)
340  throw GeographicErr("Extraneous input: " + strc);
341  // In fraction mode input is read as a distance
342  s12 = ReadDistance(ss12, !fraction && arcmode) * mult;
343  a12 = exact ?
344  le.GenPosition(arcmode, s12, outmask,
345  lat2, lon2, azi2, s12, m12, M12, M21, S12) :
346  ls.GenPosition(arcmode, s12, outmask,
347  lat2, lon2, azi2, s12, m12, M12, M21, S12);
348  } else {
349  if (!(str >> slat1 >> slon1 >> sazi1 >> ss12))
350  throw GeographicErr("Incomplete input: " + s);
351  if (str >> strc)
352  throw GeographicErr("Extraneous input: " + strc);
353  DMS::DecodeLatLon(slat1, slon1, lat1, lon1, longfirst);
354  azi1 = DMS::DecodeAzimuth(sazi1);
355  s12 = ReadDistance(ss12, arcmode);
356  a12 = exact ?
357  geode.GenDirect(lat1, lon1, azi1, arcmode, s12, outmask,
358  lat2, lon2, azi2, s12, m12, M12, M21, S12) :
359  geods.GenDirect(lat1, lon1, azi1, arcmode, s12, outmask,
360  lat2, lon2, azi2, s12, m12, M12, M21, S12);
361  }
362  if (full)
363  *output
364  << LatLonString(lat1, unroll ? lon1 : Math::AngNormalize(lon1),
365  prec, dms, dmssep, longfirst)
366  << " " << AzimuthString(azi1, prec, dms, dmssep) << " ";
367  if (azi2back)
368  azi2 += azi2 >= 0 ? -180 : 180;
369  *output << LatLonString(lat2, lon2, prec, dms, dmssep, longfirst)
370  << " " << AzimuthString(azi2, prec, dms, dmssep);
371  if (full)
372  *output << " "
373  << DistanceStrings(s12, a12, full, arcmode, prec, dms)
374  << " " << Utility::str(m12, prec)
375  << " " << Utility::str(M12, prec+7)
376  << " " << Utility::str(M21, prec+7)
377  << " " << Utility::str(S12, std::max(prec-7, 0));
378  *output << eol;
379  }
380  }
381  catch (const std::exception& e) {
382  // Write error message cout so output lines match input lines
383  *output << "ERROR: " << e.what() << "\n";
384  retval = 1;
385  }
386  }
387  return retval;
388  }
389  catch (const std::exception& e) {
390  std::cerr << "Caught exception: " << e.what() << "\n";
391  return 1;
392  }
393  catch (...) {
394  std::cerr << "Caught unknown exception\n";
395  return 1;
396  }
397 }
static T AngNormalize(T x)
Definition: Math.hpp:440
Matrix3f m
GeodesicLineExact InverseLine(real lat1, real lon1, real lat2, real lon2, unsigned caps=ALL) const
#define max(a, b)
Definition: datatypes.h:20
Header for GeographicLib::GeodesicLine class.
static Math::real DecodeAngle(const std::string &angstr)
Definition: src/DMS.cpp:281
float real
Definition: datatypes.h:10
int main(int argc, const char *const argv[])
Definition: GeodSolve.cpp:69
static const double lat
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
real GenInverse(real lat1, real lon1, real lat2, real lon2, unsigned outmask, real &s12, real &salp1, real &calp1, real &salp2, real &calp2, real &m12, real &M12, real &M21, real &S12) const
Header for GeographicLib::Utility class.
#define min(a, b)
Definition: datatypes.h:19
static T AngDiff(T x, T y, T &e)
Definition: Math.hpp:486
real ReadDistance(const std::string &s, bool arcmode)
Definition: GeodSolve.cpp:64
std::string LatLonString(real lat, real lon, int prec, bool dms, char dmssep, bool longfirst)
Definition: GeodSolve.cpp:33
Array33i a
static int extra_digits()
Definition: Math.hpp:187
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
Header for GeographicLib::Geodesic class.
static std::string Encode(real angle, component trailing, unsigned prec, flag ind=NONE, char dmssep=char(0))
Definition: src/DMS.cpp:299
static Math::real DecodeAzimuth(const std::string &azistr)
Definition: src/DMS.cpp:290
GeodesicLineExact Line(real lat1, real lon1, real azi1, unsigned caps=ALL) const
Header for GeographicLib::GeodesicLineExact class.
Math::real GenPosition(bool arcmode, real s12_a12, unsigned outmask, real &lat2, real &lon2, real &azi2, real &s12, real &m12, real &M12, real &M21, real &S12) const
#define GEOGRAPHICLIB_VERSION_STRING
Definition: Config.h:3
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Namespace for GeographicLib.
Math::real GenDistance(bool arcmode) const
Array< double, 1, 3 > e(1./3., 0.5, 2.)
RealScalar s
static std::string str(T x, int p=-1)
Definition: Utility.hpp:276
static void DecodeLatLon(const std::string &dmsa, const std::string &dmsb, real &lat, real &lon, bool longfirst=false)
Definition: src/DMS.cpp:252
static int set_digits(int ndigits=0)
Definition: Utility.cpp:48
Math::real Azimuth() const
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const ArgReturnType arg() const
Math::real GenDirect(real lat1, real lon1, real azi1, bool arcmode, real s12_a12, unsigned outmask, real &lat2, real &lon2, real &azi2, real &s12, real &m12, real &M12, real &M21, real &S12) const
Exact geodesic calculations.
Header for GeographicLib::GeodesicExact class.
GeodesicLineExact GenDirectLine(real lat1, real lon1, real azi1, bool arcmode, real s12_a12, unsigned caps=ALL) const
Exception handling for GeographicLib.
Definition: Constants.hpp:389
std::string AzimuthString(real azi, int prec, bool dms, char dmssep)
Definition: GeodSolve.cpp:45
Math::real GenPosition(bool arcmode, real s12_a12, unsigned outmask, real &lat2, real &lon2, real &azi2, real &s12, real &m12, real &M12, real &M21, real &S12) const
static const double lon
Math::real GenDistance(bool arcmode) const
std::string DistanceStrings(real s12, real a12, bool full, bool arcmode, int prec, bool dms)
Definition: GeodSolve.cpp:51
Geodesic calculations
Definition: Geodesic.hpp:172
EIGEN_DEVICE_FUNC const InverseReturnType inverse() const
Header for GeographicLib::DMS class.


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:42:08