Skip to content

Commit aaec1a7

Browse files
mfedderlymfedderly
committed
Simplest turf-isolines rewrite, update tests to check out the geojson diffs
1 parent 9abb1e1 commit aaec1a7

File tree

5 files changed

+1262
-1190
lines changed

5 files changed

+1262
-1190
lines changed

‎packages/turf-isolines/index.ts‎

Lines changed: 230 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -108,17 +108,243 @@ function createIsoLines(
108108

109109
const properties = { ...commonProperties, ...breaksProperties[i] };
110110
properties[zProperty] = threshold;
111-
// Pass options to marchingsquares lib to reproduce historical turf
112-
// behaviour.
113111
const isoline = multiLineString(isoContours(matrix, threshold), properties);
114112

115113
results.push(isoline);
116114
}
117115
return results;
118116
}
119117

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

124350
/**

0 commit comments

Comments
 (0)