[Proj] Area of a geodesic polygon
Charles Karney
ckarney at sarnoff.com
Sun Sep 12 06:35:07 EST 2010
GeographicLib now includes a utility Planimeter for computing the areas
of geodesic polygons. See
http://geographiclib.sf.net/html/utilities.html#planimeter
The basic method entails computing the area "under" a geodesic, i.e.,
the area bounded by a geodesic, 2 meridians, and the equator and is
adapted from
J. Danielsen,
The Area under the Geodesic,
Survey Review 30 (232), 61-66 (1989)
who sets up a rigorous framework for the problem. I've extended his
series so that the truncation error is less than roundoff and rearranged
the terms to give a simple trigonometric series which is amenable to
Clenshaw summation. The roundoff error is typically 0.1 m^2 per vertex
or so.
A similar formulation is given by
L.E. Sjöberg,
Determination of areas on the plane, sphere and ellipsoid,
Survey Review 38 (301), 583-593 (July, 2006)
but this results in integrals which diverge for geodesics which pass
close to a pole.
In my reformulation of Danielsen's method, I express the area under the
geodesic, S12, as
S12 = S(sigma2) - S(sigma1)
S(sigma) = c^2*alpha + e^2*a^2 * cos(alpha0)*sin(alpha0) * I4(sigma)
I4(sigma) = - integrate(
(t(ep^2) - t(k^2*sin(sigma)^2)) / (ep^2 - k^2*sin(sigma)^2),
sigma)
t(x) = sqrt(1+1/x) * asinh(sqrt(x)) + x
a = major radius
c = authalic radius (area of ellipsoid = 4*pi*c^2)
e = eccentricity
ep = second eccentricity, ep^2 = e^2/(1-e^2)
alpha = azimuth of geodesic
alpha0 = azimuth at the equator crossing
k = ep*cos(alpha0)
sigma = arc length on auxiliary sphere (measured from the equator)
This expression for S12 is exact. I then expand I4 in a joint Taylor
series in ep^2 and k^2 to give a rapidly convergent expansion similar to
Danielsen's. The coefficients of the expansion C4[i] are given in
http://geographiclib.sf.net/html/geodesic.html#geodseries
The attractive feature of this method is that computing the area of an
polygon with N sides take a time which scales as O(N). The standard
alternative method is to transform the points to an equal-area
projection and compute the plane area. However this will take a time
which scales as the length of the perimeter (so that the geodesic sides
may be faithfully represented in the projected space). This comparison
supposes, of course, that the sides of the polygon are geodesics.
--
Charles Karney <ckarney at sarnoff.com>
Sarnoff Corporation, Princeton, NJ 08543-5300
Tel: +1 609 734 2312
Fax: +1 609 734 2662
More information about the Proj
mailing list