fix(center-of-mass): calculate polygon centroids ring by ring

packages/turf-center-of-mass/index.ts

@@ -2,8 +2,7 @@ import { Feature, GeoJsonProperties, Point, Position } from "geojson";
 import { convex } from "@turf/convex";
 import { centroid } from "@turf/centroid";
 import { point } from "@turf/helpers";
-import { getType, getCoord } from "@turf/invariant";
-import { coordEach } from "@turf/meta";
+import { getType, getCoord, getCoords } from "@turf/invariant";
 
 /**
  * Takes any {@link Feature} or a {@link FeatureCollection} and returns its [center of mass](https://en.wikipedia.org/wiki/Center_of_mass) using this formula: [Centroid of Polygon](https://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon).
@@ -31,10 +30,7 @@ function centerOfMass<P extends GeoJsonProperties = GeoJsonProperties>(
     case "Point":
       return point(getCoord(geojson), options.properties);
     case "Polygon":
-      var coords: Position[] = [];
-      coordEach(geojson, function (coord) {
-        coords.push(coord);
-      });
+      var coords: Position[][] = getCoords(geojson);
 
       // First, we neutralize the feature (set it around coordinates [0,0]) to prevent rounding errors
       // We take any point to translate all the points around 0
@@ -45,31 +41,42 @@ function centerOfMass<P extends GeoJsonProperties = GeoJsonProperties>(
       var sArea = 0;
       var i, pi, pj, xi, xj, yi, yj, a;
 
-      var neutralizedPoints = coords.map(function (point) {
-        return [point[0] - translation[0], point[1] - translation[1]];
-      });
+      coords.forEach(function (ring, ringIndex) {
+        var ringArea = 0;
+        var ringSx = 0;
+        var ringSy = 0;
+        var neutralizedPoints = ring.map(function (point) {
+          return [point[0] - translation[0], point[1] - translation[1]];
+        });
 
-      for (i = 0; i < coords.length - 1; i++) {
-        // pi is the current point
-        pi = neutralizedPoints[i];
-        xi = pi[0];
-        yi = pi[1];
+        for (i = 0; i < ring.length - 1; i++) {
+          // pi is the current point
+          pi = neutralizedPoints[i];
+          xi = pi[0];
+          yi = pi[1];
 
-        // pj is the next point (pi+1)
-        pj = neutralizedPoints[i + 1];
-        xj = pj[0];
-        yj = pj[1];
+          // pj is the next point (pi+1)
+          pj = neutralizedPoints[i + 1];
+          xj = pj[0];
+          yj = pj[1];
 
-        // a is the common factor to compute the signed area and the final coordinates
-        a = xi * yj - xj * yi;
+          // a is the common factor to compute the signed area and the final coordinates
+          a = xi * yj - xj * yi;
 
-        // sArea is the sum used to compute the signed area
-        sArea += a;
+          // ringArea is the sum used to compute the signed area
+          ringArea += a;
 
-        // sx and sy are the sums used to compute the final coordinates
-        sx += (xi + xj) * a;
-        sy += (yi + yj) * a;
-      }
+          // ringSx and ringSy are the sums used to compute the final coordinates
+          ringSx += (xi + xj) * a;
+          ringSy += (yi + yj) * a;
+        }
+
+        // Add the exterior ring and subtract holes, regardless of winding order.
+        var sign = (ringIndex === 0 ? 1 : -1) * (ringArea < 0 ? -1 : 1);
+        sArea += sign * ringArea;
+        sx += sign * ringSx;
+        sy += sign * ringSy;
+      });
 
       // Shape has no area: fallback on turf.centroid
       if (sArea === 0) {

packages/turf-center-of-mass/package.json

@@ -3,6 +3,9 @@
   "version": "7.4.0",
   "description": "Finds the “balance point” of irregular shapes, like a country outline.",
   "author": "Turf Authors",
+  "contributors": [
+    "vera-rubin <@vera-rubin>"
+  ],
   "license": "MIT",
   "bugs": {
     "url": "https://github.com/Turfjs/turf/issues"

packages/turf-center-of-mass/test.ts

@@ -81,6 +81,41 @@ test("center of mass -- no area", (t) => {
   t.end();
 });
 
+test("center of mass -- polygon with a hole", (t) => {
+  const outer = [
+    [0, 0],
+    [10, 0],
+    [10, 10],
+    [0, 10],
+    [0, 0],
+  ];
+  const hole = [
+    [6, 2],
+    [6, 4],
+    [8, 4],
+    [8, 2],
+    [6, 2],
+  ];
+  const expected = point([59 / 12, 61 / 12]);
+
+  for (const reverseOuter of [false, true]) {
+    for (const reverseHole of [false, true]) {
+      const poly = polygon([
+        reverseOuter ? outer.slice().reverse() : outer,
+        reverseHole ? hole.slice().reverse() : hole,
+      ]);
+      const winding = `outer reversed: ${reverseOuter}, hole reversed: ${reverseHole}`;
+      t.deepEqual(centerOfMass(poly), expected, `feature, ${winding}`);
+      t.deepEqual(
+        centerOfMass(poly.geometry),
+        expected,
+        `geometry, ${winding}`
+      );
+    }
+  }
+  t.end();
+});
+
 test("center of mass -- point", (t) => {
   const p = point([0, 0]);
   const pt = centerOfMass(p);