GRASS Programmer's Manual  6.5.svn(2014)-r66266
pole_in_poly.c
Go to the documentation of this file.
1 #include <grass/gis.h>
2
3 /**********************************************************
4  * G_pole_in_polygon(x, y, n)
5  * double *x, *y, n;
6  *
7  * For lat-lon coordinates, this routine determines if the polygon
8  * defined by the n verticies x,y contain one of the poles
9  *
10  * returns
11  * -1 if it contains the south pole,
12  * 1 if it contains the north pole,
13  * 0 no pole
14  *
15  * Note: don't use this routine if the projection isn't PROJECTION_LL
16  * no check is made by this routine for valid projection
17  ***********************************************************/
18
19 static int mystats(double, double, double, double, double *, double *);
20
21
38 int G_pole_in_polygon(const double *x, const double *y, int n)
39 {
40  int i;
41  double len, area, total_len, total_area;
42
43  if (n <= 1)
44  return 0;
45
46  mystats(x[n - 1], y[n - 1], x[0], y[0], &total_len, &total_area);
47  for (i = 1; i < n; i++) {
48  mystats(x[i - 1], y[i - 1], x[i], y[i], &len, &area);
49  total_len += len;
50  total_area += area;
51  }
52
53  /* if polygon contains a pole then the x-coordinate length of
54  * the perimeter should compute to 0, otherwise it should be about 360
55  * (or -360, depending on the direction of perimeter traversal)
56  *
57  * instead of checking for exactly 0, check from -1 to 1 to avoid
58  * roundoff error.
59  */
60  if (total_len < 1.0 && total_len > -1.0)
61  return 0;
62
64 }
65
66 static int mystats(double x0, double y0, double x1, double y1, double *len,
67  double *area)
68 {
69  if (x1 > x0)
70  while (x1 - x0 > 180)
71  x0 += 360;
72  else if (x0 > x1)
73  while (x0 - x1 > 180)
74  x0 -= 360;
75
76  *len = x0 - x1;
77
78  if (x0 > x1)
79  *area = (x0 - x1) * (y0 + y1) / 2.0;
80  else
81  *area = (x1 - x0) * (y1 + y0) / 2.0;
82
83  return 0;
84 }
int y
Definition: plot.c:34
int G_pole_in_polygon(const double *x, const double *y, int n)
pole in polygon
Definition: pole_in_poly.c:38
double x