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