WRAPPER_API double PolygonArea(const CoordPointCollection &points)
The PolygonArea method returns the area of a polygon on the surface of the Ellipsoid. The only argument is a CoordPointCollection, which should contain only GeodeticPoints. The method returns the polygon's area in square meters.
These points indicate the locations of the polygon's vertices, in order. The first point should not be duplicated at the end of the point collection (closure of the set is assumed). The polygon must be contiguous and cannot have holes. Areas of polygons with multiple rings or holes must be found through multiple calls to this function. Similarly, polygons whose boundaries are self-intersecting may produce unexpected results. Clockwise oriented sets of points will be reported as having positive area, while counter-clockwise sets will have negative area. The function does not modify the input points.
This function calculates the polygon's area directly from GeodeticPoints. The polygon's sides are assumed to be the shortest-path geodesic between consecutive vertices. It does not account for differences between the geoid and the Ellipsoid. Similarly, elevation values are not considered (all points are treated as lying on the surface of the Ellipsoid). For polygons where the geodesic interpretation of sides is reasonable, geodetic area calculation will typically produce more stable results than simply projecting the vertices and calculating a planar area. This is because a projection will distort either area or direction (or both) to some extent. Obviously, area distortion limits the effectiveness of projective area calculation. If there is significant separation between the vertices, directional distortion will may also have a noticeable impact on projective area calculation. This is because the projection will, in effect, draw the wrong line between vertices, altering the shape and thus the area of the polygon. For extremely dense polygons, the directional distortion becomes less of an issue.
The type of geodetic area calculation used in this function is susceptible to increased error for some point configurations. There are two characteristics in particular that can increase the expected area. The first characteristic is difference between the latitude range of the input points and latitude range of the polygon itself. This characteristic is most pronounced in "polar polygons" that encircle the north or south pole (the boundary of Antarctica for example), but is also present to some extent with widely spaced consecutive vertices near its upper or lower edge. The second characteristic is having area that is unusually weighted toward or away from the equator within the range of latitude spanned by the polygon. Put another way, does the polygon more fully fill some strips of latitude than others? For example, consider a polygon at the equator, with a very long skinny extension reaching up to a pole. Note that these characteristics may co-exist to some extent. As an illustration, image a "polar dumbell" consisting of two near-circles around each pole connected by a very thin strip along a line of latitude. This polygon would span all latitude values, even though its boundary points may never be very near to a pole; and its area would be strongly concentrated toward the poles and away from the equator, even though it would span all latitudes. The worst case theoretical error introduced by these considerations varies according to a complex formula based on the Ellipsoid's semi-major axis and flattening, but would be less than half a percent for the WGS84 ellipsoid. For most polygons, however, this theoretical error should be quite small.
void Ellipsoid_PolygonArea(GEOCALCPBW_NAMESPACE::DataSource & data)
{
GEOCALCPBW_NAMESPACE::Ellipsoid ell;
GEOCALCPBW_NAMESPACE::LinearValue semiMaj;
semiMaj.set_InMeters(6378137);
ell.set_SemiMajor(semiMaj);
ell.set_InvFlatDefinitive(true);
ell.set_InverseFlattening(298.257223563);
GEOCALCPBW_NAMESPACE::GeodeticPoint * pt = data.GetGeodeticPoint(L"BMG", L"GEODETIC_POINT_DEGREES");
GEOCALCPBW_NAMESPACE::CoordPointCollection vertices;
vertices.set_PointStyle(*pt);
pt->set_InUnits(0, 0);
vertices.AddPoint(*pt);
pt->set_InUnits(-5, 5);
vertices.AddPoint(*pt);
pt->set_InUnits(0, 10);
vertices.AddPoint(*pt);
pt->set_InUnits(6, 6);
vertices.AddPoint(*pt);
pt->set_InUnits(5, 2);
vertices.AddPoint(*pt);
double area = ell.PolygonArea(vertices);
delete pt;
}