GRASS Programmer's Manual  6.5.svn(2014)-r66266
area_poly1.c
Go to the documentation of this file.
1
16 #include <math.h>
17 #include <grass/gis.h>
18 #include "pi.h"
19
20 #define TWOPI M_PI + M_PI
21
22 static double QA, QB, QC;
23 static double QbarA, QbarB, QbarC, QbarD;
24
25 static double AE;
27 static double Qp;
29 static double E;
32 static double Q(double x)
33 {
34  double sinx, sinx2;
35
36  sinx = sin(x);
37  sinx2 = sinx * sinx;
38
39  return sinx * (1 + sinx2 * (QA + sinx2 * (QB + sinx2 * QC)));
40 }
41
42 static double Qbar(double x)
43 {
44  double cosx, cosx2;
45
46  cosx = cos(x);
47  cosx2 = cosx * cosx;
48
49  return cosx * (QbarA + cosx2 * (QbarB + cosx2 * (QbarC + cosx2 * QbarD)));
50 }
51
52
65 int G_begin_ellipsoid_polygon_area(double a, double e2)
66 {
67  double e4, e6;
68
69  e4 = e2 * e2;
70  e6 = e4 * e2;
71
72  AE = a * a * (1 - e2);
73
74  QA = (2.0 / 3.0) * e2;
75  QB = (3.0 / 5.0) * e4;
76  QC = (4.0 / 7.0) * e6;
77
78  QbarA = -1.0 - (2.0 / 3.0) * e2 - (3.0 / 5.0) * e4 - (4.0 / 7.0) * e6;
79  QbarB = (2.0 / 9.0) * e2 + (2.0 / 5.0) * e4 + (4.0 / 7.0) * e6;
80  QbarC = -(3.0 / 25.0) * e4 - (12.0 / 35.0) * e6;
81  QbarD = (4.0 / 49.0) * e6;
82
83  Qp = Q(M_PI_2);
84  E = 4 * M_PI * Qp * AE;
85  if (E < 0.0)
86  E = -E;
87
88  return 0;
89 }
90
91
137 double G_ellipsoid_polygon_area(const double *lon, const double *lat, int n)
138 {
139  double x1, y1, x2, y2, dx, dy;
140  double Qbar1, Qbar2;
141  double area;
142
143  x2 = Radians(lon[n - 1]);
144  y2 = Radians(lat[n - 1]);
145  Qbar2 = Qbar(y2);
146
147  area = 0.0;
148
149  while (--n >= 0) {
150  x1 = x2;
151  y1 = y2;
152  Qbar1 = Qbar2;
153
156  Qbar2 = Qbar(y2);
157
158  if (x1 > x2)
159  while (x1 - x2 > M_PI)
160  x2 += TWOPI;
161  else if (x2 > x1)
162  while (x2 - x1 > M_PI)
163  x1 += TWOPI;
164
165  dx = x2 - x1;
166  area += dx * (Qp - Q(y2));
167
168  if ((dy = y2 - y1) != 0.0)
169  area += dx * Q(y2) - (dx / dy) * (Qbar2 - Qbar1);
170  }
171  if ((area *= AE) < 0.0)
172  area = -area;
173
174  /* kludge - if polygon circles the south pole the area will be
175  * computed as if it cirlced the north pole. The correction is
176  * the difference between total surface area of the earth and
177  * the "north pole" area.
178  */
179  if (area > E)
180  area = E;
181  if (area > E / 2)
182  area = E - area;
183
184  return area;
185 }
Definition: pi.h:6
int G_begin_ellipsoid_polygon_area(double a, double e2)
Begin area calculations.
Definition: area_poly1.c:65
double G_ellipsoid_polygon_area(const double *lon, const double *lat, int n)
Area of lat-long polygon.
Definition: area_poly1.c:137
#define TWOPI
Definition: area_poly1.c:20
int n