-
Notifications
You must be signed in to change notification settings - Fork 0
/
point-in-polygon.js
72 lines (65 loc) · 2.51 KB
/
point-in-polygon.js
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
import "geo";
import "area";
import "cartesian";
import "../math/abs";
import "../math/trigonometry";
function d3_geo_pointInPolygon(point, polygon) {
if(d3.geo.area(polygon) ==0 ) return false; //a point cannot be in a zero sized polygon
var meridian = point[0],
parallel = point[1],
meridianNormal = [Math.sin(meridian), -Math.cos(meridian), 0],
polarAngle = 0,
winding = 0;
d3_geo_areaRingSum.reset();
for (var i = 0, n = polygon.length; i < n; ++i) {
var ring = polygon[i],
m = ring.length;
if (!m) continue;
var point0 = ring[0],
λ0 = point0[0],
φ0 = point0[1] / 2 + π / 4,
sinφ0 = Math.sin(φ0),
cosφ0 = Math.cos(φ0),
j = 1;
while (true) {
if (j === m) j = 0;
point = ring[j];
var λ = point[0],
φ = point[1] / 2 + π / 4,
sinφ = Math.sin(φ),
cosφ = Math.cos(φ),
dλ = λ - λ0,
sdλ = dλ >= 0 ? 1 : -1,
adλ = sdλ * dλ,
antimeridian = adλ > π,
k = sinφ0 * sinφ;
d3_geo_areaRingSum.add(Math.atan2(k * sdλ * Math.sin(adλ), cosφ0 * cosφ + k * Math.cos(adλ)));
polarAngle += antimeridian ? dλ + sdλ * τ : dλ;
// Are the longitudes either side of the point's meridian, and are the
// latitudes smaller than the parallel?
if (antimeridian ^ λ0 >= meridian ^ λ >= meridian) {
var arc = d3_geo_cartesianCross(d3_geo_cartesian(point0), d3_geo_cartesian(point));
d3_geo_cartesianNormalize(arc);
var intersection = d3_geo_cartesianCross(meridianNormal, arc);
d3_geo_cartesianNormalize(intersection);
var φarc = (antimeridian ^ dλ >= 0 ? -1 : 1) * d3_asin(intersection[2]);
if (parallel > φarc || parallel === φarc && (arc[0] || arc[1])) {
winding += antimeridian ^ dλ >= 0 ? 1 : -1;
}
}
if (!j++) break;
λ0 = λ, sinφ0 = sinφ, cosφ0 = cosφ, point0 = point;
}
}
// First, determine whether the South pole is inside or outside:
//
// It is inside if:
// * the polygon winds around it in a clockwise direction.
// * the polygon does not (cumulatively) wind around it, but has a negative
// (counter-clockwise) area.
//
// Second, count the (signed) number of times a segment crosses a meridian
// from the point to the South pole. If it is zero, then the point is the
// same side as the South pole.
return (polarAngle < -ε || polarAngle < ε && d3_geo_areaRingSum < 0) ^ (winding & 1);
}