GRASS Programmer's Manual
6.4.1(2011)
|
00001 #include <grass/gis.h> 00002 00003 /********************************************************** 00004 * G_pole_in_polygon(x, y, n) 00005 * double *x, *y, n; 00006 * 00007 * For lat-lon coordinates, this routine determines if the polygon 00008 * defined by the n verticies x,y contain one of the poles 00009 * 00010 * returns 00011 * -1 if it contains the south pole, 00012 * 1 if it contains the north pole, 00013 * 0 no pole 00014 * 00015 * Note: don't use this routine if the projection isn't PROJECTION_LL 00016 * no check is made by this routine for valid projection 00017 ***********************************************************/ 00018 00019 static int mystats(double, double, double, double, double *, double *); 00020 00021 00038 int G_pole_in_polygon(const double *x, const double *y, int n) 00039 { 00040 int i; 00041 double len, area, total_len, total_area; 00042 00043 if (n <= 1) 00044 return 0; 00045 00046 mystats(x[n - 1], y[n - 1], x[0], y[0], &total_len, &total_area); 00047 for (i = 1; i < n; i++) { 00048 mystats(x[i - 1], y[i - 1], x[i], y[i], &len, &area); 00049 total_len += len; 00050 total_area += area; 00051 } 00052 00053 /* if polygon contains a pole then the x-coordinate length of 00054 * the perimeter should compute to 0, otherwise it should be about 360 00055 * (or -360, depending on the direction of perimeter traversal) 00056 * 00057 * instead of checking for exactly 0, check from -1 to 1 to avoid 00058 * roundoff error. 00059 */ 00060 if (total_len < 1.0 && total_len > -1.0) 00061 return 0; 00062 00063 return total_area >= 0.0 ? 1 : -1; 00064 } 00065 00066 static int mystats(double x0, double y0, double x1, double y1, double *len, 00067 double *area) 00068 { 00069 if (x1 > x0) 00070 while (x1 - x0 > 180) 00071 x0 += 360; 00072 else if (x0 > x1) 00073 while (x0 - x1 > 180) 00074 x0 -= 360; 00075 00076 *len = x0 - x1; 00077 00078 if (x0 > x1) 00079 *area = (x0 - x1) * (y0 + y1) / 2.0; 00080 else 00081 *area = (x1 - x0) * (y1 + y0) / 2.0; 00082 00083 return 0; 00084 }