[Proj] Area of a geodesic polygon

[Charles Karney]

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