-
-
Notifications
You must be signed in to change notification settings - Fork 302
/
pole_in_poly.c
79 lines (65 loc) · 1.98 KB
/
pole_in_poly.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
/*!
\file lib/gis/pole_in_poly.c
\brief GIS Library - Pole in polygon
(C) 2001-2009 by the GRASS Development Team
This program is free software under the GNU General Public License
(>=v2). Read the file COPYING that comes with GRASS for details.
\author CERL
*/
#include <grass/gis.h>
static void mystats(double, double, double, double, double *, double *);
/*!
* \brief Check if pole is in polygon
*
* For latitude-longitude coordinates, this routine determines if the polygon
* defined by the <i>n</i> coordinate vertices <i>x,y</i> contains one of the
* poles.
*
* <b>Note:</b> Use this routine only if the projection is PROJECTION_LL.
*
* \param x array of x coordinates
* \param y array of y coordinates
* \param n number of coordinates
*
* \return -1 if it contains the south pole
* \return 1 if it contains the north pole
* \return 0 if it contains neither pole.
*/
int G_pole_in_polygon(const double *x, const double *y, int n)
{
int i;
double len, area, total_len, total_area;
if (n <= 1)
return 0;
mystats(x[n - 1], y[n - 1], x[0], y[0], &total_len, &total_area);
for (i = 1; i < n; i++) {
mystats(x[i - 1], y[i - 1], x[i], y[i], &len, &area);
total_len += len;
total_area += area;
}
/* if polygon contains a pole then the x-coordinate length of
* the perimeter should compute to 0, otherwise it should be about 360
* (or -360, depending on the direction of perimeter traversal)
*
* instead of checking for exactly 0, check from -1 to 1 to avoid
* roundoff error.
*/
if (total_len < 1.0 && total_len > -1.0)
return 0;
return total_area >= 0.0 ? 1 : -1;
}
static void mystats(double x0, double y0, double x1, double y1, double *len,
double *area)
{
if (x1 > x0)
while (x1 - x0 > 180)
x0 += 360;
else if (x0 > x1)
while (x0 - x1 > 180)
x0 -= 360;
*len = x0 - x1;
if (x0 > x1)
*area = (x0 - x1) * (y0 + y1) / 2.0;
else
*area = (x1 - x0) * (y1 + y0) / 2.0;
}