GRASS GIS 7 Programmer's Manual  7.0.3(2016)-r00000
area_poly1.c
Go to the documentation of this file.
1 
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  double AE;
24  double Qp;
25  double E;
26 } state;
27 
28 static struct state *st = &state;
29 
30 static 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 
40 static 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 
61 void 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 
129 double 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 
135  x2 = Radians(lon[n - 1]);
136  y2 = Radians(lat[n - 1]);
137  Qbar2 = Qbar(y2);
138 
139  area = 0.0;
140 
141  while (--n >= 0) {
142  x1 = x2;
143  y1 = y2;
144  Qbar1 = Qbar2;
145 
146  x2 = Radians(*lon++);
147  y2 = Radians(*lat++);
148  Qbar2 = Qbar(y2);
149 
150  if (x1 > x2)
151  while (x1 - x2 > M_PI)
152  x2 += TWOPI;
153  else if (x2 > x1)
154  while (x2 - x1 > M_PI)
155  x1 += TWOPI;
156 
157  dx = x2 - x1;
158  area += dx * (st->Qp - Q(y2));
159 
160  if ((dy = y2 - y1) != 0.0)
161  area += dx * Q(y2) - (dx / dy) * (Qbar2 - Qbar1);
162  }
163  if ((area *= st->AE) < 0.0)
164  area = -area;
165 
166  /* kludge - if polygon circles the south pole the area will be
167  * computed as if it cirlced the north pole. The correction is
168  * the difference between total surface area of the earth and
169  * the "north pole" area.
170  */
171  if (area > st->E)
172  area = st->E;
173  if (area > st->E / 2)
174  area = st->E - area;
175 
176  return area;
177 }
struct state * st
Definition: parser.c:101
#define Radians(x)
Definition: pi.h:6
double G_ellipsoid_polygon_area(const double *lon, const double *lat, int n)
Area of lat-long polygon.
Definition: area_poly1.c:129
#define TWOPI
Definition: area_poly1.c:18
struct state state
Definition: parser.c:100
void G_begin_ellipsoid_polygon_area(double a, double e2)
Begin area calculations.
Definition: area_poly1.c:61