import "geo"; import "area"; import "cartesian"; import "../math/abs"; import "../math/trigonometry"; function d3_geo_pointInPolygon(point, 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 < -ε) ^ (winding & 1); }