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