249 lines
8.5 KiB
Plaintext
249 lines
8.5 KiB
Plaintext
=head1 NAME
|
|
|
|
Planimeter -- compute the area of geodesic polygons
|
|
|
|
=head1 SYNOPSIS
|
|
|
|
B<Planimeter> [ B<-r> ] [ B<-s> ] [ B<-l> ] [ B<-e> I<a> I<f> ]
|
|
[ B<-w> ] [ B<-p> I<prec> ] [ B<-G> | B<-E> | B<-Q> | B<-R> ]
|
|
[ B<--geoconvert-input> ]
|
|
[ B<--comment-delimiter> I<commentdelim> ]
|
|
[ B<--version> | B<-h> | B<--help> ]
|
|
[ B<--input-file> I<infile> | B<--input-string> I<instring> ]
|
|
[ B<--line-separator> I<linesep> ]
|
|
[ B<--output-file> I<outfile> ]
|
|
|
|
=head1 DESCRIPTION
|
|
|
|
Measure the area of a geodesic polygon. Reads polygon vertices from
|
|
standard input, one per line. Vertices are be given as latitude and
|
|
longitude. By default latitude precedes longitude, however this
|
|
convention is reversed with the B<-w> flag and a hemisphere designator
|
|
(I<N>, I<S>, I<E>, I<W>) can be used to disambiguate the coordinates.
|
|
The end of input, a blank line, or a line which can't be interpreted
|
|
as a vertex signals the end of one polygon and the start of the next.
|
|
For each polygon print a summary line with the number of points, the
|
|
perimeter (in meters), and the area (in meters^2).
|
|
|
|
The edges of the polygon are given by the I<shortest> geodesic (or
|
|
rhumb line) between consecutive vertices. In certain cases, there may
|
|
be two or many such shortest path, and in that case, the polygon is
|
|
not uniquely specified by its vertices. For geodesics, this only
|
|
happens with very long edges (for the WGS84 ellipsoid, any edge
|
|
shorter than 19970 km is uniquely specified by its end points). In
|
|
such cases, insert an additional vertex near the middle of the long
|
|
edge to define the boundary of the polygon.
|
|
|
|
By default, polygons traversed in a counter-clockwise direction return a
|
|
positive area and those traversed in a clockwise direction return a
|
|
negative area. This sign convention is reversed if the B<-r> option is
|
|
given.
|
|
|
|
Of course, encircling an area in the clockwise direction is equivalent
|
|
to encircling the rest of the ellipsoid in the counter-clockwise
|
|
direction. The default interpretation used by B<Planimeter> is the one
|
|
that results in a smaller magnitude of area; i.e., the magnitude of the
|
|
area is less than or equal to one half the total area of the ellipsoid.
|
|
If the B<-s> option is given, then the interpretation used is the one
|
|
that results in a positive area; i.e., the area is positive and less
|
|
than the total area of the ellipsoid.
|
|
|
|
Arbitrarily complex polygons are allowed. In the case of
|
|
self-intersecting polygons the area is accumulated "algebraically",
|
|
e.g., the areas of the 2 loops in a figure-8 polygon will partially
|
|
cancel. Polygons may include one or both poles. There is no need to
|
|
close the polygon.
|
|
|
|
=head1 OPTIONS
|
|
|
|
=over
|
|
|
|
=item B<-r>
|
|
|
|
toggle whether counter-clockwise traversal of the polygon returns a
|
|
positive (the default) or negative result.
|
|
|
|
=item B<-s>
|
|
|
|
toggle whether to return a signed result (the default) or not.
|
|
|
|
=item B<-l>
|
|
|
|
toggle whether the vertices represent a polygon (the default) or a
|
|
polyline. For a polyline, the number of points and the length of the
|
|
path joining them is returned; the path is not closed and the area is
|
|
not reported.
|
|
|
|
=item B<-e> I<a> I<f>
|
|
|
|
specify the ellipsoid via the equatorial radius, I<a> and
|
|
the flattening, I<f>. Setting I<f> = 0 results in a sphere. Specify
|
|
I<f> E<lt> 0 for a prolate ellipsoid. A simple fraction, e.g., 1/297,
|
|
is allowed for I<f>. By default, the WGS84 ellipsoid is used, I<a> =
|
|
6378137 m, I<f> = 1/298.257223563. If entering vertices as UTM/UPS or
|
|
MGRS coordinates, use the default ellipsoid, since the conversion of
|
|
these coordinates to latitude and longitude always uses the WGS84
|
|
parameters.
|
|
|
|
=item B<-w>
|
|
|
|
toggle the longitude first flag (it starts off); if the flag is on, then
|
|
when reading geographic coordinates, longitude precedes latitude (this
|
|
can be overridden by a hemisphere designator, I<N>, I<S>, I<E>, I<W>).
|
|
|
|
=item B<-p> I<prec>
|
|
|
|
set the output precision to I<prec> (default 6); the perimeter is given
|
|
(in meters) with I<prec> digits after the decimal point; the area is
|
|
given (in meters^2) with (I<prec> - 5) digits after the decimal point.
|
|
|
|
=item B<-G>
|
|
|
|
use the series formulation for the geodesics. This is the default
|
|
option and is recommended for terrestrial applications. This option,
|
|
B<-G>, and the following three options, B<-E>, B<-Q>, and B<-R>, are
|
|
mutually exclusive.
|
|
|
|
=item B<-E>
|
|
|
|
use "exact" algorithms (based on elliptic integrals) for the geodesic
|
|
calculations; the area is computed by applying discrete sine
|
|
transforms to the integrand in the expression for the area. These are
|
|
more accurate, albeit slower, than the (default) series expansions for
|
|
|I<f>| E<gt> 0.02 and will give high accuracy for -99 E<le> I<f> E<le>
|
|
0.99.
|
|
|
|
=item B<-Q>
|
|
|
|
perform the calculation on the authalic sphere. The area calculation is
|
|
accurate even if the flattening is large, I<provided> the edges are
|
|
sufficiently short. The perimeter calculation is not accurate.
|
|
|
|
=item B<-R>
|
|
|
|
The lines joining the vertices are rhumb lines instead of geodesics.
|
|
|
|
=item B<--geoconvert-input>
|
|
|
|
The input lines are interpreted in the same way as GeoConvert(1)
|
|
allowing the coordinates for the vertices to be given as UTM/UPS or
|
|
MGRS coordinates, as well as latitude and longitude. B<CAUTION>:
|
|
GeoConvert assumes the coordinates refer to the WGS84 ellipsoid
|
|
(disregarding the B<-e> flag) and MGRS coordinates signify the center
|
|
of the corresponding MGRS square.
|
|
|
|
=item B<--comment-delimiter> I<commentdelim>
|
|
|
|
set the comment delimiter to I<commentdelim> (e.g., "#" or "//"). If
|
|
set, the input lines will be scanned for this delimiter and, if found,
|
|
the delimiter and the rest of the line will be removed prior to
|
|
processing. For a given polygon, the last such string found will be
|
|
appended to the output line (separated by a space).
|
|
|
|
=item B<--version>
|
|
|
|
print version and exit.
|
|
|
|
=item B<-h>
|
|
|
|
print usage and exit.
|
|
|
|
=item B<--help>
|
|
|
|
print full documentation and exit.
|
|
|
|
=item B<--input-file> I<infile>
|
|
|
|
read input from the file I<infile> instead of from standard input; a file
|
|
name of "-" stands for standard input.
|
|
|
|
=item B<--input-string> I<instring>
|
|
|
|
read input from the string I<instring> instead of from standard input.
|
|
All occurrences of the line separator character (default is a semicolon)
|
|
in I<instring> are converted to newlines before the reading begins.
|
|
|
|
=item B<--line-separator> I<linesep>
|
|
|
|
set the line separator character to I<linesep>. By default this is a
|
|
semicolon.
|
|
|
|
=item B<--output-file> I<outfile>
|
|
|
|
write output to the file I<outfile> instead of to standard output; a
|
|
file name of "-" stands for standard output.
|
|
|
|
=back
|
|
|
|
=head1 EXAMPLES
|
|
|
|
Example (the area of the 100km MGRS square 18SWK)
|
|
|
|
Planimeter --geoconvert-input <<EOF
|
|
18n 500000 4400000
|
|
18n 600000 4400000
|
|
18n 600000 4500000
|
|
18n 500000 4500000
|
|
EOF
|
|
=> 4 400139.532959 10007388597.2
|
|
|
|
The following code takes the output from gdalinfo and reports the area
|
|
covered by the data (assuming the edges of the image are geodesics).
|
|
|
|
#! /bin/sh
|
|
egrep '^((Upper|Lower) (Left|Right)|Center) ' |
|
|
sed -e 's/d /d/g' -e "s/' /'/g" | tr -s '(),\r\t' ' ' | awk '{
|
|
if ($1 $2 == "UpperLeft")
|
|
ul = $6 " " $5;
|
|
else if ($1 $2 == "LowerLeft")
|
|
ll = $6 " " $5;
|
|
else if ($1 $2 == "UpperRight")
|
|
ur = $6 " " $5;
|
|
else if ($1 $2 == "LowerRight")
|
|
lr = $6 " " $5;
|
|
else if ($1 == "Center") {
|
|
printf "%s\n%s\n%s\n%s\n\n", ul, ll, lr, ur;
|
|
ul = ll = ur = lr = "";
|
|
}
|
|
}
|
|
' | Planimeter | cut -f3 -d' '
|
|
|
|
=head1 ACCURACY
|
|
|
|
Using the B<-G> option (the default), the accuracy was estimated by
|
|
computing the error in the area for 10^7 approximately regular
|
|
polygons on the WGS84 ellipsoid. The centers and the orientations of
|
|
the polygons were uniformly distributed, the number of vertices was
|
|
log-uniformly distributed in [3, 300], and the center to vertex
|
|
distance log-uniformly distributed in [0.1 m, 9000 km].
|
|
|
|
The maximum error in the perimeter was 200 nm, and the maximum error
|
|
in the area was
|
|
|
|
0.0013 m^2 for perimeter < 10 km
|
|
0.0070 m^2 for perimeter < 100 km
|
|
0.070 m^2 for perimeter < 1000 km
|
|
0.11 m^2 for all perimeters
|
|
|
|
=head1 SEE ALSO
|
|
|
|
GeoConvert(1), GeodSolve(1).
|
|
|
|
An online version of this utility is availbable at
|
|
L<https://geographiclib.sourceforge.io/cgi-bin/Planimeter>.
|
|
|
|
The algorithm for the area of geodesic polygon is
|
|
given in Section 6 of C. F. F. Karney, I<Algorithms for geodesics>,
|
|
J. Geodesy 87, 43-55 (2013);
|
|
DOI L<https://doi.org/10.1007/s00190-012-0578-z>;
|
|
addenda: L<https://geographiclib.sourceforge.io/geod-addenda.html>.
|
|
|
|
=head1 AUTHOR
|
|
|
|
B<Planimeter> was written by Charles Karney.
|
|
|
|
=head1 HISTORY
|
|
|
|
B<Planimeter> was added to GeographicLib,
|
|
L<https://geographiclib.sourceforge.io>, in version 1.4.
|