Horizon / node_modules /d3-geo /src /polygonContains.js

Upload 6270 files

4083225 verified 6 months ago

2.79 kB

import {Adder} from "d3-array";
import {cartesian, cartesianCross, cartesianNormalizeInPlace} from "./cartesian.js";
import {abs, asin, atan2, cos, epsilon, epsilon2, halfPi, pi, quarterPi, sign, sin, tau} from "./math.js";

function longitude(point) {
  return abs(point[0]) <= pi ? point[0] : sign(point[0]) * ((abs(point[0]) + pi) % tau - pi);
}

export default function(polygon, point) {
  var lambda = longitude(point),
      phi = point[1],
      sinPhi = sin(phi),
      normal = [sin(lambda), -cos(lambda), 0],
      angle = 0,
      winding = 0;

  var sum = new Adder();

  if (sinPhi === 1) phi = halfPi + epsilon;
  else if (sinPhi === -1) phi = -halfPi - epsilon;

  for (var i = 0, n = polygon.length; i < n; ++i) {
    if (!(m = (ring = polygon[i]).length)) continue;
    var ring,
        m,
        point0 = ring[m - 1],
        lambda0 = longitude(point0),
        phi0 = point0[1] / 2 + quarterPi,
        sinPhi0 = sin(phi0),
        cosPhi0 = cos(phi0);

    for (var j = 0; j < m; ++j, lambda0 = lambda1, sinPhi0 = sinPhi1, cosPhi0 = cosPhi1, point0 = point1) {
      var point1 = ring[j],
          lambda1 = longitude(point1),
          phi1 = point1[1] / 2 + quarterPi,
          sinPhi1 = sin(phi1),
          cosPhi1 = cos(phi1),
          delta = lambda1 - lambda0,
          sign = delta >= 0 ? 1 : -1,
          absDelta = sign * delta,
          antimeridian = absDelta > pi,
          k = sinPhi0 * sinPhi1;

      sum.add(atan2(k * sign * sin(absDelta), cosPhi0 * cosPhi1 + k * cos(absDelta)));
      angle += antimeridian ? delta + sign * tau : delta;

      // Are the longitudes either side of the point’s meridian (lambda),
      // and are the latitudes smaller than the parallel (phi)?
      if (antimeridian ^ lambda0 >= lambda ^ lambda1 >= lambda) {
        var arc = cartesianCross(cartesian(point0), cartesian(point1));
        cartesianNormalizeInPlace(arc);
        var intersection = cartesianCross(normal, arc);
        cartesianNormalizeInPlace(intersection);
        var phiArc = (antimeridian ^ delta >= 0 ? -1 : 1) * asin(intersection[2]);
        if (phi > phiArc || phi === phiArc && (arc[0] || arc[1])) {
          winding += antimeridian ^ delta >= 0 ? 1 : -1;
        }
      }
    }
  }

  // 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 lambda
  // from the point to the South pole.  If it is zero, then the point is the
  // same side as the South pole.

  return (angle < -epsilon || angle < epsilon && sum < -epsilon2) ^ (winding & 1);
}

	import {Adder} from "d3-array";
	import {cartesian, cartesianCross, cartesianNormalizeInPlace} from "./cartesian.js";
	import {abs, asin, atan2, cos, epsilon, epsilon2, halfPi, pi, quarterPi, sign, sin, tau} from "./math.js";

	function longitude(point) {
	return abs(point[0]) <= pi ? point[0] : sign(point[0]) * ((abs(point[0]) + pi) % tau - pi);
	}

	export default function(polygon, point) {
	var lambda = longitude(point),
	phi = point[1],
	sinPhi = sin(phi),
	normal = [sin(lambda), -cos(lambda), 0],
	angle = 0,
	winding = 0;

	var sum = new Adder();

	if (sinPhi === 1) phi = halfPi + epsilon;
	else if (sinPhi === -1) phi = -halfPi - epsilon;

	for (var i = 0, n = polygon.length; i < n; ++i) {
	if (!(m = (ring = polygon[i]).length)) continue;
	var ring,
	m,
	point0 = ring[m - 1],
	lambda0 = longitude(point0),
	phi0 = point0[1] / 2 + quarterPi,
	sinPhi0 = sin(phi0),
	cosPhi0 = cos(phi0);

	for (var j = 0; j < m; ++j, lambda0 = lambda1, sinPhi0 = sinPhi1, cosPhi0 = cosPhi1, point0 = point1) {
	var point1 = ring[j],
	lambda1 = longitude(point1),
	phi1 = point1[1] / 2 + quarterPi,
	sinPhi1 = sin(phi1),
	cosPhi1 = cos(phi1),
	delta = lambda1 - lambda0,
	sign = delta >= 0 ? 1 : -1,
	absDelta = sign * delta,
	antimeridian = absDelta > pi,
	k = sinPhi0 * sinPhi1;

	sum.add(atan2(k * sign * sin(absDelta), cosPhi0 * cosPhi1 + k * cos(absDelta)));
	angle += antimeridian ? delta + sign * tau : delta;

	// Are the longitudes either side of the point’s meridian (lambda),
	// and are the latitudes smaller than the parallel (phi)?
	if (antimeridian ^ lambda0 >= lambda ^ lambda1 >= lambda) {
	var arc = cartesianCross(cartesian(point0), cartesian(point1));
	cartesianNormalizeInPlace(arc);
	var intersection = cartesianCross(normal, arc);
	cartesianNormalizeInPlace(intersection);
	var phiArc = (antimeridian ^ delta >= 0 ? -1 : 1) * asin(intersection[2]);
	if (phi > phiArc \|\| phi === phiArc && (arc[0] \|\| arc[1])) {
	winding += antimeridian ^ delta >= 0 ? 1 : -1;
	}
	}
	}
	}

	// 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 lambda
	// from the point to the South pole. If it is zero, then the point is the
	// same side as the South pole.

	return (angle < -epsilon \|\| angle < epsilon && sum < -epsilon2) ^ (winding & 1);
	}