Add area-weighted geographic polygon centroids.

This generalises the centre of mass formula given in:

  J. E. Brock, The Inertia Tensor for a Spherical Triangle,
  J. Applied Mechanics 42, 239 (1975)

In addition, centroids of zero-extent features now fall back to the
centroid of the feature treated as a lower dimension, i.e. polygons →
lines → points.

Specifically, if the centroid of a polygon or multipolygon feature is
ambiguous, it will treat it as a set of lines.  Similarly (but not
exactly the same), if a line or set of lines has zero length, it will be
treated as a set of points.

This also applies to d3.geo.path.centroid, except that the fallback
always occurs if the feature has zero extent (area or length).  For
geographic polygons, it was simpler to fall back if the centroid was
undefined, as area is not tracked in a simple fashion.

Fixes #1011.

Optimised d3_acos, and used d3_asin and d3_acos where appropriate
instead of Math.min and Math.max.
This commit is contained in:
Jason Davies 2013-01-15 22:50:42 +00:00
Родитель 915ad7f4c1
Коммит 492215ad42
10 изменённых файлов: 291 добавлений и 191 удалений

152
d3.js поставляемый
Просмотреть файл

@ -1298,12 +1298,12 @@ d3 = function() {
}
return d3_rgb(vv(h + 120), vv(h), vv(h - 120));
}
var π = Math.PI, ε = 1e-6, d3_radians = π / 180, d3_degrees = 180 / π;
var π = Math.PI, ε = 1e-6, ε2 = ε * ε, d3_radians = π / 180, d3_degrees = 180 / π;
function d3_sgn(x) {
return x > 0 ? 1 : x < 0 ? -1 : 0;
}
function d3_acos(x) {
return Math.acos(Math.max(-1, Math.min(1, x)));
return x > 1 ? 0 : x < -1 ? π : Math.acos(x);
}
function d3_asin(x) {
return x > 1 ? π / 2 : x < -1 ? -π / 2 : Math.asin(x);
@ -2170,7 +2170,7 @@ d3 = function() {
d[2] /= l;
}
function d3_geo_spherical(cartesian) {
return [ Math.atan2(cartesian[1], cartesian[0]), Math.asin(Math.max(-1, Math.min(1, cartesian[2]))) ];
return [ Math.atan2(cartesian[1], cartesian[0]), d3_asin(cartesian[2]) ];
}
function d3_geo_sphericalEqual(a, b) {
return Math.abs(a[0] - b[0]) < ε && Math.abs(a[1] - b[1]) < ε;
@ -2304,27 +2304,24 @@ d3 = function() {
};
}();
d3.geo.centroid = function(object) {
d3_geo_centroidDimension = d3_geo_centroidW = d3_geo_centroidX = d3_geo_centroidY = d3_geo_centroidZ = 0;
d3_geo_centroidW0 = d3_geo_centroidW1 = d3_geo_centroidX0 = d3_geo_centroidY0 = d3_geo_centroidZ0 = d3_geo_centroidX1 = d3_geo_centroidY1 = d3_geo_centroidZ1 = d3_geo_centroidX2 = d3_geo_centroidY2 = d3_geo_centroidZ2 = 0;
d3.geo.stream(object, d3_geo_centroid);
var m;
return d3_geo_centroidW && Math.abs(m = Math.sqrt(d3_geo_centroidX * d3_geo_centroidX + d3_geo_centroidY * d3_geo_centroidY + d3_geo_centroidZ * d3_geo_centroidZ)) > ε ? [ Math.atan2(d3_geo_centroidY, d3_geo_centroidX) * d3_degrees, Math.asin(Math.max(-1, Math.min(1, d3_geo_centroidZ / m))) * d3_degrees ] : [ NaN, NaN ];
var x = d3_geo_centroidX2, y = d3_geo_centroidY2, z = d3_geo_centroidZ2, m = x * x + y * y + z * z;
if (m < ε2) {
x = d3_geo_centroidX1, y = d3_geo_centroidY1, z = d3_geo_centroidZ1;
if (d3_geo_centroidW1 < ε) x = d3_geo_centroidX0, y = d3_geo_centroidY0, z = d3_geo_centroidZ0;
m = x * x + y * y + z * z;
if (m < ε2) return [ NaN, NaN ];
}
return [ Math.atan2(y, x) * d3_degrees, d3_asin(z / Math.sqrt(m)) * d3_degrees ];
};
var d3_geo_centroidDimension, d3_geo_centroidW, d3_geo_centroidX, d3_geo_centroidY, d3_geo_centroidZ;
var d3_geo_centroidW0, d3_geo_centroidW1, d3_geo_centroidX0, d3_geo_centroidY0, d3_geo_centroidZ0, d3_geo_centroidX1, d3_geo_centroidY1, d3_geo_centroidZ1, d3_geo_centroidX2, d3_geo_centroidY2, d3_geo_centroidZ2;
var d3_geo_centroid = {
sphere: function() {
if (d3_geo_centroidDimension < 2) {
d3_geo_centroidDimension = 2;
d3_geo_centroidW = d3_geo_centroidX = d3_geo_centroidY = d3_geo_centroidZ = 0;
}
},
sphere: d3_noop,
point: d3_geo_centroidPoint,
lineStart: d3_geo_centroidLineStart,
lineEnd: d3_geo_centroidLineEnd,
polygonStart: function() {
if (d3_geo_centroidDimension < 2) {
d3_geo_centroidDimension = 2;
d3_geo_centroidW = d3_geo_centroidX = d3_geo_centroidY = d3_geo_centroidZ = 0;
}
d3_geo_centroid.lineStart = d3_geo_centroidRingStart;
},
polygonEnd: function() {
@ -2332,36 +2329,18 @@ d3 = function() {
}
};
function d3_geo_centroidPoint(λ, φ) {
if (d3_geo_centroidDimension) return;
++d3_geo_centroidW;
λ *= d3_radians;
var cosφ = Math.cos(φ *= d3_radians);
d3_geo_centroidX += (cosφ * Math.cos(λ) - d3_geo_centroidX) / d3_geo_centroidW;
d3_geo_centroidY += (cosφ * Math.sin(λ) - d3_geo_centroidY) / d3_geo_centroidW;
d3_geo_centroidZ += (Math.sin(φ) - d3_geo_centroidZ) / d3_geo_centroidW;
d3_geo_centroidPointXYZ(cosφ * Math.cos(λ), cosφ * Math.sin(λ), Math.sin(φ));
}
function d3_geo_centroidRingStart() {
var λ00, φ00;
d3_geo_centroidDimension = 1;
d3_geo_centroidLineStart();
d3_geo_centroidDimension = 2;
var linePoint = d3_geo_centroid.point;
d3_geo_centroid.point = function(λ, φ) {
linePoint(λ00 = λ, φ00 = φ);
};
d3_geo_centroid.lineEnd = function() {
d3_geo_centroid.point(λ00, φ00);
d3_geo_centroidLineEnd();
d3_geo_centroid.lineEnd = d3_geo_centroidLineEnd;
};
function d3_geo_centroidPointXYZ(x, y, z) {
++d3_geo_centroidW0;
d3_geo_centroidX0 += (x - d3_geo_centroidX0) / d3_geo_centroidW0;
d3_geo_centroidY0 += (y - d3_geo_centroidY0) / d3_geo_centroidW0;
d3_geo_centroidZ0 += (z - d3_geo_centroidZ0) / d3_geo_centroidW0;
}
function d3_geo_centroidLineStart() {
var x0, y0, z0;
if (d3_geo_centroidDimension > 1) return;
if (d3_geo_centroidDimension < 1) {
d3_geo_centroidDimension = 1;
d3_geo_centroidW = d3_geo_centroidX = d3_geo_centroidY = d3_geo_centroidZ = 0;
}
d3_geo_centroid.point = function(λ, φ) {
λ *= d3_radians;
var cosφ = Math.cos(φ *= d3_radians);
@ -2369,19 +2348,51 @@ d3 = function() {
y0 = cosφ * Math.sin(λ);
z0 = Math.sin(φ);
d3_geo_centroid.point = nextPoint;
d3_geo_centroidPointXYZ(x0, y0, z0);
};
function nextPoint(λ, φ) {
λ *= d3_radians;
var cosφ = Math.cos(φ *= d3_radians), x = cosφ * Math.cos(λ), y = cosφ * Math.sin(λ), z = Math.sin(φ), w = Math.atan2(Math.sqrt((w = y0 * z - z0 * y) * w + (w = z0 * x - x0 * z) * w + (w = x0 * y - y0 * x) * w), x0 * x + y0 * y + z0 * z);
d3_geo_centroidW += w;
d3_geo_centroidX += w * (x0 + (x0 = x));
d3_geo_centroidY += w * (y0 + (y0 = y));
d3_geo_centroidZ += w * (z0 + (z0 = z));
d3_geo_centroidW1 += w;
d3_geo_centroidX1 += w * (x0 + (x0 = x));
d3_geo_centroidY1 += w * (y0 + (y0 = y));
d3_geo_centroidZ1 += w * (z0 + (z0 = z));
d3_geo_centroidPointXYZ(x0, y0, z0);
}
}
function d3_geo_centroidLineEnd() {
d3_geo_centroid.point = d3_geo_centroidPoint;
}
function d3_geo_centroidRingStart() {
var λ00, φ00, x0, y0, z0;
d3_geo_centroid.point = function(λ, φ) {
λ00 = λ, φ00 = φ;
d3_geo_centroid.point = nextPoint;
λ *= d3_radians;
var cosφ = Math.cos(φ *= d3_radians);
x0 = cosφ * Math.cos(λ);
y0 = cosφ * Math.sin(λ);
z0 = Math.sin(φ);
d3_geo_centroidPointXYZ(x0, y0, z0);
};
d3_geo_centroid.lineEnd = function() {
nextPoint(λ00, φ00);
d3_geo_centroid.lineEnd = d3_geo_centroidLineEnd;
d3_geo_centroid.point = d3_geo_centroidPoint;
};
function nextPoint(λ, φ) {
λ *= d3_radians;
var cosφ = Math.cos(φ *= d3_radians), x = cosφ * Math.cos(λ), y = cosφ * Math.sin(λ), z = Math.sin(φ), cx = y0 * z - z0 * y, cy = z0 * x - x0 * z, cz = x0 * y - y0 * x, m = Math.sqrt(cx * cx + cy * cy + cz * cz), u = x0 * x + y0 * y + z0 * z, v = m && -d3_acos(u) / m, w = Math.atan2(m, u);
d3_geo_centroidX2 += v * cx;
d3_geo_centroidY2 += v * cy;
d3_geo_centroidZ2 += v * cz;
d3_geo_centroidW1 += w;
d3_geo_centroidX1 += w * (x0 + (x0 = x));
d3_geo_centroidY1 += w * (y0 + (y0 = y));
d3_geo_centroidZ1 += w * (z0 + (z0 = z));
d3_geo_centroidPointXYZ(x0, y0, z0);
}
}
function d3_true() {
return true;
}
@ -3146,29 +3157,22 @@ d3 = function() {
}
};
function d3_geo_pathCentroidPoint(x, y) {
if (d3_geo_centroidDimension) return;
d3_geo_centroidX += x;
d3_geo_centroidY += y;
++d3_geo_centroidZ;
d3_geo_centroidX0 += x;
d3_geo_centroidY0 += y;
++d3_geo_centroidZ0;
}
function d3_geo_pathCentroidLineStart() {
var x0, y0;
if (d3_geo_centroidDimension !== 1) {
if (d3_geo_centroidDimension < 1) {
d3_geo_centroidDimension = 1;
d3_geo_centroidX = d3_geo_centroidY = d3_geo_centroidZ = 0;
} else return;
}
d3_geo_pathCentroid.point = function(x, y) {
d3_geo_pathCentroid.point = nextPoint;
x0 = x, y0 = y;
d3_geo_pathCentroidPoint(x0 = x, y0 = y);
};
function nextPoint(x, y) {
var dx = x - x0, dy = y - y0, z = Math.sqrt(dx * dx + dy * dy);
d3_geo_centroidX += z * (x0 + x) / 2;
d3_geo_centroidY += z * (y0 + y) / 2;
d3_geo_centroidZ += z;
x0 = x, y0 = y;
d3_geo_centroidX1 += z * (x0 + x) / 2;
d3_geo_centroidY1 += z * (y0 + y) / 2;
d3_geo_centroidZ1 += z;
d3_geo_pathCentroidPoint(x0 = x, y0 = y);
}
}
function d3_geo_pathCentroidLineEnd() {
@ -3176,20 +3180,20 @@ d3 = function() {
}
function d3_geo_pathCentroidRingStart() {
var x00, y00, x0, y0;
if (d3_geo_centroidDimension < 2) {
d3_geo_centroidDimension = 2;
d3_geo_centroidX = d3_geo_centroidY = d3_geo_centroidZ = 0;
}
d3_geo_pathCentroid.point = function(x, y) {
d3_geo_pathCentroid.point = nextPoint;
x00 = x0 = x, y00 = y0 = y;
d3_geo_pathCentroidPoint(x00 = x0 = x, y00 = y0 = y);
};
function nextPoint(x, y) {
var z = y0 * x - x0 * y;
d3_geo_centroidX += z * (x0 + x);
d3_geo_centroidY += z * (y0 + y);
d3_geo_centroidZ += z * 3;
x0 = x, y0 = y;
var dx = x - x0, dy = y - y0, z = Math.sqrt(dx * dx + dy * dy);
d3_geo_centroidX1 += z * (x0 + x) / 2;
d3_geo_centroidY1 += z * (y0 + y) / 2;
d3_geo_centroidZ1 += z;
z = y0 * x - x0 * y;
d3_geo_centroidX2 += z * (x0 + x);
d3_geo_centroidY2 += z * (y0 + y);
d3_geo_centroidZ2 += z * 3;
d3_geo_pathCentroidPoint(x0 = x, y0 = y);
}
d3_geo_pathCentroid.lineEnd = function() {
nextPoint(x00, y00);
@ -3328,9 +3332,9 @@ d3 = function() {
return d3_geo_pathAreaSum;
};
path.centroid = function(object) {
d3_geo_centroidDimension = d3_geo_centroidX = d3_geo_centroidY = d3_geo_centroidZ = 0;
d3_geo_centroidX0 = d3_geo_centroidY0 = d3_geo_centroidZ0 = d3_geo_centroidX1 = d3_geo_centroidY1 = d3_geo_centroidZ1 = d3_geo_centroidX2 = d3_geo_centroidY2 = d3_geo_centroidZ2 = 0;
d3.geo.stream(object, projectStream(d3_geo_pathCentroid));
return d3_geo_centroidZ ? [ d3_geo_centroidX / d3_geo_centroidZ, d3_geo_centroidY / d3_geo_centroidZ ] : undefined;
return d3_geo_centroidZ2 ? [ d3_geo_centroidX2 / d3_geo_centroidZ2, d3_geo_centroidY2 / d3_geo_centroidZ2 ] : d3_geo_centroidZ1 ? [ d3_geo_centroidX1 / d3_geo_centroidZ1, d3_geo_centroidY1 / d3_geo_centroidZ1 ] : d3_geo_centroidZ0 ? [ d3_geo_centroidX0 / d3_geo_centroidZ0, d3_geo_centroidY0 / d3_geo_centroidZ0 ] : [ NaN, NaN ];
};
path.bounds = function(object) {
d3_geo_pathBoundsX1 = d3_geo_pathBoundsY1 = -(d3_geo_pathBoundsX0 = d3_geo_pathBoundsY0 = Infinity);
@ -3531,11 +3535,11 @@ d3 = function() {
var cosδφ = Math.cos(δφ), sinδφ = Math.sin(δφ), cosδγ = Math.cos(δγ), sinδγ = Math.sin(δγ);
function rotation(λ, φ) {
var cosφ = Math.cos(φ), x = Math.cos(λ) * cosφ, y = Math.sin(λ) * cosφ, z = Math.sin(φ), k = z * cosδφ + x * sinδφ;
return [ Math.atan2(y * cosδγ - k * sinδγ, x * cosδφ - z * sinδφ), Math.asin(Math.max(-1, Math.min(1, k * cosδγ + y * sinδγ))) ];
return [ Math.atan2(y * cosδγ - k * sinδγ, x * cosδφ - z * sinδφ), d3_asin(k * cosδγ + y * sinδγ) ];
}
rotation.invert = function(λ, φ) {
var cosφ = Math.cos(φ), x = Math.cos(λ) * cosφ, y = Math.sin(λ) * cosφ, z = Math.sin(φ), k = z * cosδγ - y * sinδγ;
return [ Math.atan2(y * cosδγ + z * sinδγ, x * cosδφ + k * sinδφ), Math.asin(Math.max(-1, Math.min(1, k * cosδφ - x * sinδφ))) ];
return [ Math.atan2(y * cosδγ + z * sinδγ, x * cosδφ + k * sinδφ), d3_asin(k * cosδφ - x * sinδφ) ];
};
return rotation;
}

10
d3.min.js поставляемый

Различия файлов скрыты, потому что одна или несколько строк слишком длинны

Просмотреть файл

@ -1,37 +1,51 @@
import "../core/noop";
import "../math/trigonometry";
import "geo";
import "stream";
d3.geo.centroid = function(object) {
d3_geo_centroidDimension = d3_geo_centroidW = d3_geo_centroidX = d3_geo_centroidY = d3_geo_centroidZ = 0;
d3_geo_centroidW0 = d3_geo_centroidW1 =
d3_geo_centroidX0 = d3_geo_centroidY0 = d3_geo_centroidZ0 =
d3_geo_centroidX1 = d3_geo_centroidY1 = d3_geo_centroidZ1 =
d3_geo_centroidX2 = d3_geo_centroidY2 = d3_geo_centroidZ2 = 0;
d3.geo.stream(object, d3_geo_centroid);
var m;
return d3_geo_centroidW && Math.abs(m = Math.sqrt(d3_geo_centroidX * d3_geo_centroidX + d3_geo_centroidY * d3_geo_centroidY + d3_geo_centroidZ * d3_geo_centroidZ)) > ε
? [Math.atan2(d3_geo_centroidY, d3_geo_centroidX) * d3_degrees, Math.asin(Math.max(-1, Math.min(1, d3_geo_centroidZ / m))) * d3_degrees]
: [NaN, NaN];
var x = d3_geo_centroidX2,
y = d3_geo_centroidY2,
z = d3_geo_centroidZ2,
m = x * x + y * y + z * z;
// If the area-weighted centroid is undefined, fall back to length-weighted centroid.
if (m < ε2) {
x = d3_geo_centroidX1, y = d3_geo_centroidY1, z = d3_geo_centroidZ1;
// If the feature has zero length, fall back to arithmetic mean of point vectors.
if (d3_geo_centroidW1 < ε) x = d3_geo_centroidX0, y = d3_geo_centroidY0, z = d3_geo_centroidZ0;
m = x * x + y * y + z * z;
// If the feature still has an undefined centroid, then return.
if (m < ε2) return [NaN, NaN];
}
return [Math.atan2(y, x) * d3_degrees, d3_asin(z / Math.sqrt(m)) * d3_degrees];
};
var d3_geo_centroidDimension,
d3_geo_centroidW,
d3_geo_centroidX,
d3_geo_centroidY,
d3_geo_centroidZ;
var d3_geo_centroidW0,
d3_geo_centroidW1,
d3_geo_centroidX0,
d3_geo_centroidY0,
d3_geo_centroidZ0,
d3_geo_centroidX1,
d3_geo_centroidY1,
d3_geo_centroidZ1,
d3_geo_centroidX2,
d3_geo_centroidY2,
d3_geo_centroidZ2;
var d3_geo_centroid = {
sphere: function() {
if (d3_geo_centroidDimension < 2) {
d3_geo_centroidDimension = 2;
d3_geo_centroidW = d3_geo_centroidX = d3_geo_centroidY = d3_geo_centroidZ = 0;
}
},
sphere: d3_noop,
point: d3_geo_centroidPoint,
lineStart: d3_geo_centroidLineStart,
lineEnd: d3_geo_centroidLineEnd,
polygonStart: function() {
if (d3_geo_centroidDimension < 2) {
d3_geo_centroidDimension = 2;
d3_geo_centroidW = d3_geo_centroidX = d3_geo_centroidY = d3_geo_centroidZ = 0;
}
d3_geo_centroid.lineStart = d3_geo_centroidRingStart;
},
polygonEnd: function() {
@ -41,42 +55,21 @@ var d3_geo_centroid = {
// Arithmetic mean of Cartesian vectors.
function d3_geo_centroidPoint(λ, φ) {
if (d3_geo_centroidDimension) return;
++d3_geo_centroidW;
λ *= d3_radians;
var cosφ = Math.cos(φ *= d3_radians);
d3_geo_centroidX += (cosφ * Math.cos(λ) - d3_geo_centroidX) / d3_geo_centroidW;
d3_geo_centroidY += (cosφ * Math.sin(λ) - d3_geo_centroidY) / d3_geo_centroidW;
d3_geo_centroidZ += (Math.sin(φ) - d3_geo_centroidZ) / d3_geo_centroidW;
d3_geo_centroidPointXYZ(cosφ * Math.cos(λ), cosφ * Math.sin(λ), Math.sin(φ));
}
function d3_geo_centroidRingStart() {
var λ00, φ00; // first point
d3_geo_centroidDimension = 1;
d3_geo_centroidLineStart();
d3_geo_centroidDimension = 2;
var linePoint = d3_geo_centroid.point;
d3_geo_centroid.point = function(λ, φ) {
linePoint(λ00 = λ, φ00 = φ);
};
d3_geo_centroid.lineEnd = function() {
d3_geo_centroid.point(λ00, φ00);
d3_geo_centroidLineEnd();
d3_geo_centroid.lineEnd = d3_geo_centroidLineEnd;
};
function d3_geo_centroidPointXYZ(x, y, z) {
++d3_geo_centroidW0;
d3_geo_centroidX0 += (x - d3_geo_centroidX0) / d3_geo_centroidW0;
d3_geo_centroidY0 += (y - d3_geo_centroidY0) / d3_geo_centroidW0;
d3_geo_centroidZ0 += (z - d3_geo_centroidZ0) / d3_geo_centroidW0;
}
function d3_geo_centroidLineStart() {
var x0, y0, z0; // previous point
if (d3_geo_centroidDimension > 1) return;
if (d3_geo_centroidDimension < 1) {
d3_geo_centroidDimension = 1;
d3_geo_centroidW = d3_geo_centroidX = d3_geo_centroidY = d3_geo_centroidZ = 0;
}
d3_geo_centroid.point = function(λ, φ) {
λ *= d3_radians;
var cosφ = Math.cos(φ *= d3_radians);
@ -84,6 +77,7 @@ function d3_geo_centroidLineStart() {
y0 = cosφ * Math.sin(λ);
z0 = Math.sin(φ);
d3_geo_centroid.point = nextPoint;
d3_geo_centroidPointXYZ(x0, y0, z0);
};
function nextPoint(λ, φ) {
@ -95,13 +89,61 @@ function d3_geo_centroidLineStart() {
w = Math.atan2(
Math.sqrt((w = y0 * z - z0 * y) * w + (w = z0 * x - x0 * z) * w + (w = x0 * y - y0 * x) * w),
x0 * x + y0 * y + z0 * z);
d3_geo_centroidW += w;
d3_geo_centroidX += w * (x0 + (x0 = x));
d3_geo_centroidY += w * (y0 + (y0 = y));
d3_geo_centroidZ += w * (z0 + (z0 = z));
d3_geo_centroidW1 += w;
d3_geo_centroidX1 += w * (x0 + (x0 = x));
d3_geo_centroidY1 += w * (y0 + (y0 = y));
d3_geo_centroidZ1 += w * (z0 + (z0 = z));
d3_geo_centroidPointXYZ(x0, y0, z0);
}
}
function d3_geo_centroidLineEnd() {
d3_geo_centroid.point = d3_geo_centroidPoint;
}
// See J. E. Brock, The Inertia Tensor for a Spherical Triangle,
// J. Applied Mechanics 42, 239 (1975).
function d3_geo_centroidRingStart() {
var λ00, φ00, // first point
x0, y0, z0; // previous point
d3_geo_centroid.point = function(λ, φ) {
λ00 = λ, φ00 = φ;
d3_geo_centroid.point = nextPoint;
λ *= d3_radians;
var cosφ = Math.cos(φ *= d3_radians);
x0 = cosφ * Math.cos(λ);
y0 = cosφ * Math.sin(λ);
z0 = Math.sin(φ);
d3_geo_centroidPointXYZ(x0, y0, z0);
};
d3_geo_centroid.lineEnd = function() {
nextPoint(λ00, φ00);
d3_geo_centroid.lineEnd = d3_geo_centroidLineEnd;
d3_geo_centroid.point = d3_geo_centroidPoint;
};
function nextPoint(λ, φ) {
λ *= d3_radians;
var cosφ = Math.cos(φ *= d3_radians),
x = cosφ * Math.cos(λ),
y = cosφ * Math.sin(λ),
z = Math.sin(φ),
cx = y0 * z - z0 * y,
cy = z0 * x - x0 * z,
cz = x0 * y - y0 * x,
m = Math.sqrt(cx * cx + cy * cy + cz * cz),
u = x0 * x + y0 * y + z0 * z,
v = m && -d3_acos(u) / m, // area weight
w = Math.atan2(m, u); // line weight
d3_geo_centroidX2 += v * cx;
d3_geo_centroidY2 += v * cy;
d3_geo_centroidZ2 += v * cz;
d3_geo_centroidW1 += w;
d3_geo_centroidX1 += w * (x0 + (x0 = x));
d3_geo_centroidY1 += w * (y0 + (y0 = y));
d3_geo_centroidZ1 += w * (z0 + (z0 = z));
d3_geo_centroidPointXYZ(x0, y0, z0);
}
}

Просмотреть файл

@ -22,33 +22,25 @@ var d3_geo_pathCentroid = {
};
function d3_geo_pathCentroidPoint(x, y) {
if (d3_geo_centroidDimension) return;
d3_geo_centroidX += x;
d3_geo_centroidY += y;
++d3_geo_centroidZ;
d3_geo_centroidX0 += x;
d3_geo_centroidY0 += y;
++d3_geo_centroidZ0;
}
function d3_geo_pathCentroidLineStart() {
var x0, y0;
if (d3_geo_centroidDimension !== 1) {
if (d3_geo_centroidDimension < 1) {
d3_geo_centroidDimension = 1;
d3_geo_centroidX = d3_geo_centroidY = d3_geo_centroidZ = 0;
} else return;
}
d3_geo_pathCentroid.point = function(x, y) {
d3_geo_pathCentroid.point = nextPoint;
x0 = x, y0 = y;
d3_geo_pathCentroidPoint(x0 = x, y0 = y);
};
function nextPoint(x, y) {
var dx = x - x0, dy = y - y0, z = Math.sqrt(dx * dx + dy * dy);
d3_geo_centroidX += z * (x0 + x) / 2;
d3_geo_centroidY += z * (y0 + y) / 2;
d3_geo_centroidZ += z;
x0 = x, y0 = y;
d3_geo_centroidX1 += z * (x0 + x) / 2;
d3_geo_centroidY1 += z * (y0 + y) / 2;
d3_geo_centroidZ1 += z;
d3_geo_pathCentroidPoint(x0 = x, y0 = y);
}
}
@ -59,24 +51,24 @@ function d3_geo_pathCentroidLineEnd() {
function d3_geo_pathCentroidRingStart() {
var x00, y00, x0, y0;
if (d3_geo_centroidDimension < 2) {
d3_geo_centroidDimension = 2;
d3_geo_centroidX = d3_geo_centroidY = d3_geo_centroidZ = 0;
}
// For the first point, …
d3_geo_pathCentroid.point = function(x, y) {
d3_geo_pathCentroid.point = nextPoint;
x00 = x0 = x, y00 = y0 = y;
d3_geo_pathCentroidPoint(x00 = x0 = x, y00 = y0 = y);
};
// For subsequent points, …
function nextPoint(x, y) {
var z = y0 * x - x0 * y;
d3_geo_centroidX += z * (x0 + x);
d3_geo_centroidY += z * (y0 + y);
d3_geo_centroidZ += z * 3;
x0 = x, y0 = y;
var dx = x - x0, dy = y - y0, z = Math.sqrt(dx * dx + dy * dy);
d3_geo_centroidX1 += z * (x0 + x) / 2;
d3_geo_centroidY1 += z * (y0 + y) / 2;
d3_geo_centroidZ1 += z;
z = y0 * x - x0 * y;
d3_geo_centroidX2 += z * (x0 + x);
d3_geo_centroidY2 += z * (y0 + y);
d3_geo_centroidZ2 += z * 3;
d3_geo_pathCentroidPoint(x0 = x, y0 = y);
}
// For the last point, return to the start.

Просмотреть файл

@ -38,9 +38,14 @@ d3.geo.path = function() {
};
path.centroid = function(object) {
d3_geo_centroidDimension = d3_geo_centroidX = d3_geo_centroidY = d3_geo_centroidZ = 0;
d3_geo_centroidX0 = d3_geo_centroidY0 = d3_geo_centroidZ0 =
d3_geo_centroidX1 = d3_geo_centroidY1 = d3_geo_centroidZ1 =
d3_geo_centroidX2 = d3_geo_centroidY2 = d3_geo_centroidZ2 = 0;
d3.geo.stream(object, projectStream(d3_geo_pathCentroid));
return d3_geo_centroidZ ? [d3_geo_centroidX / d3_geo_centroidZ, d3_geo_centroidY / d3_geo_centroidZ] : undefined;
return d3_geo_centroidZ2 ? [d3_geo_centroidX2 / d3_geo_centroidZ2, d3_geo_centroidY2 / d3_geo_centroidZ2]
: d3_geo_centroidZ1 ? [d3_geo_centroidX1 / d3_geo_centroidZ1, d3_geo_centroidY1 / d3_geo_centroidZ1]
: d3_geo_centroidZ0 ? [d3_geo_centroidX0 / d3_geo_centroidZ0, d3_geo_centroidY0 / d3_geo_centroidZ0]
: [NaN, NaN];
};
path.bounds = function(object) {

Просмотреть файл

@ -52,7 +52,7 @@ function d3_geo_rotationφγ(δφ, δγ) {
k = z * cosδφ + x * sinδφ;
return [
Math.atan2(y * cosδγ - k * sinδγ, x * cosδφ - z * sinδφ),
Math.asin(Math.max(-1, Math.min(1, k * cosδγ + y * sinδγ)))
d3_asin(k * cosδγ + y * sinδγ)
];
}
@ -64,7 +64,7 @@ function d3_geo_rotationφγ(δφ, δγ) {
k = z * cosδγ - y * sinδγ;
return [
Math.atan2(y * cosδγ + z * sinδγ, x * cosδφ + k * sinδφ),
Math.asin(Math.max(-1, Math.min(1, k * cosδφ - x * sinδφ)))
d3_asin(k * cosδφ - x * sinδφ)
];
};

Просмотреть файл

@ -3,7 +3,7 @@ import "../math/trigonometry";
function d3_geo_spherical(cartesian) {
return [
Math.atan2(cartesian[1], cartesian[0]),
Math.asin(Math.max(-1, Math.min(1, cartesian[2])))
d3_asin(cartesian[2])
];
}

Просмотреть файл

@ -1,5 +1,6 @@
var π = Math.PI,
ε = 1e-6,
ε2 = ε * ε,
d3_radians = π / 180,
d3_degrees = 180 / π;
@ -8,7 +9,7 @@ function d3_sgn(x) {
}
function d3_acos(x) {
return Math.acos(Math.max(-1, Math.min(1, x)));
return x > 1 ? 0 : x < -1 ? π : Math.acos(x);
}
function d3_asin(x) {

Просмотреть файл

@ -52,28 +52,86 @@ suite.addBatch({
assert.inDelta(centroid({type: "MultiLineString", coordinates: [[[0, 0], [0, 2]]]}), [0, 1], 1e-6);
},
"an empty line is treated as a point": function(centroid) {
assert.inDelta(centroid({type: "LineString", coordinates: [[1, 1], [1, 1]]}), [1, 1], 1e-6); // TODO
assert.inDelta(centroid({type: "GeometryCollection", geometries: [{type: "Point", coordinates: [0, 0]}, {type: "LineString", coordinates: [[1, 2], [1, 2]]}]}), [0.499847, 1.000038], 1e-6);
"a line of zero length is treated as points": function(centroid) {
assert.inDelta(centroid({type: "LineString", coordinates: [[1, 1], [1, 1]]}), [1, 1], 1e-6);
assert.inDelta(centroid({type: "GeometryCollection", geometries: [{type: "Point", coordinates: [0, 0]}, {type: "LineString", coordinates: [[1, 2], [1, 2]]}]}), [0.666534, 1.333408], 1e-6);
},
"an empty polygon is treated as a point": function(centroid) {
assert.inDelta(centroid({type: "Polygon", coordinates: [[[1, 1], [1, 1], [1, 1], [1, 1]]]}), [1, 1], 1e-6); // TODO
assert.inDelta(centroid({type: "GeometryCollection", geometries: [{type: "Point", coordinates: [0, 0]}, {type: "Polygon", coordinates: [[[1, 2], [1, 2], [1, 2], [1, 2]]]}]}), [0.499847, 1.000038], 1e-6);
"an empty polygon with non-zero extent is treated as a line": function(centroid) {
assert.inDelta(centroid({type: "Polygon", coordinates: [[[1, 1], [2, 1], [3, 1], [2, 1], [1, 1]]]}), [2, 1.000076], 1e-6);
assert.inDelta(centroid({type: "GeometryCollection", geometries: [{type: "Point", coordinates: [0, 0]}, {type: "Polygon", coordinates: [[[1, 2], [1, 2], [1, 2], [1, 2]]]}]}), [0.799907, 1.600077], 1e-6);
},
"an empty polygon with zero extent is treated as a point": function(centroid) {
assert.inDelta(centroid({type: "Polygon", coordinates: [[[1, 1], [1, 1], [1, 1], [1, 1]]]}), [1, 1], 1e-6);
assert.inDelta(centroid({type: "GeometryCollection", geometries: [{type: "Point", coordinates: [0, 0]}, {type: "Polygon", coordinates: [[[1, 2], [1, 2], [1, 2], [1, 2]]]}]}), [0.799907, 1.600077], 1e-6);
},
"the centroid of the equator is ambiguous": function(centroid) {
assert.ok(centroid({type: "LineString", coordinates: [[0, 0], [120, 0], [-120, 0], [0, 0]]}).every(isNaN));
},
// TODO Dont treat a polygon as a line string.
"the centroid of a polygon is the (spherical) average of its surface": function(centroid) {
assert.inDelta(centroid({type: "Polygon", coordinates: [[[0, -90], [0, 0], [0, 90], [1, 0], [0, -90]]]}), [.5, 0], 1e-6);
assert.inDelta(centroid(_.geo.circle().angle(5).origin([0, 45])()), [0, 45], 1e-6);
assert.inDelta(centroid({type: "Polygon", coordinates: [_.range(-180, 180 + 1 / 2, 1).map(function(x) { return [x, -60]; })]})[1], -90, 1e-6);
assert.inDelta(centroid({type: "Polygon", coordinates: [[[0, -10], [0, 10], [10, 10], [10, -10], [0, -10]]]}), [5, 0], 1e-6);
},
// TODO Actually test multiple polygons here.
"the centroid of a set of polygons is the (spherical) average of its surface": function(centroid) {
assert.inDelta(centroid({type: "MultiPolygon", coordinates: [[[[0, -90], [0, 0], [0, 90], [1, 0], [0, -90]]]]}), [.5, 0], 1e-6);
assert.inDelta(centroid({type: "GeometryCollection", geometries: [{type: "Polygon", coordinates: [[[0, -90], [0, 0], [0, 90], [1, 0], [0, -90]]]}]}), [.5, 0], 1e-6);
var circle = _.geo.circle();
assert.inDelta(centroid({
type: "MultiPolygon",
coordinates: [
circle.angle(45).origin([0, 0])().coordinates,
circle.angle(60).origin([180, 0])().coordinates
]
}), [180, 0], 1e-6);
},
"the centroid of a lune is the (spherical) average of its surface": function(centroid) {
assert.inDelta(centroid({type: "Polygon", coordinates: [[[0, -90], [0, 0], [0, 90], [1, 0], [0, -90]]]}), [.5, 0], 1e-6);
},
"the centroid of a small circle is its origin": {
"5°": function(centroid) {
assert.inDelta(centroid(_.geo.circle().angle(5).origin([30, 45])()), [30, 45], 1e-6);
},
"135°": function(centroid) {
assert.inDelta(centroid(_.geo.circle().angle(135).origin([30, 45])()), [30, 45], 1e-6);
},
"South Pole": function(centroid) {
assert.equal(centroid({type: "Polygon", coordinates: [_.range(-180, 180 + 1 / 2, 1).map(function(x) { return [x, -60]; })]})[1], -90);
},
"equator": function(centroid) {
assert.inDelta(centroid({type: "Polygon", coordinates: [[[0, -10], [0, 10], [10, 10], [10, -10], [0, -10]]]}), [5, 0], 1e-6);
},
"equator with coincident points": function(centroid) {
assert.inDelta(centroid({type: "Polygon", coordinates: [[[0, -10], [0, 10], [0, 10], [10, 10], [10, -10], [0, -10]]]}), [5, 0], 1e-6);
},
"other": function(centroid) {
assert.inDelta(centroid({type: "Polygon", coordinates: [[[-180, 0], [-180, 10], [-179, 10], [-179, 0], [-180, 0]]]}), [-179.5, 4.987448], 1e-6);
},
"concentric rings": function(centroid) {
var circle = _.geo.circle().origin([0, 45]),
coordinates = circle.angle(60)().coordinates;
coordinates.push(circle.angle(45)().coordinates[0].reverse());
assert.inDelta(centroid({type: "Polygon", coordinates: coordinates}), [0, 45], 1e-6);
}
},
"the centroid of a spherical square on the equator": function(centroid) {
assert.inDelta(centroid({type: "Polygon", coordinates: [[[0, -10], [0, 10], [10, 10], [10, -10], [0, -10]]]}), [5, 0], 1e-6);
},
"the centroid of a spherical square touching the antimeridian": function(centroid) {
assert.inDelta(centroid({type: "Polygon", coordinates: [[[-180, 0], [-180, 10], [-179, 10], [-179, 0], [-180, 0]]]}), [-179.5, 4.987448], 1e-6);
},
"concentric rings": function(centroid) {
var circle = _.geo.circle().origin([0, 45]),
coordinates = circle.angle(60)().coordinates;
coordinates.push(circle.angle(45)().coordinates[0].reverse());
assert.inDelta(centroid({type: "Polygon", coordinates: coordinates}), [0, 45], 1e-6);
},
"the centroid of a sphere is ambiguous": function(centroid) {
@ -81,7 +139,7 @@ suite.addBatch({
},
"the centroid of a feature is the centroid of its constituent geometry": function(centroid) {
assert.inDelta(centroid({type: "Feature", geometry: {type: "LineString", coordinates: [[1, 1], [1, 1]]}}), [1, 1], 1e-6); // TODO
assert.inDelta(centroid({type: "Feature", geometry: {type: "LineString", coordinates: [[1, 1], [1, 1]]}}), [1, 1], 1e-6);
assert.inDelta(centroid({type: "Feature", geometry: {type: "Point", coordinates: [1, 1]}}), [1, 1], 1e-6);
assert.inDelta(centroid({type: "Feature", geometry: {type: "Polygon", coordinates: [[[0, -90], [0, 0], [0, 90], [1, 0], [0, -90]]]}}), [.5, 0], 1e-6);
},
@ -105,23 +163,23 @@ suite.addBatch({
{type: "Polygon", coordinates: [[[-180, 0], [-180, 1], [-179, 1], [-179, 0], [-180, 0]]]},
{type: "LineString", coordinates: [[179, 0], [180, 0]]},
{type: "Point", coordinates: [0, 0]}
]}), [-179.5, 0.5], 1e-6);
]}), [-179.5, 0.500006], 1e-6);
assert.inDelta(centroid({type: "GeometryCollection", geometries: [
{type: "Point", coordinates: [0, 0]},
{type: "LineString", coordinates: [[179, 0], [180, 0]]},
{type: "Polygon", coordinates: [[[-180, 0], [-180, 1], [-179, 1], [-179, 0], [-180, 0]]]}
]}), [-179.5, 0.5], 1e-6);
]}), [-179.5, 0.500006], 1e-6);
},
"the centroid of the sphere and a point only considers the sphere": function(centroid) {
assert.ok(centroid({type: "GeometryCollection", geometries: [
"the centroid of the sphere and a point is the point": function(centroid) {
assert.deepEqual(centroid({type: "GeometryCollection", geometries: [
{type: "Sphere"},
{type: "Point", coordinates: [0, 0]}
]}).every(isNaN));
assert.ok(centroid({type: "GeometryCollection", geometries: [
]}), [0, 0]);
assert.deepEqual(centroid({type: "GeometryCollection", geometries: [
{type: "Point", coordinates: [0, 0]},
{type: "Sphere"}
]}).every(isNaN));
]}), [0, 0]);
}
}
});

Просмотреть файл

@ -130,7 +130,7 @@ suite.addBatch({
assert.deepEqual(centroid({type: "Point", coordinates: [0, 0]}), [480, 250]);
},
"of an empty multipoint": function(centroid) {
assert.isUndefined(centroid({type: "MultiPoint", coordinates: []}));
assert.ok(centroid({type: "MultiPoint", coordinates: []}).every(isNaN));
},
"of a singleton multipoint": function(centroid) {
assert.deepEqual(centroid({type: "MultiPoint", coordinates: [[0, 0]]}), [480, 250]);
@ -139,15 +139,15 @@ suite.addBatch({
assert.deepEqual(centroid({type: "MultiPoint", coordinates: [[-122, 37], [-74, 40]]}), [-10, 57.5]);
},
"of an empty linestring": function(centroid) {
assert.isUndefined(centroid({type: "LineString", coordinates: []}));
assert.ok(centroid({type: "LineString", coordinates: []}).every(isNaN));
},
"of a linestring with two points": function(centroid) {
assert.deepEqual(centroid({type: "LineString", coordinates: [[100, 0], [0, 0]]}), [730, 250]);
assert.deepEqual(centroid({type: "LineString", coordinates: [[0, 0], [100, 0], [101, 0]]}), [732.5, 250]);
},
"of a linestring with two points, one unique": function(centroid) {
assert.isUndefined(centroid({type: "LineString", coordinates: [[-122, 37], [-122, 37]]}));
assert.isUndefined(centroid({type: "LineString", coordinates: [[ -74, 40], [ -74, 40]]}));
assert.deepEqual(centroid({type: "LineString", coordinates: [[-122, 37], [-122, 37]]}), [-130, 65]);
assert.deepEqual(centroid({type: "LineString", coordinates: [[-74, 40], [-74, 40]]}), [110, 50]);
},
"of a linestring with three points; two unique": function(centroid) {
assert.deepEqual(centroid({type: "LineString", coordinates: [[-122, 37], [-74, 40], [-74, 40]]}), [-10, 57.5]);
@ -162,7 +162,7 @@ suite.addBatch({
assert.deepEqual(centroid({type: "Polygon", coordinates: [[[100, 0], [100, 1], [101, 1], [101, 0], [100, 0]]]}), [982.5, 247.5]);
},
"of a zero-area polygon": function(centroid) {
assert.isUndefined(centroid({type: "Polygon", coordinates: [[[1, 0], [2, 0], [3, 0], [1, 0]]]}));
assert.deepEqual(centroid({type: "Polygon", coordinates: [[[1, 0], [2, 0], [3, 0], [1, 0]]]}), [490, 250]);
},
"of a polygon with two rings, one with zero area": function(centroid) {
assert.deepEqual(centroid({type: "Polygon", coordinates: [
@ -180,7 +180,7 @@ suite.addBatch({
}), [479.642857, 250], 1e-6);
},
"of an empty multipolygon": function(centroid) {
assert.isUndefined(centroid({type: "MultiPolygon", coordinates: []}));
assert.ok(centroid({type: "MultiPolygon", coordinates: []}).every(isNaN));
},
"of a singleton multipolygon": function(centroid) {
assert.deepEqual(centroid({type: "MultiPolygon", coordinates: [[[[100, 0], [100, 1], [101, 1], [101, 0], [100, 0]]]]}), [982.5, 247.5]);
@ -294,9 +294,7 @@ suite.addBatch({
[[109.378, 189.584], [797.758, 504.660]], 1e-3);
},
"centroid of a line string": function(p) {
var centroid = p.centroid({type: "LineString", coordinates: [[-122, 37], [-74, 40], [-100, 0]]});
assert.inDelta(centroid[0], 545.131, 1e-3);
assert.inDelta(centroid[1], 253.860, 1e-3);
assert.inDelta(p.centroid({type: "LineString", coordinates: [[-122, 37], [-74, 40], [-100, 0]]}), [545.131, 253.860], 1e-3);
}
},