GeographicLib 2.5
Planimeter.cpp
Go to the documentation of this file.
1/**
2 * \file Planimeter.cpp
3 * \brief Command line utility for measuring the area of geodesic polygons
4 *
5 * Copyright (c) Charles Karney (2010-2023) <karney@alum.mit.edu> and licensed
6 * under the MIT/X11 License. For more information, see
7 * https://geographiclib.sourceforge.io/
8 *
9 * See the <a href="Planimeter.1.html">man page</a> for usage information.
10 **********************************************************************/
11
12#include <iostream>
13#include <string>
14#include <sstream>
15#include <fstream>
17#include <GeographicLib/DMS.hpp>
21
22#include "Planimeter.usage"
23
24int main(int argc, const char* const argv[]) {
25 try {
26 using namespace GeographicLib;
27 typedef Math::real real;
28 Utility::set_digits();
29 enum { GEODESIC, AUTHALIC, RHUMB };
30 real
31 a = Constants::WGS84_a(),
32 f = Constants::WGS84_f();
33 bool reverse = false, sign = true, polyline = false, longfirst = false,
34 exact = false, geoconvert_compat = false;
35 int linetype = GEODESIC;
36 int prec = 6;
37 std::string istring, ifile, ofile, cdelim;
38 char lsep = ';';
39
40 for (int m = 1; m < argc; ++m) {
41 std::string arg(argv[m]);
42 if (arg == "-r")
43 reverse = !reverse;
44 else if (arg == "-s")
45 sign = !sign;
46 else if (arg == "-l")
47 polyline = !polyline;
48 else if (arg == "-e") {
49 if (m + 2 >= argc) return usage(1, true);
50 try {
51 a = Utility::val<real>(std::string(argv[m + 1]));
52 f = Utility::fract<real>(std::string(argv[m + 2]));
53 }
54 catch (const std::exception& e) {
55 std::cerr << "Error decoding arguments of -e: " << e.what() << "\n";
56 return 1;
57 }
58 m += 2;
59 } else if (arg == "-w")
60 longfirst = !longfirst;
61 else if (arg == "-p") {
62 if (++m == argc) return usage(1, true);
63 try {
64 prec = Utility::val<int>(std::string(argv[m]));
65 }
66 catch (const std::exception&) {
67 std::cerr << "Precision " << argv[m] << " is not a number\n";
68 return 1;
69 }
70 } else if (arg == "-G")
71 linetype = GEODESIC;
72 else if (arg == "-Q")
73 linetype = AUTHALIC;
74 else if (arg == "-R")
75 linetype = RHUMB;
76 else if (arg == "-E")
77 exact = true;
78 else if (arg == "--geoconvert-input")
79 geoconvert_compat = true;
80 else if (arg == "--input-string") {
81 if (++m == argc) return usage(1, true);
82 istring = argv[m];
83 } else if (arg == "--input-file") {
84 if (++m == argc) return usage(1, true);
85 ifile = argv[m];
86 } else if (arg == "--output-file") {
87 if (++m == argc) return usage(1, true);
88 ofile = argv[m];
89 } else if (arg == "--line-separator") {
90 if (++m == argc) return usage(1, true);
91 if (std::string(argv[m]).size() != 1) {
92 std::cerr << "Line separator must be a single character\n";
93 return 1;
94 }
95 lsep = argv[m][0];
96 } else if (arg == "--comment-delimiter") {
97 if (++m == argc) return usage(1, true);
98 cdelim = argv[m];
99 } else if (arg == "--version") {
100 std::cout << argv[0] << ": GeographicLib version "
101 << GEOGRAPHICLIB_VERSION_STRING << "\n";
102 return 0;
103 } else
104 return usage(!(arg == "-h" || arg == "--help"), arg != "--help");
105 }
106
107 if (!ifile.empty() && !istring.empty()) {
108 std::cerr << "Cannot specify --input-string and --input-file together\n";
109 return 1;
110 }
111 if (ifile == "-") ifile.clear();
112 std::ifstream infile;
113 std::istringstream instring;
114 if (!ifile.empty()) {
115 infile.open(ifile.c_str());
116 if (!infile.is_open()) {
117 std::cerr << "Cannot open " << ifile << " for reading\n";
118 return 1;
119 }
120 } else if (!istring.empty()) {
121 std::string::size_type m = 0;
122 while (true) {
123 m = istring.find(lsep, m);
124 if (m == std::string::npos)
125 break;
126 istring[m] = '\n';
127 }
128 instring.str(istring);
129 }
130 std::istream* input = !ifile.empty() ? &infile :
131 (!istring.empty() ? &instring : &std::cin);
132
133 std::ofstream outfile;
134 if (ofile == "-") ofile.clear();
135 if (!ofile.empty()) {
136 outfile.open(ofile.c_str());
137 if (!outfile.is_open()) {
138 std::cerr << "Cannot open " << ofile << " for writing\n";
139 return 1;
140 }
141 }
142 std::ostream* output = !ofile.empty() ? &outfile : &std::cout;
143
144 // Linetype is one of GEODESIC, AUTHALIC, RHUMB
145 const AuxLatitude ellip(a, f);
146 if (linetype == AUTHALIC) {
147 // adjusting a and f to correspond to the authalic sphere
148 using std::sqrt;
149 a = sqrt(ellip.AuthalicRadiusSquared(exact));
150 f = 0;
151 }
152 const Geodesic geod(a, linetype != RHUMB ? f : 0,
153 exact && linetype != RHUMB);
154 const Rhumb rhumb(a, linetype == RHUMB ? f : 0,
155 exact && linetype == RHUMB);
156 PolygonArea poly(geod, polyline);
157 PolygonAreaRhumb polyr(rhumb, polyline);
158 GeoCoords p;
159
160 // Max precision = 10: 0.1 nm in distance, 10^-15 deg (= 0.11 nm),
161 // 10^-11 sec (= 0.3 nm).
162 prec = std::min(10 + Math::extra_digits(), std::max(0, prec));
163 std::string s, eol("\n");
164 real perimeter, area;
165 unsigned num;
166 std::istringstream str;
167 std::string slat, slon, junk;
168 real lat = 0, lon = 0;
169 while (std::getline(*input, s)) {
170 if (!cdelim.empty()) {
171 std::string::size_type m = s.find(cdelim);
172 if (m != std::string::npos) {
173 eol = " " + s.substr(m) + "\n";
174 s = s.substr(0, m);
175 }
176 }
177 bool endpoly = s.empty();
178 if (!endpoly) {
179 try {
180 using std::isnan;
181 if (geoconvert_compat) {
182 p.Reset(s, true, longfirst);
183 lat = p.Latitude(); lon = p.Longitude();
184 } else {
185 str.clear(); str.str(s);
186 if (!(str >> slat >> slon))
187 throw GeographicErr("incomplete input");
188 if (str >> junk)
189 throw GeographicErr("extra input");
190 DMS::DecodeLatLon(slat, slon, lat, lon, longfirst);
191 }
192 if (isnan(lat) || isnan(lon))
193 endpoly = true;
194 }
195 catch (const GeographicErr&) {
196 endpoly = true;
197 }
198 }
199 if (endpoly) {
200 num =
201 linetype == RHUMB ? polyr.Compute(reverse, sign, perimeter, area) :
202 poly.Compute(reverse, sign, perimeter, area); // geodesic + authalic
203 if (num > 0) {
204 *output << num << " " << Utility::str(perimeter, prec);
205 if (!polyline) {
206 *output << " " << Utility::str(area, std::max(0, prec - 5));
207 }
208 *output << eol;
209 }
210 linetype == RHUMB ? polyr.Clear() : poly.Clear();
211 eol = "\n";
212 } else {
213 linetype == RHUMB ? polyr.AddPoint(lat, lon) :
214 poly.AddPoint
215 (linetype == AUTHALIC ?
216 ellip.Convert(AuxLatitude::PHI, AuxLatitude::XI, lat, exact) : lat,
217 lon);
218 }
219 }
220 num =
221 linetype == RHUMB ? polyr.Compute(reverse, sign, perimeter, area) :
222 poly.Compute(reverse, sign, perimeter, area);
223 if (num > 0) {
224 *output << num << " " << Utility::str(perimeter, prec);
225 if (!polyline) {
226 *output << " " << Utility::str(area, std::max(0, prec - 5));
227 }
228 *output << eol;
229 }
230 linetype == RHUMB ? polyr.Clear() : poly.Clear();
231 eol = "\n";
232 return 0;
233 }
234 catch (const std::exception& e) {
235 std::cerr << "Caught exception: " << e.what() << "\n";
236 return 1;
237 }
238 catch (...) {
239 std::cerr << "Caught unknown exception\n";
240 return 1;
241 }
242}
Header for the GeographicLib::AuxLatitude class.
Header for GeographicLib::DMS class.
Header for GeographicLib::GeoCoords class.
GeographicLib::Math::real real
Definition GeodSolve.cpp:28
int main(int argc, const char *const argv[])
Header for GeographicLib::PolygonAreaT class.
Header for GeographicLib::Utility class.
Conversions between auxiliary latitudes.
AuxAngle Convert(int auxin, int auxout, const AuxAngle &zeta, bool exact=false) const
Math::real AuthalicRadiusSquared(bool exact=false) const
Conversion between geographic coordinates.
Definition GeoCoords.hpp:49
void Reset(const std::string &s, bool centerp=true, bool longfirst=false)
Definition GeoCoords.cpp:19
Math::real Latitude() const
Math::real Longitude() const
Geodesic calculations
Definition Geodesic.hpp:175
Exception handling for GeographicLib.
Solve of the direct and inverse rhumb problems.
Definition Rhumb.hpp:80
Namespace for GeographicLib.