GRASS Programmer's Manual  6.4.1(2011)
pole_in_poly.c
Go to the documentation of this file.
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 }
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Defines