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