GRASS GIS 8 Programmer's Manual 8.2.1RC1(2022)-exported
area_poly1.c
Go to the documentation of this file.
1/*!
2 * \file lib/gis/area_poly1.c
3 *
4 * \brief GIS Library - Polygon area calculation routines.
5 *
6 * (C) 2001-2013 by the GRASS Development Team
7 *
8 * This program is free software under the GNU General Public License
9 * (>=v2). Read the file COPYING that comes with GRASS for details.
10 *
11 * \author Original author CERL
12 */
13
14#include <math.h>
15#include <grass/gis.h>
16#include "pi.h"
17
18#define TWOPI M_PI + M_PI
19
20static struct state {
21 double QA, QB, QC;
22 double QbarA, QbarB, QbarC, QbarD;
23 double AE; /** a^2(1-e^2) */
24 double Qp; /** Q at the north pole */
25 double E; /** Area of the earth */
26} state;
27
28static struct state *st = &state;
29
30static double Q(double x)
31{
32 double sinx, sinx2;
33
34 sinx = sin(x);
35 sinx2 = sinx * sinx;
36
37 return sinx * (1 + sinx2 * (st->QA + sinx2 * (st->QB + sinx2 * st->QC)));
38}
39
40static double Qbar(double x)
41{
42 double cosx, cosx2;
43
44 cosx = cos(x);
45 cosx2 = cosx * cosx;
46
47 return cosx * (st->QbarA + cosx2 * (st->QbarB + cosx2 * (st->QbarC + cosx2 * st->QbarD)));
48}
49
50/*!
51 * \brief Begin area calculations.
52 *
53 * This initializes the polygon area calculations for the
54 * ellipsoid with semi-major axis <i>a</i> (in meters) and ellipsoid
55 * eccentricity squared <i>e2</i>.
56 *
57 * \param a semi-major axis
58 * \param e2 ellipsoid eccentricity squared
59 */
60
61void G_begin_ellipsoid_polygon_area(double a, double e2)
62{
63 double e4, e6;
64
65 e4 = e2 * e2;
66 e6 = e4 * e2;
67
68 st->AE = a * a * (1 - e2);
69
70 st->QA = (2.0 / 3.0) * e2;
71 st->QB = (3.0 / 5.0) * e4;
72 st->QC = (4.0 / 7.0) * e6;
73
74 st->QbarA = -1.0 - (2.0 / 3.0) * e2 - (3.0 / 5.0) * e4 - (4.0 / 7.0) * e6;
75 st->QbarB = (2.0 / 9.0) * e2 + (2.0 / 5.0) * e4 + (4.0 / 7.0) * e6;
76 st->QbarC = -(3.0 / 25.0) * e4 - (12.0 / 35.0) * e6;
77 st->QbarD = (4.0 / 49.0) * e6;
78
79 st->Qp = Q(M_PI_2);
80 st->E = 4 * M_PI * st->Qp * st->AE;
81 if (st->E < 0.0)
82 st->E = -st->E;
83}
84
85/*!
86 * \brief Area of lat-long polygon.
87 *
88 * Returns the area in square meters of the polygon described by the
89 * <i>n</i> pairs of <i>lat,long</i> vertices for latitude-longitude
90 * grids.
91 *
92 * <b>Note:</b> This routine computes the area of a polygon on the
93 * ellipsoid. The sides of the polygon are rhumb lines and, in general,
94 * not geodesics. Each side is actually defined by a linear relationship
95 * between latitude and longitude, i.e., on a rectangular/equidistant
96 * cylindrical/Plate Carr{'e}e grid, the side would appear as a
97 * straight line. For two consecutive vertices of the polygon,
98 * (lat_1, long1) and (lat_2,long_2), the line joining them (i.e., the
99 * polygon's side) is defined by:
100 *
101 \verbatim
102 lat_2 - lat_1
103 lat = lat_1 + (long - long_1) * ---------------
104 long_2 - long_1
105 \endverbatim
106 *
107 * where long_1 < long < long_2.
108 * The values of QbarA, etc., are determined by the integration of
109 * the Q function. Into www.integral-calculator.com, paste this
110 * expression :
111 *
112 \verbatim
113 sin(x)+ (2/3)e^2(sin(x))^3 + (3/5)e^4(sin(x))^5 + (4/7)e^6(sin(x))^7
114 \endverbatim
115 *
116 * and you'll get their values. (Last checked 30 Oct 2013).
117 *
118 * This function correctly computes (within the limits of the series
119 * approximation) the area of a quadrilateral on the ellipsoid when
120 * two of its sides run along meridians and the other two sides run
121 * along parallels of latitude.
122 *
123 * \param lon array of longitudes
124 * \param lat array of latitudes
125 * \param n number of lat,lon pairs
126 *
127 * \return area in square meters
128 */
129double G_ellipsoid_polygon_area(const double *lon, const double *lat, int n)
130{
131 double x1, y1, x2, y2, dx, dy;
132 double Qbar1, Qbar2;
133 double area;
134 double thresh = 1e-6; /* threshold for dy, should be between 1e-4 and 1e-7 */
135
136 x2 = Radians(lon[n - 1]);
137 y2 = Radians(lat[n - 1]);
138 Qbar2 = Qbar(y2);
139
140 area = 0.0;
141
142 while (--n >= 0) {
143 x1 = x2;
144 y1 = y2;
145 Qbar1 = Qbar2;
146
147 x2 = Radians(*lon++);
148 y2 = Radians(*lat++);
149 Qbar2 = Qbar(y2);
150
151 if (x1 > x2)
152 while (x1 - x2 > M_PI)
153 x2 += TWOPI;
154 else if (x2 > x1)
155 while (x2 - x1 > M_PI)
156 x1 += TWOPI;
157
158 dx = x2 - x1;
159 dy = y2 - y1;
160
161 if (fabs(dy) > thresh) {
162 /* account for different latitudes y1, y2 */
163 area += dx * (st->Qp - (Qbar2 - Qbar1) / dy);
164 /* original:
165 * area += dx * st->Qp - (dx / dy) * (Qbar2 - Qbar1);
166 */
167 }
168 else {
169 /* latitudes y1, y2 are (nearly) identical */
170 /* if y2 becomes similar to y1, i.e. y2 -> y1
171 * Qbar2 - Qbar1 -> 0 and dy -> 0
172 * (Qbar2 - Qbar1) / dy -> ?
173 * (Qbar2 - Qbar1) / dy should approach Q((y1 + y2) / 2)
174 * Metz 2017
175 */
176 area += dx * (st->Qp - Q((y1 + y2) / 2));
177 }
178 }
179 if ((area *= st->AE) < 0.0)
180 area = -area;
181
182 /* kludge - if polygon circles the south pole the area will be
183 * computed as if it cirlced the north pole. The correction is
184 * the difference between total surface area of the earth and
185 * the "north pole" area.
186 */
187 if (area > st->E)
188 area = st->E;
189 if (area > st->E / 2)
190 area = st->E - area;
191
192 return area;
193}
#define TWOPI
Definition: area_poly1.c:18
void G_begin_ellipsoid_polygon_area(double a, double e2)
Begin area calculations.
Definition: area_poly1.c:61
double G_ellipsoid_polygon_area(const double *lon, const double *lat, int n)
Area of lat-long polygon.
Definition: area_poly1.c:129
struct state state
Definition: parser.c:103
struct state * st
Definition: parser.c:104
#define Radians(x)
Definition: pi.h:6
#define x