Rewrite @turf/isolines (#2918)

* Reimplement @turf/isolines to clear up the licensing concerns around marchingsquares.

---------

Co-authored-by: James Beard <[email protected]>

Author: mfedderly

eslint.config.mjs

@@ -71,17 +71,5 @@ export default tsEslint.config(
       },
     },
   },
-  {
-    files: [
-      "packages/turf-isobands/lib/marchingsquares-isobands.js",
-      "packages/turf-isolines/lib/marchingsquares-isocontours.js",
-    ],
-
-    languageOptions: {
-      globals: {
-        ...globals.browser,
-      },
-    },
-  },
   prettierRecommended
 );

packages/turf-isobands/package.json

@@ -4,7 +4,8 @@
   "description": "Takes a grid of values (GeoJSON format) and a set of threshold ranges. It outputs polygons that group areas within those ranges, effectively creating filled contour isobands.",
   "author": "Turf Authors",
   "contributors": [
-    "Stefano Borghi <@stebogit>"
+    "Stefano Borghi <@stebogit>",
+    "Matt Fedderly <@mfedderly>"
   ],
   "license": "MIT",
   "bugs": {

packages/turf-isolines/index.ts

@@ -2,15 +2,14 @@ import { bbox } from "@turf/bbox";
 import { coordEach } from "@turf/meta";
 import { collectionOf } from "@turf/invariant";
 import { multiLineString, featureCollection, isObject } from "@turf/helpers";
-// @ts-expect-error Legacy JS library with no types defined
-import { isoContours } from "marchingsquares";
 import { gridToMatrix } from "./lib/grid-to-matrix.js";
 import {
   FeatureCollection,
   Point,
   MultiLineString,
   Feature,
   GeoJsonProperties,
+  Position,
 } from "geojson";
 
 /**
@@ -70,6 +69,10 @@ function isolines(
   // Isoline methods
   const matrix = gridToMatrix(pointGrid, { zProperty: zProperty, flip: true });
 
+  // A quick note on what 'top' and 'bottom' mean in coordinate system of `matrix`:
+  // Remember that the southern hemisphere is represented by negative numbers,
+  // so a matrix Y of 0 is actually the *bottom*, and a Y of dy - 1 is the *top*.
+
   // check that the resulting matrix has consistent x and y dimensions and
   // has at least a 2x2 size so that we can actually build grid squares
   const dx = matrix[0].length;
@@ -122,18 +125,226 @@ function createIsoLines(
 
     const properties = { ...commonProperties, ...breaksProperties[i] };
     properties[zProperty] = threshold;
-    // Pass options to marchingsquares lib to reproduce historical turf
-    // behaviour.
-    const isoline = multiLineString(
-      isoContours(matrix, threshold, { linearRing: false, noFrame: true }),
-      properties
-    );
+    const isoline = multiLineString(isoContours(matrix, threshold), properties);
 
     results.push(isoline);
   }
   return results;
 }
 
+function isoContours(
+  matrix: ReadonlyArray<ReadonlyArray<number>>,
+  threshold: number
+): Position[][] {
+  // see https://en.wikipedia.org/wiki/Marching_squares
+  const segments: [Position, Position][] = [];
+
+  const dy = matrix.length;
+  const dx = matrix[0].length;
+
+  for (let y = 0; y < dy - 1; y++) {
+    for (let x = 0; x < dx - 1; x++) {
+      const tr = matrix[y + 1][x + 1];
+      const br = matrix[y][x + 1];
+      const bl = matrix[y][x];
+      const tl = matrix[y + 1][x];
+
+      let grid =
+        (tl >= threshold ? 8 : 0) |
+        (tr >= threshold ? 4 : 0) |
+        (br >= threshold ? 2 : 0) |
+        (bl >= threshold ? 1 : 0);
+
+      switch (grid) {
+        case 0:
+          continue;
+        case 1:
+          segments.push([
+            [x + frac(bl, br), y],
+            [x, y + frac(bl, tl)],
+          ]);
+          break;
+        case 2:
+          segments.push([
+            [x + 1, y + frac(br, tr)],
+            [x + frac(bl, br), y],
+          ]);
+          break;
+        case 3:
+          segments.push([
+            [x + 1, y + frac(br, tr)],
+            [x, y + frac(bl, tl)],
+          ]);
+          break;
+        case 4:
+          segments.push([
+            [x + frac(tl, tr), y + 1],
+            [x + 1, y + frac(br, tr)],
+          ]);
+          break;
+        case 5: {
+          // use the average of the 4 corners to differentiate the saddle case and correctly honor the counter-clockwise winding
+          const avg = (tl + tr + br + bl) / 4;
+          const above = avg >= threshold;
+
+          if (above) {
+            segments.push(
+              [
+                [x + frac(tl, tr), y + 1],
+                [x, y + frac(bl, tl)],
+              ],
+              [
+                [x + frac(bl, br), y],
+                [x + 1, y + frac(br, tr)],
+              ]
+            );
+          } else {
+            segments.push(
+              [
+                [x + frac(tl, tr), y + 1],
+                [x + 1, y + frac(br, tr)],
+              ],
+              [
+                [x + frac(bl, br), y],
+                [x, y + frac(bl, tl)],
+              ]
+            );
+          }
+          break;
+        }
+        case 6:
+          segments.push([
+            [x + frac(tl, tr), y + 1],
+            [x + frac(bl, br), y],
+          ]);
+          break;
+        case 7:
+          segments.push([
+            [x + frac(tl, tr), y + 1],
+            [x, y + frac(bl, tl)],
+          ]);
+          break;
+        case 8:
+          segments.push([
+            [x, y + frac(bl, tl)],
+            [x + frac(tl, tr), y + 1],
+          ]);
+          break;
+        case 9:
+          segments.push([
+            [x + frac(bl, br), y],
+            [x + frac(tl, tr), y + 1],
+          ]);
+          break;
+        case 10: {
+          const avg = (tl + tr + br + bl) / 4;
+          const above = avg >= threshold;
+
+          if (above) {
+            segments.push(
+              [
+                [x, y + frac(bl, tl)],
+                [x + frac(bl, br), y],
+              ],
+              [
+                [x + 1, y + frac(br, tr)],
+                [x + frac(tl, tr), y + 1],
+              ]
+            );
+          } else {
+            segments.push(
+              [
+                [x, y + frac(bl, tl)],
+                [x + frac(tl, tr), y + 1],
+              ],
+              [
+                [x + 1, y + frac(br, tr)],
+                [x + frac(bl, br), y],
+              ]
+            );
+          }
+          break;
+        }
+        case 11:
+          segments.push([
+            [x + 1, y + frac(br, tr)],
+            [x + frac(tl, tr), y + 1],
+          ]);
+          break;
+        case 12:
+          segments.push([
+            [x, y + frac(bl, tl)],
+            [x + 1, y + frac(br, tr)],
+          ]);
+          break;
+        case 13:
+          segments.push([
+            [x + frac(bl, br), y],
+            [x + 1, y + frac(br, tr)],
+          ]);
+          break;
+        case 14:
+          segments.push([
+            [x, y + frac(bl, tl)],
+            [x + frac(bl, br), y],
+          ]);
+          break;
+        case 15:
+          // all above
+          continue;
+      }
+    }
+  }
+
+  const contours: Position[][] = [];
+
+  while (segments.length > 0) {
+    const contour: Position[] = [...segments.shift()!];
+    contours.push(contour);
+
+    let found: boolean;
+    do {
+      found = false;
+      for (let i = 0; i < segments.length; i++) {
+        const segment = segments[i];
+        // add the segment's end point to the end of the contour
+        if (
+          segment[0][0] === contour[contour.length - 1][0] &&
+          segment[0][1] === contour[contour.length - 1][1]
+        ) {
+          found = true;
+          contour.push(segment[1]);
+          segments.splice(i, 1);
+          break;
+        }
+        // add the segment's start point to the start of the contour
+        if (
+          segment[1][0] === contour[0][0] &&
+          segment[1][1] === contour[0][1]
+        ) {
+          found = true;
+          contour.unshift(segment[0]);
+          segments.splice(i, 1);
+          break;
+        }
+      }
+    } while (found);
+  }
+
+  return contours;
+
+  // get the linear interpolation fraction of how far z is between z0 and z1
+  // See https://github.com/fschutt/marching-squares/blob/master/src/lib.rs
+  function frac(z0: number, z1: number): number {
+    if (z0 === z1) {
+      return 0.5;
+    }
+
+    let t = (threshold - z0) / (z1 - z0);
+    return t > 1 ? 1 : t < 0 ? 0 : t;
+  }
+}
+
 /**
  * Translates and scales isolines
  *

packages/turf-isolines/package.json

@@ -4,7 +4,8 @@
   "description": "Generate contour lines from a grid of data.",
   "author": "Turf Authors",
   "contributors": [
-    "Stefano Borghi <@stebogit>"
+    "Stefano Borghi <@stebogit>",
+    "Matt Fedderly <@mfedderly>"
   ],
   "license": "MIT",
   "bugs": {
@@ -79,7 +80,6 @@
     "@turf/invariant": "workspace:*",
     "@turf/meta": "workspace:*",
     "@types/geojson": "^7946.0.10",
-    "marchingsquares": "^1.3.3",
     "tslib": "^2.8.1"
   }
 }