1=head1 NAME
2
3Planimeter -- compute the area of geodesic polygons
4
5=head1 SYNOPSIS
6
7B<Planimeter> [ B<-r> ] [ B<-s> ] [ B<-l> ] [ B<-e> I<a> I<f> ]
8[ B<-w> ] [ B<-p> I<prec> ] [ B<-G> | B<-E> | B<-Q> | B<-R> ]
9[ B<--comment-delimiter> I<commentdelim> ]
10[ B<--version> | B<-h> | B<--help> ]
11[ B<--input-file> I<infile> | B<--input-string> I<instring> ]
12[ B<--line-separator> I<linesep> ]
13[ B<--output-file> I<outfile> ]
14
15=head1 DESCRIPTION
16
17Measure the area of a geodesic polygon.  Reads polygon vertices from
18standard input, one per line.  Vertices may be given as latitude and
19longitude, UTM/UPS, or MGRS coordinates, interpreted in the same way as
20GeoConvert(1).  (MGRS coordinates signify the center of the
21corresponding MGRS square.)  The end of input, a blank line, or a line
22which can't be interpreted as a vertex signals the end of one polygon
23and the start of the next.  For each polygon print a summary line with
24the number of points, the perimeter (in meters), and the area (in
25meters^2).
26
27The edges of the polygon are given by the I<shortest> geodesic between
28consecutive vertices.  In certain cases, there may be two or many such
29shortest geodesics, and in that case, the polygon is not uniquely
30specified by its vertices.  This only happens with very long edges (for
31the WGS84 ellipsoid, any edge shorter than 19970 km is uniquely
32specified by its end points).  In such cases, insert an additional
33vertex near the middle of the long edge to define the boundary of the
34polygon.
35
36By default, polygons traversed in a counter-clockwise direction return a
37positive area and those traversed in a clockwise direction return a
38negative area.  This sign convention is reversed if the B<-r> option is
39given.
40
41Of course, encircling an area in the clockwise direction is equivalent
42to encircling the rest of the ellipsoid in the counter-clockwise
43direction.  The default interpretation used by B<Planimeter> is the one
44that results in a smaller magnitude of area; i.e., the magnitude of the
45area is less than or equal to one half the total area of the ellipsoid.
46If the B<-s> option is given, then the interpretation used is the one
47that results in a positive area; i.e., the area is positive and less
48than the total area of the ellipsoid.
49
50Arbitrarily complex polygons are allowed.  In the case of
51self-intersecting polygons the area is accumulated "algebraically",
52e.g., the areas of the 2 loops in a figure-8 polygon will partially
53cancel.  Polygons may include one or both poles.  There is no need to
54close the polygon.
55
56=head1 OPTIONS
57
58=over
59
60=item B<-r>
61
62toggle whether counter-clockwise traversal of the polygon returns a
63positive (the default) or negative result.
64
65=item B<-s>
66
67toggle whether to return a signed result (the default) or not.
68
69=item B<-l>
70
71toggle whether the vertices represent a polygon (the default) or a
72polyline.  For a polyline, the number of points and the length of the
73path joining them is returned; the path is not closed and the area is
74not reported.
75
76=item B<-e> I<a> I<f>
77
78specify the ellipsoid via the equatorial radius, I<a> and
79the flattening, I<f>.  Setting I<f> = 0 results in a sphere.  Specify
80I<f> E<lt> 0 for a prolate ellipsoid.  A simple fraction, e.g., 1/297,
81is allowed for I<f>.  By default, the WGS84 ellipsoid is used, I<a> =
826378137 m, I<f> = 1/298.257223563.  If entering vertices as UTM/UPS or
83MGRS coordinates, use the default ellipsoid, since the conversion of
84these coordinates to latitude and longitude always uses the WGS84
85parameters.
86
87=item B<-w>
88
89toggle the longitude first flag (it starts off); if the flag is on, then
90when reading geographic coordinates, longitude precedes latitude (this
91can be overridden by a hemisphere designator, I<N>, I<S>, I<E>, I<W>).
92
93=item B<-p> I<prec>
94
95set the output precision to I<prec> (default 6); the perimeter is given
96(in meters) with I<prec> digits after the decimal point; the area is
97given (in meters^2) with (I<prec> - 5) digits after the decimal point.
98
99=item B<-G>
100
101use the series formulation for the geodesics.  This is the default
102option and is recommended for terrestrial applications.  This option,
103B<-G>, and the following three options, B<-E>, B<-Q>, and B<-R>, are
104mutually exclusive.
105
106=item B<-E>
107
108use "exact" algorithms (based on elliptic integrals) for the geodesic
109calculations.  These are more accurate than the (default) series
110expansions for |I<f>| E<gt> 0.02.  (But note that the implementation of
111areas in GeodesicExact uses a high order series and this is only
112accurate for modest flattenings.)
113
114=item B<-Q>
115
116perform the calculation on the authalic sphere.  The area calculation is
117accurate even if the flattening is large, I<provided> the edges are
118sufficiently short.  The perimeter calculation is not accurate.
119
120=item B<-R>
121
122The lines joining the vertices are rhumb lines instead of geodesics.
123
124=item B<--comment-delimiter> I<commentdelim>
125
126set the comment delimiter to I<commentdelim> (e.g., "#" or "//").  If
127set, the input lines will be scanned for this delimiter and, if found,
128the delimiter and the rest of the line will be removed prior to
129processing.  For a given polygon, the last such string found will be
130appended to the output line (separated by a space).
131
132=item B<--version>
133
134print version and exit.
135
136=item B<-h>
137
138print usage and exit.
139
140=item B<--help>
141
142print full documentation and exit.
143
144=item B<--input-file> I<infile>
145
146read input from the file I<infile> instead of from standard input; a file
147name of "-" stands for standard input.
148
149=item B<--input-string> I<instring>
150
151read input from the string I<instring> instead of from standard input.
152All occurrences of the line separator character (default is a semicolon)
153in I<instring> are converted to newlines before the reading begins.
154
155=item B<--line-separator> I<linesep>
156
157set the line separator character to I<linesep>.  By default this is a
158semicolon.
159
160=item B<--output-file> I<outfile>
161
162write output to the file I<outfile> instead of to standard output; a
163file name of "-" stands for standard output.
164
165=back
166
167=head1 EXAMPLES
168
169Example (the area of the 100km MGRS square 18SWK)
170
171   Planimeter <<EOF
172   18n 500000 4400000
173   18n 600000 4400000
174   18n 600000 4500000
175   18n 500000 4500000
176   EOF
177   => 4 400139.53295860 10007388597.1913
178
179The following code takes the output from gdalinfo and reports the area
180covered by the data (assuming the edges of the image are geodesics).
181
182   #! /bin/sh
183   egrep '^((Upper|Lower) (Left|Right)|Center) ' |
184   sed -e 's/d /d/g' -e "s/' /'/g" | tr -s '(),\r\t' ' ' | awk '{
185       if ($1 $2 == "UpperLeft")
186           ul = $6 " " $5;
187       else if ($1 $2 == "LowerLeft")
188           ll = $6 " " $5;
189       else if ($1 $2 == "UpperRight")
190           ur = $6 " " $5;
191       else if ($1 $2 == "LowerRight")
192           lr = $6 " " $5;
193       else if ($1 == "Center") {
194           printf "%s\n%s\n%s\n%s\n\n", ul, ll, lr, ur;
195           ul = ll = ur = lr = "";
196       }
197   }
198   ' | Planimeter | cut -f3 -d' '
199
200=head1 ACCURACY
201
202Using the B<-G> option (the default), the accuracy was estimated by
203computing the error in the area for 10^7 approximately regular
204polygons on the WGS84 ellipsoid.  The centers and the orientations of
205the polygons were uniformly distributed, the number of vertices was
206log-uniformly distributed in [3, 300], and the center to vertex
207distance log-uniformly distributed in [0.1 m, 9000 km].
208
209The maximum error in the perimeter was 200 nm, and the maximum error
210in the area was
211
212   0.0013 m^2 for perimeter < 10 km
213   0.0070 m^2 for perimeter < 100 km
214   0.070 m^2 for perimeter < 1000 km
215   0.11 m^2 for all perimeters
216
217=head1 SEE ALSO
218
219GeoConvert(1), GeodSolve(1).
220
221An online version of this utility is availbable at
222L<https://geographiclib.sourceforge.io/cgi-bin/Planimeter>.
223
224The algorithm for the area of geodesic polygon is
225given in Section 6 of C. F. F. Karney, I<Algorithms for geodesics>,
226J. Geodesy 87, 43-55 (2013);
227DOI L<https://doi.org/10.1007/s00190-012-0578-z>;
228addenda: L<https://geographiclib.sourceforge.io/geod-addenda.html>.
229
230=head1 AUTHOR
231
232B<Planimeter> was written by Charles Karney.
233
234=head1 HISTORY
235
236B<Planimeter> was added to GeographicLib,
237L<https://geographiclib.sourceforge.io>, in version 1.4.
238