GRASS Programmer's Manual 6.4.1(2011)
|
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 }