Skip to content

Commit b46a74b

Browse files
mfedderlymfedderly
committed
Reimplement @turf/isolines to clear up the licensing concerns around marchingsquares.
1 parent 05494d8 commit b46a74b

File tree

8 files changed

+1004
-983
lines changed

8 files changed

+1004
-983
lines changed

‎eslint.config.mjs‎

Lines changed: 0 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -71,17 +71,5 @@ export default tsEslint.config(
7171
},
7272
},
7373
},
74-
{
75-
files: [
76-
"packages/turf-isobands/lib/marchingsquares-isobands.js",
77-
"packages/turf-isolines/lib/marchingsquares-isocontours.js",
78-
],
79-
80-
languageOptions: {
81-
globals: {
82-
...globals.browser,
83-
},
84-
},
85-
},
8674
prettierRecommended
8775
);

‎packages/turf-isolines/index.ts‎

Lines changed: 220 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -2,15 +2,14 @@ import { bbox } from "@turf/bbox";
22
import { coordEach } from "@turf/meta";
33
import { collectionOf } from "@turf/invariant";
44
import { multiLineString, featureCollection, isObject } from "@turf/helpers";
5-
// @ts-expect-error Legacy JS library with no types defined
6-
import { isoContours } from "marchingsquares";
75
import { gridToMatrix } from "./lib/grid-to-matrix.js";
86
import {
97
FeatureCollection,
108
Point,
119
MultiLineString,
1210
Feature,
1311
GeoJsonProperties,
12+
Position,
1413
} from "geojson";
1514

1615
/**
@@ -69,6 +68,11 @@ function isolines(
6968

7069
// Isoline methods
7170
const matrix = gridToMatrix(pointGrid, { zProperty: zProperty, flip: true });
71+
72+
// A quick note on what 'top' and 'bottom' mean in coordinate system of `matrix`:
73+
// Remember that the southern hemisphere is represented by negative numbers,
74+
// so a matrix Y of 0 is actually the *bottom*, and a Y of dy - 1 is the *top*.
75+
7276
const createdIsoLines = createIsoLines(
7377
matrix,
7478
breaks,
@@ -109,18 +113,226 @@ function createIsoLines(
109113

110114
const properties = { ...commonProperties, ...breaksProperties[i] };
111115
properties[zProperty] = threshold;
112-
// Pass options to marchingsquares lib to reproduce historical turf
113-
// behaviour.
114-
const isoline = multiLineString(
115-
isoContours(matrix, threshold, { linearRing: false, noFrame: true }),
116-
properties
117-
);
116+
const isoline = multiLineString(isoContours(matrix, threshold), properties);
118117

119118
results.push(isoline);
120119
}
121120
return results;
122121
}
123122

123+
function isoContours(
124+
matrix: ReadonlyArray<ReadonlyArray<number>>,
125+
threshold: number
126+
): Position[][] {
127+
// see https://en.wikipedia.org/wiki/Marching_squares
128+
const segments: [Position, Position][] = [];
129+
130+
const dy = matrix.length;
131+
const dx = matrix[0].length;
132+
133+
for (let y = 0; y < dy - 1; y++) {
134+
for (let x = 0; x < dx - 1; x++) {
135+
const tr = matrix[y + 1][x + 1];
136+
const br = matrix[y][x + 1];
137+
const bl = matrix[y][x];
138+
const tl = matrix[y + 1][x];
139+
140+
let grid =
141+
(tl >= threshold ? 8 : 0) |
142+
(tr >= threshold ? 4 : 0) |
143+
(br >= threshold ? 2 : 0) |
144+
(bl >= threshold ? 1 : 0);
145+
146+
switch (grid) {
147+
case 0:
148+
continue;
149+
case 1:
150+
segments.push([
151+
[x + frac(bl, br), y],
152+
[x, y + frac(bl, tl)],
153+
]);
154+
break;
155+
case 2:
156+
segments.push([
157+
[x + 1, y + frac(br, tr)],
158+
[x + frac(bl, br), y],
159+
]);
160+
break;
161+
case 3:
162+
segments.push([
163+
[x + 1, y + frac(br, tr)],
164+
[x, y + frac(bl, tl)],
165+
]);
166+
break;
167+
case 4:
168+
segments.push([
169+
[x + frac(tl, tr), y + 1],
170+
[x + 1, y + frac(br, tr)],
171+
]);
172+
break;
173+
case 5: {
174+
// use the average of the 4 corners to differentiate the saddle case and correctly honor the counter-clockwise winding
175+
const avg = (tl + tr + br + bl) / 4;
176+
const above = avg >= threshold;
177+
178+
if (above) {
179+
segments.push(
180+
[
181+
[x + frac(tl, tr), y + 1],
182+
[x, y + frac(bl, tl)],
183+
],
184+
[
185+
[x + frac(bl, br), y],
186+
[x + 1, y + frac(br, tr)],
187+
]
188+
);
189+
} else {
190+
segments.push(
191+
[
192+
[x + frac(tl, tr), y + 1],
193+
[x + 1, y + frac(br, tr)],
194+
],
195+
[
196+
[x + frac(bl, br), y],
197+
[x, y + frac(bl, tl)],
198+
]
199+
);
200+
}
201+
break;
202+
}
203+
case 6:
204+
segments.push([
205+
[x + frac(tl, tr), y + 1],
206+
[x + frac(bl, br), y],
207+
]);
208+
break;
209+
case 7:
210+
segments.push([
211+
[x + frac(tl, tr), y + 1],
212+
[x, y + frac(bl, tl)],
213+
]);
214+
break;
215+
case 8:
216+
segments.push([
217+
[x, y + frac(bl, tl)],
218+
[x + frac(tl, tr), y + 1],
219+
]);
220+
break;
221+
case 9:
222+
segments.push([
223+
[x + frac(bl, br), y],
224+
[x + frac(tl, tr), y + 1],
225+
]);
226+
break;
227+
case 10: {
228+
const avg = (tl + tr + br + bl) / 4;
229+
const above = avg >= threshold;
230+
231+
if (above) {
232+
segments.push(
233+
[
234+
[x, y + frac(bl, tl)],
235+
[x + frac(bl, br), y],
236+
],
237+
[
238+
[x + 1, y + frac(br, tr)],
239+
[x + frac(tl, tr), y + 1],
240+
]
241+
);
242+
} else {
243+
segments.push(
244+
[
245+
[x, y + frac(bl, tl)],
246+
[x + frac(tl, tr), y + 1],
247+
],
248+
[
249+
[x + 1, y + frac(br, tr)],
250+
[x + frac(bl, br), y],
251+
]
252+
);
253+
}
254+
break;
255+
}
256+
case 11:
257+
segments.push([
258+
[x + 1, y + frac(br, tr)],
259+
[x + frac(tl, tr), y + 1],
260+
]);
261+
break;
262+
case 12:
263+
segments.push([
264+
[x, y + frac(bl, tl)],
265+
[x + 1, y + frac(br, tr)],
266+
]);
267+
break;
268+
case 13:
269+
segments.push([
270+
[x + frac(bl, br), y],
271+
[x + 1, y + frac(br, tr)],
272+
]);
273+
break;
274+
case 14:
275+
segments.push([
276+
[x, y + frac(bl, tl)],
277+
[x + frac(bl, br), y],
278+
]);
279+
break;
280+
case 15:
281+
// all above
282+
continue;
283+
}
284+
}
285+
}
286+
287+
const contours: Position[][] = [];
288+
289+
while (segments.length > 0) {
290+
const contour: Position[] = [...segments.shift()!];
291+
contours.push(contour);
292+
293+
let found: boolean;
294+
do {
295+
found = false;
296+
for (let i = 0; i < segments.length; i++) {
297+
const segment = segments[i];
298+
// add the segment's end point to the end of the contour
299+
if (
300+
segment[0][0] === contour[contour.length - 1][0] &&
301+
segment[0][1] === contour[contour.length - 1][1]
302+
) {
303+
found = true;
304+
contour.push(segment[1]);
305+
segments.splice(i, 1);
306+
break;
307+
}
308+
// add the segment's start point to the start of the contour
309+
if (
310+
segment[1][0] === contour[0][0] &&
311+
segment[1][1] === contour[0][1]
312+
) {
313+
found = true;
314+
contour.unshift(segment[0]);
315+
segments.splice(i, 1);
316+
break;
317+
}
318+
}
319+
} while (found);
320+
}
321+
322+
return contours;
323+
324+
// get the linear interpolation fraction of how far z is between z0 and z1
325+
// See https://github.com/fschutt/marching-squares/blob/master/src/lib.rs
326+
function frac(z0: number, z1: number): number {
327+
if (z0 === z1) {
328+
return 0.5;
329+
}
330+
331+
let t = (threshold - z0) / (z1 - z0);
332+
return t > 1 ? 1 : t < 0 ? 0 : t;
333+
}
334+
}
335+
124336
/**
125337
* Translates and scales isolines
126338
*

‎packages/turf-isolines/package.json‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -79,7 +79,6 @@
7979
"@turf/invariant": "workspace:*",
8080
"@turf/meta": "workspace:*",
8181
"@types/geojson": "^7946.0.10",
82-
"marchingsquares": "^1.3.3",
8382
"tslib": "^2.8.1"
8483
}
8584
}

0 commit comments

Comments
 (0)