[OpenLayers-Trac] [OpenLayers] #3427: Area of a spherical polygon

[OpenLayers]

#3427: Area of a spherical polygon
-------------------------+--------------------------------------------------
 Reporter:  karney       |       Owner:              
     Type:  feature      |      Status:  new         
 Priority:  minor        |   Milestone:  2.12 Release
Component:  general      |     Version:  2.10        
 Keywords:  sphere area  |       State:              
-------------------------+--------------------------------------------------
 Ticket #3426 includes a patch to give you the area of a geodesic polygon
 on the ellipsoid and points out the limitations of the existing
 getGeodesicArea function.

 If you want to retain a simple function to compute the spherical area
 (getSphericalArea maybe), then use the following code.  This computes
 the area of a spherical polygon whose edges are great circles.  This is
 an improvement over the current getGeodesicArea where the edges are some
 weirder lines.  It keeps the limitation of getGeodesicArea in that the
 polygon can't encircle a pole.  (This could be fixed...  but then you'd
 be better off using the general ellipsoidal routine.)
 {{{
     /**
      * APIMethod: getSphericalArea
      * Calculate the area of the polygon were it projected onto a
      *     spherical earth.  Note that this area will be positive if
      *     ring is oriented clockwise, otherwise it will be negative.
      *     The polygon should not encircle a pole.
      *
      * Parameters:
      * projection - {<OpenLayers.Projection>} The spatial reference system
      *     for the geometry coordinates.  If not provided,
 Geographic/WGS84 is
      *     assumed.
      *
      * Reference:
      * F. W. Bessel, "The calculation of longitude and latitude from
      *     geodesic mesurements", Astron. Nachr. 331(8) 852-861 (2010),
      *     Sec. 11.  Preprint http://arxiv.org/abs/0908.1824
      *
      * An earlier version of this routine (misleadingly called
      * getGeodesicArea) used a formula from Chamberlain and Duquette.
      * However this returned the area of the polygon in a cylindrical
      * equal area projection (and great circles don't project to
      * straight lines in that projection).
      *
      * Returns:
      * {float} The signed spherical area of the polygon in square meters.
      */
     getSphericalArea: function(projection) {
         var ring = this;  // so we can work with a clone if needed
         if(projection) {
             var gg = new OpenLayers.Projection("EPSG:4326");
             if(!gg.equals(projection)) {
                 ring = this.clone().transform(projection, gg);
             }
         }
         var area = 0.0;
         var len = ring.components && ring.components.length;
         if(len > 2) {
             var p1, p2;
             // The term added to area here is 0.5 * spherical excess for
             // the quadrilateral bounded by the great circle, 2
             // meridians, and the equator.  The excess is the difference
             // of the two forward azimuths alpha2 - alpha1.  Use
             // Bessel's equation for cot((alpha' - alpha)/2)
             // substituting alpha' = alpha1, and alpha = alpha2 + pi
             // (see Figure in the paper).  Summing this over all the
             // edges of a polygon given the 0.5 * spherical area
             // provided the polygon doesn't encircle a pole.
             for(var i=0; i<len-1; i++) {
                 p1 = ring.components[i];
                 p2 = ring.components[i+1];
                 area += Math.atan(
                     Math.tan(OpenLayers.Util.rad(p2.x - p1.x) / 2) *
                     Math.sin(OpenLayers.Util.rad(p2.y + p1.y) / 2) /
                     Math.cos(OpenLayers.Util.rad(p2.y - p1.y) / 2) );
             }
             // Factor of 2 accounts for the 0.5 in the explanation above.
             area = area * 2 * 6378137.0 * 6378137.0;
         }
         return area;
     },
 }}}

-- 
Ticket URL: <http://trac.openlayers.org/ticket/3427>
OpenLayers <http://openlayers.org/>
A free AJAX map viewer