GRASS Programmer's Manual  6.4.3(2013)-r
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Macros Pages
area_poly1.c
Go to the documentation of this file.
1 
17 #include <math.h>
18 #include <grass/gis.h>
19 #include "pi.h"
20 
21 
22 #define TWOPI M_PI + M_PI
23 
24 static double QA, QB, QC;
25 static double QbarA, QbarB, QbarC, QbarD;
26 
27 static double AE;
29 static double Qp;
31 static double E;
34 static double Q(double x)
35 {
36  double sinx, sinx2;
37 
38  sinx = sin(x);
39  sinx2 = sinx * sinx;
40 
41  return sinx * (1 + sinx2 * (QA + sinx2 * (QB + sinx2 * QC)));
42 }
43 
44 static double Qbar(double x)
45 {
46  double cosx, cosx2;
47 
48  cosx = cos(x);
49  cosx2 = cosx * cosx;
50 
51  return cosx * (QbarA + cosx2 * (QbarB + cosx2 * (QbarC + cosx2 * QbarD)));
52 }
53 
54 
67 int G_begin_ellipsoid_polygon_area(double a, double e2)
68 {
69  double e4, e6;
70 
71  e4 = e2 * e2;
72  e6 = e4 * e2;
73 
74  AE = a * a * (1 - e2);
75 
76  QA = (2.0 / 3.0) * e2;
77  QB = (3.0 / 5.0) * e4;
78  QC = (4.0 / 7.0) * e6;
79 
80  QbarA = -1.0 - (2.0 / 3.0) * e2 - (3.0 / 5.0) * e4 - (4.0 / 7.0) * e6;
81  QbarB = (2.0 / 9.0) * e2 + (2.0 / 5.0) * e4 + (4.0 / 7.0) * e6;
82  QbarC = -(3.0 / 25.0) * e4 - (12.0 / 35.0) * e6;
83  QbarD = (4.0 / 49.0) * e6;
84 
85  Qp = Q(M_PI_2);
86  E = 4 * M_PI * Qp * AE;
87  if (E < 0.0)
88  E = -E;
89 
90  return 0;
91 }
92 
93 
110 double G_ellipsoid_polygon_area(const double *lon, const double *lat, int n)
111 {
112  double x1, y1, x2, y2, dx, dy;
113  double Qbar1, Qbar2;
114  double area;
115 
116  x2 = Radians(lon[n - 1]);
117  y2 = Radians(lat[n - 1]);
118  Qbar2 = Qbar(y2);
119 
120  area = 0.0;
121 
122  while (--n >= 0) {
123  x1 = x2;
124  y1 = y2;
125  Qbar1 = Qbar2;
126 
127  x2 = Radians(*lon++);
128  y2 = Radians(*lat++);
129  Qbar2 = Qbar(y2);
130 
131  if (x1 > x2)
132  while (x1 - x2 > M_PI)
133  x2 += TWOPI;
134  else if (x2 > x1)
135  while (x2 - x1 > M_PI)
136  x1 += TWOPI;
137 
138  dx = x2 - x1;
139  area += dx * (Qp - Q(y2));
140 
141  if ((dy = y2 - y1) != 0.0)
142  area += dx * Q(y2) - (dx / dy) * (Qbar2 - Qbar1);
143  }
144  if ((area *= AE) < 0.0)
145  area = -area;
146 
147  /* kludge - if polygon circles the south pole the area will be
148  * computed as if it cirlced the north pole. The correction is
149  * the difference between total surface area of the earth and
150  * the "north pole" area.
151  */
152  if (area > E)
153  area = E;
154  if (area > E / 2)
155  area = E - area;
156 
157  return area;
158 }
#define Radians(x)
Definition: pi.h:6
int G_begin_ellipsoid_polygon_area(double a, double e2)
Begin area calculations.
Definition: area_poly1.c:67
double G_ellipsoid_polygon_area(const double *lon, const double *lat, int n)
Area of lat-long polygon.
Definition: area_poly1.c:110
#define TWOPI
Definition: area_poly1.c:22
int n
Definition: dataquad.c:291