-
-
Notifications
You must be signed in to change notification settings - Fork 294
/
area_poly1.c
193 lines (166 loc) · 4.98 KB
/
area_poly1.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
/*!
* \file lib/gis/area_poly1.c
*
* \brief GIS Library - Polygon area calculation routines.
*
* (C) 2001-2013 by the GRASS Development Team
*
* This program is free software under the GNU General Public License
* (>=v2). Read the file COPYING that comes with GRASS for details.
*
* \author Original author CERL
*/
#include <math.h>
#include <grass/gis.h>
#include "pi.h"
#define TWOPI M_PI + M_PI
static struct state {
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 */
} state;
static struct state *st = &state;
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)));
}
/*!
* \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)
{
double e4, e6;
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;
}
/*!
* \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 x1, y1, x2, y2, dx, dy;
double Qbar1, Qbar2;
double area;
double thresh = 1e-6; /* threshold for dy, should be between 1e-4 and 1e-7 */
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 cirlced 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;
return area;
}