area_poly1.c

Go to the documentation of this file.

/*!
 * \file lib/gis/area_poly1.c
 *
 * \brief GIS Library - Polygon area calculation routines.
 *
 * SPDX-FileCopyrightText: 2001-2013 GRASS Development Team
 * SPDX-License-Identifier: GPL-2.0-or-later
 *
 * \author Original author CERL
 */
 
#include <math.h>
#include <grass/gis.h>
#include <geodesic.h>
#include "pi.h"
 
#define TWOPI             M_PI + M_PI
 
#define USE_GEOGRAPHICLIB 1
 
static struct state {
#ifndef USE_GEOGRAPHICLIB
    double QA, QB, QC;
    double QbarA, QbarB, QbarC, QbarD;
 
    double AE; /** a^2(1-e^2) */
 
    double Qp; /** Q at the north pole */
 
    double E; /** Area of the earth */
 
#else
    struct geod_geodesic g; /** GeographicLib */
#endif
} state;
 
static struct state *st = &state;
 
#ifndef USE_GEOGRAPHICLIB
static double Q(double x)
{
    double sinx, sinx2;
 
    sinx = sin(x);
    sinx2 = sinx * sinx;
 
    return sinx * (1 + sinx2 * (st->QA + sinx2 * (st->QB + sinx2 * st->QC)));
}
 
static double Qbar(double x)
{
    double cosx, cosx2;
 
    cosx = cos(x);
    cosx2 = cosx * cosx;
 
    return cosx *
           (st->QbarA +
            cosx2 * (st->QbarB + cosx2 * (st->QbarC + cosx2 * st->QbarD)));
}
#endif
 
/*!
 * \brief Begin area calculations.
 *
 * This initializes the polygon area calculations for the
 * ellipsoid with semi-major axis <i>a</i> (in meters) and ellipsoid
 * eccentricity squared <i>e2</i>.
 *
 * \param a semi-major axis
 * \param e2 ellipsoid eccentricity squared
 */
void G_begin_ellipsoid_polygon_area(double a, double e2)
{
#ifdef USE_GEOGRAPHICLIB
    double f;
 
    /* GeographicLib */
    G_get_ellipsoid_parameters(&a, &e2);
    f = 1.0 - sqrt(1.0 - e2);
    geod_init(&st->g, a, f);
 
#else
    double e4, e6;
 
    /* GRASS native */
    e4 = e2 * e2;
    e6 = e4 * e2;
 
    st->AE = a * a * (1 - e2);
 
    st->QA = (2.0 / 3.0) * e2;
    st->QB = (3.0 / 5.0) * e4;
    st->QC = (4.0 / 7.0) * e6;
 
    st->QbarA = -1.0 - (2.0 / 3.0) * e2 - (3.0 / 5.0) * e4 - (4.0 / 7.0) * e6;
    st->QbarB = (2.0 / 9.0) * e2 + (2.0 / 5.0) * e4 + (4.0 / 7.0) * e6;
    st->QbarC = -(3.0 / 25.0) * e4 - (12.0 / 35.0) * e6;
    st->QbarD = (4.0 / 49.0) * e6;
 
    st->Qp = Q(M_PI_2);
    st->E = 4 * M_PI * st->Qp * st->AE;
    if (st->E < 0.0)
        st->E = -st->E;
#endif
}
 
/*!
 * \brief Area of lat-long polygon.
 *
 * Returns the area in square meters of the polygon described by the
 * <i>n</i> pairs of <i>lat,long</i> vertices for latitude-longitude
 * grids.
 *
 * <b>Note:</b> This routine computes the area of a polygon on the
 * ellipsoid.  The sides of the polygon are rhumb lines and, in general,
 * not geodesics.  Each side is actually defined by a linear relationship
 * between latitude and longitude, i.e., on a rectangular/equidistant
 * cylindrical/Plate Carr{'e}e grid, the side would appear as a
 * straight line.  For two consecutive vertices of the polygon,
 * (lat_1, long1) and (lat_2,long_2), the line joining them (i.e., the
 * polygon's side) is defined by:
 *
 \verbatim
 lat_2  -  lat_1
 lat = lat_1 + (long - long_1) * ---------------
 long_2 - long_1
 \endverbatim
 *
 * where long_1 < long < long_2.
 *   The values of QbarA, etc., are determined by the integration of
 * the Q function.  Into www.integral-calculator.com, paste this
 * expression :
 *
 \verbatim
 sin(x)+ (2/3)e^2(sin(x))^3 + (3/5)e^4(sin(x))^5 + (4/7)e^6(sin(x))^7
 \endverbatim
 *
 * and you'll get their values.  (Last checked 30 Oct 2013).
 *
 * This function correctly computes (within the limits of the series
 * approximation) the area of a quadrilateral on the ellipsoid when
 * two of its sides run along meridians and the other two sides run
 * along parallels of latitude.
 *
 * \param lon array of longitudes
 * \param lat array of latitudes
 * \param n number of lat,lon pairs
 *
 * \return area in square meters
 */
double G_ellipsoid_polygon_area(const double *lon, const double *lat, int n)
{
    double area;
 
#ifdef USE_GEOGRAPHICLIB
    /* GeographicLib */
    struct geod_polygon p;
    double pP;
    int i;
 
    geod_polygon_init(&p, FALSE);
    /* GeographicLib does not need a closed ring,
     * see example for geod_polygonarea() in geodesic.h */
    /* add points in reverse order */
    i = n;
    while (--i)
        geod_polygon_addpoint(&st->g, &p, (double)lat[i], (double)lon[i]);
    /* The area returned is signed
     * with counter-clockwise traversal being treated as positive. */
    geod_polygon_compute(&st->g, &p, FALSE, TRUE, &area, &pP);
    area = fabs(area);
 
#else
    double x1, y1, x2, y2, dx, dy;
    double Qbar1, Qbar2;
    /* threshold for dy, should be between 1e-4 and 1e-7 */
    double thresh = 1e-6;
 
    /* GRASS native */
    x2 = Radians(lon[n - 1]);
    y2 = Radians(lat[n - 1]);
    Qbar2 = Qbar(y2);
 
    area = 0.0;
 
    while (--n >= 0) {
        x1 = x2;
        y1 = y2;
        Qbar1 = Qbar2;
 
        x2 = Radians(*lon++);
        y2 = Radians(*lat++);
        Qbar2 = Qbar(y2);
 
        if (x1 > x2)
            while (x1 - x2 > M_PI)
                x2 += TWOPI;
        else if (x2 > x1)
            while (x2 - x1 > M_PI)
                x1 += TWOPI;
 
        dx = x2 - x1;
        dy = y2 - y1;
 
        if (fabs(dy) > thresh) {
            /* account for different latitudes y1, y2 */
            area += dx * (st->Qp - (Qbar2 - Qbar1) / dy);
            /* original:
             * area += dx * st->Qp - (dx / dy) * (Qbar2 - Qbar1);
             */
        }
        else {
            /* latitudes y1, y2 are (nearly) identical */
            /* if y2 becomes similar to y1, i.e. y2 -> y1
             * Qbar2 - Qbar1 -> 0 and dy -> 0
             * (Qbar2 - Qbar1) / dy -> ?
             * (Qbar2 - Qbar1) / dy should approach Q((y1 + y2) / 2)
             * Metz 2017
             */
            area += dx * (st->Qp - Q((y1 + y2) / 2));
        }
    }
    if ((area *= st->AE) < 0.0)
        area = -area;
 
    /* kludge - if polygon circles the south pole the area will be
     * computed as if it circled the north pole. The correction is
     * the difference between total surface area of the earth and
     * the "north pole" area.
     */
    if (area > st->E)
        area = st->E;
    if (area > st->E / 2)
        area = st->E - area;
#endif
 
    return area;
}