Area of Polygon on a Sphere - Algorithm
rgc
newsgroup in an algorithm for computing the area of a simple polygon on
a sphere. It turns out to be quite elegant. The description below
appears lengthy only because I have been annoyingly explicit and
complete.
PRELIMINARIES
Given a polygon with N edges, each of which is a segment of a great
circle. The polygon must be simply connected (a single piece), without
holes, and no edge may cross another. Neither pole may be inside the
polygon.
The polygon is described in a counterclockwise direction by a succession
of vertices, numbered from 0 to N-1. A vertex N, if given, is identical
to vertex 0. Note that "inside" and "outside" are defined by the
requirement that the polygon be counterclockwise. (If the polygon is
described in a clockwise direction, the area computed will be negative.)
The location of each vertex is given by latitude and longitude. In this
algorithm, latitude and longitude are expressed in radians.
Latitude is zero at the equator; north is positive, south is negative.
The latitude of point i is denoted by lat.i or lat.(i)
Longitude is zero at Greenwich: east is positive, west is negative. The
longitude of point i is denoted by lon.i or lon.(i)
With the radius of the Earth denoted by R, the area of the polygon,
expressed in the square of the units used for R, is given by:
THE ALGORITHM ITSELF
Area = - ( R^2 /2 ) * ( SUM )
where
SUM = ( lon.1 - lon.(N-1) ) * sin( lat.0 )
+ ( lon.2 - lon.0 ) * sin( lat.1 )
+ ( lon.3 - lon.1 ) * sin( lat.2 )
+ ( lon.4 - lon.2 ) * sin( lat.3 )
+ ...
+ ( lon.(N-1) - lon.(N-3) ) * sin( lat.(N-2) )
+ ( lon.0 - lon.(N-2) ) * sin( lat.(N-1) )
DISCLAIMERS
Other than that the Earth is assumed to be satisfactorily
approximated by a sphere throughout the area covered by the polygon, the
only simplifying assumption required is that the sine of the latitude of
each edge can be approximated throughout its length by the average of
the sines of the latitudes of its end points.
I worked out the formula for the exact solution on the sphere, but it
is horrendously complicated, involving inverse tangents of square roots
of products of tangents of sums of arc lengths that, in turn, require
the application of the haversine distance formula.
TEST RESULTS
We (actually, Will Duquette of JPL) conducted some tests, comparing
the areas of 1-degree and 0.01-degree "squares" on the sphere and on the
plane. The planar "square" was a rectangle with width equal to the side
arc length (in radians) times the cosine of the average latitude times
the radius of the Earth. Its height was the side arc length times the
radius of the Earth.
A one-degree "square" at the equator
Area by the above algorithm: 12,363.683990261115 sq km
Area of the planar "square": 12,363.840904434479 sq km
A one-degree "square" at latitude 35N
Area by the above algorithm: 10,065.850277260293 sq km
Area of the planar "square": 10,065.934954143779 sq km
A 0.01-degree "square" at the equator
Area by the above algorithm: 1.236431164871572 sq km
Area of the planar "square": 1.2364311664408989 sq km
A 0.01-degree "square" at latitude 35N
Area by the above algorithm: 1.0127632280177459 sq km
Area of the planar "square": 1.0127632288804711 sq km
These results show several things, only the first of which is a
surprise. (1) I apparently did not make an egregious mistake in algebra
during the derivation. (2) There is a larger relative difference (2
parts in a million vs 2 parts in a billion) between the areas of the two
kinds of "square" when the polygons are larger. (3) There is very little
difference even when the polygons are rather large (in my context) -
because their sides are still quite small compared to the radius of the
Earth.
Since we did not calculate the exact area, these results do not
actually indicate the error caused by the constant-sine approximation.
It could well be that the above algorithm gives results even closer to
the exact answer than it does to the planar approximation.
PLANAR FORMULA
Inspired by the simplicity of the above algorithm, I then worked out
the planar case to facilitate clarity in my write-up. The result is
indeed very similar.
Given: A planar polygon of N vertices, described counterclockwise by
successive points whose coordinates are (x.0,y.0), (x.1,y.1), (x.2,y.2),
..., (x.(N-1),y.(N-1)).
Area = - ( 1/2 ) * ( SUM )
where
SUM = ( x.1 - x.(N-1) ) * y.0
+ ( x.2 - x.0 ) * y.1
+ ( x.3 - x.1 ) * y.2
+ ( x.4 - x.2 ) * y.3
+ ...
+ ( x.(N-1) - x.(N-3) ) * y.(N-2)
+ ( x.0 - x.(N-2) ) * y.(N-1)
THE DERIVATIONS
If anybody is interested in the derivations, email me directly and I
will send you a copy.
Bob
Bob (Robert G.) Chamberlain |
r...@jpl.nasa.gov | If this were our worst problem,
Opinions & quips are mine | we could all sleep easily.
- or public - not JPL's | :-)