GRASS Programmer's Manual 6.4.1(2011)
area_poly1.c
Go to the documentation of this file.
00001 
00017 #include <math.h>
00018 #include <grass/gis.h>
00019 #include "pi.h"
00020 
00021 
00022 #define TWOPI M_PI + M_PI
00023 
00024 static double QA, QB, QC;
00025 static double QbarA, QbarB, QbarC, QbarD;
00026 
00027 static double AE;  
00029 static double Qp;  
00031 static double E;   
00034 static double Q(double x)
00035 {
00036     double sinx, sinx2;
00037 
00038     sinx = sin(x);
00039     sinx2 = sinx * sinx;
00040 
00041     return sinx * (1 + sinx2 * (QA + sinx2 * (QB + sinx2 * QC)));
00042 }
00043 
00044 static double Qbar(double x)
00045 {
00046     double cosx, cosx2;
00047 
00048     cosx = cos(x);
00049     cosx2 = cosx * cosx;
00050 
00051     return cosx * (QbarA + cosx2 * (QbarB + cosx2 * (QbarC + cosx2 * QbarD)));
00052 }
00053 
00054 
00067 int G_begin_ellipsoid_polygon_area(double a, double e2)
00068 {
00069     double e4, e6;
00070 
00071     e4 = e2 * e2;
00072     e6 = e4 * e2;
00073 
00074     AE = a * a * (1 - e2);
00075 
00076     QA = (2.0 / 3.0) * e2;
00077     QB = (3.0 / 5.0) * e4;
00078     QC = (4.0 / 7.0) * e6;
00079 
00080     QbarA = -1.0 - (2.0 / 3.0) * e2 - (3.0 / 5.0) * e4 - (4.0 / 7.0) * e6;
00081     QbarB = (2.0 / 9.0) * e2 + (2.0 / 5.0) * e4 + (4.0 / 7.0) * e6;
00082     QbarC = -(3.0 / 25.0) * e4 - (12.0 / 35.0) * e6;
00083     QbarD = (4.0 / 49.0) * e6;
00084 
00085     Qp = Q(M_PI_2);
00086     E = 4 * M_PI * Qp * AE;
00087     if (E < 0.0)
00088         E = -E;
00089 
00090     return 0;
00091 }
00092 
00093 
00110 double G_ellipsoid_polygon_area(const double *lon, const double *lat, int n)
00111 {
00112     double x1, y1, x2, y2, dx, dy;
00113     double Qbar1, Qbar2;
00114     double area;
00115 
00116     x2 = Radians(lon[n - 1]);
00117     y2 = Radians(lat[n - 1]);
00118     Qbar2 = Qbar(y2);
00119 
00120     area = 0.0;
00121 
00122     while (--n >= 0) {
00123         x1 = x2;
00124         y1 = y2;
00125         Qbar1 = Qbar2;
00126 
00127         x2 = Radians(*lon++);
00128         y2 = Radians(*lat++);
00129         Qbar2 = Qbar(y2);
00130 
00131         if (x1 > x2)
00132             while (x1 - x2 > M_PI)
00133                 x2 += TWOPI;
00134         else if (x2 > x1)
00135             while (x2 - x1 > M_PI)
00136                 x1 += TWOPI;
00137 
00138         dx = x2 - x1;
00139         area += dx * (Qp - Q(y2));
00140 
00141         if ((dy = y2 - y1) != 0.0)
00142             area += dx * Q(y2) - (dx / dy) * (Qbar2 - Qbar1);
00143     }
00144     if ((area *= AE) < 0.0)
00145         area = -area;
00146 
00147     /* kludge - if polygon circles the south pole the area will be
00148      * computed as if it cirlced the north pole. The correction is
00149      * the difference between total surface area of the earth and
00150      * the "north pole" area.
00151      */
00152     if (area > E)
00153         area = E;
00154     if (area > E / 2)
00155         area = E - area;
00156 
00157     return area;
00158 }
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Defines