contour_isobands/
lib.rs

1#![cfg_attr(docsrs, feature(doc_auto_cfg))]
2//! Compute isobands *(i.e. contour polygons which enclose
3//! all the points of a grid included between two chosen values)*
4//! by applying marching squares to an array of values.
5//!
6//! Use the [`ContourBuilder`] to create the contour polygons.
7//! Output is a Vec of [`Band`], each [`Band`] is characterized by its minimum value,
8//! its maximum value, and a [`MultiPolygon`] geometry. Each band can be serialised
9//! to a GeoJSON Feature (using the `Band::to_geojson` method that is enabled if compiling
10//! this crate with the optional `geojson` feature flag).
11//! #### Example:
12#![cfg_attr(feature = "geojson", doc = "```")]
13#![cfg_attr(not(feature = "geojson"), doc = "```ignore")]
14//! # use contour_isobands::ContourBuilder;
15//! let c = ContourBuilder::new(10, 10); // x dim., y dim.
16//! let res = c.contours(&vec![
17//!     0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
18//!     0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
19//!     0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
20//!     0., 0., 0., 1., 1., 1., 0., 0., 0., 0.,
21//!     0., 0., 0., 1., 1., 1., 0., 0., 0., 0.,
22//!     0., 0., 0., 1., 1., 1., 0., 0., 0., 0.,
23//!     0., 0., 0., 1., 1., 1., 0., 0., 0., 0.,
24//!     0., 0., 0., 1., 1., 1., 0., 0., 0., 0.,
25//!     0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
26//!     0., 0., 0., 0., 0., 0., 0., 0., 0., 0.
27//! ], &[0.5, 1.]).unwrap(); // values, thresholds
28//!
29//! let expected_output = serde_json::json!({
30//!   "type": "Feature",
31//!   "geometry": {
32//!     "type": "MultiPolygon",
33//!     "coordinates": [[[
34//!         [3.0, 2.5], [2.5, 3.0], [2.5, 4.0], [2.5, 5.0],
35//!         [2.5, 6.0], [2.5, 7.0], [3.0, 7.5], [4.0, 7.5],
36//!         [5.0, 7.5], [5.5, 7.0], [5.5, 6.0], [5.5, 5.0],
37//!         [5.5, 4.0], [5.5, 3.0], [5.0, 2.5], [4.0, 2.5],
38//!         [3.0, 2.5]
39//!     ]]]
40//!   },
41//!   "properties": {"min_v": 0.5, "max_v": 1.0}
42//! });
43//!
44//! assert_eq!(
45//!     res[0].to_geojson(),
46//!     std::convert::TryFrom::try_from(expected_output).unwrap(),
47//! );
48//! ```
49//! [`MultiPolygon`]: ../geo_types/geometry/struct.MultiPolygon.html
50#![cfg_attr(debug_assertions, allow(dead_code))]
51mod area;
52mod contains;
53mod errors;
54mod grid;
55mod isobands;
56mod polygons;
57mod quadtree;
58mod shape_coordinates;
59mod utils;
60
61pub use crate::isobands::{isobands, Band, BandRaw, ContourBuilder};
62
63#[cfg(test)]
64mod tests {
65    use crate::isobands::isobands;
66    use geo_types::Point;
67
68    fn make_grid_from2d_vec(data: &[Vec<f64>]) -> (Vec<f64>, usize, usize) {
69        let width = data[0].len();
70        let height = data.len();
71        let mut grid: Vec<f64> = Vec::with_capacity(width * height);
72        for row in data {
73            grid.extend_from_slice(row);
74        }
75        (grid, width, height)
76    }
77
78    #[test]
79    fn isobands_err_matrix_empty() {
80        let matrix: Vec<Vec<f64>> = vec![vec![]];
81        let (matrix, width, height) = make_grid_from2d_vec(&matrix);
82        let lower_band = 1.;
83        let bandwidth = 2.;
84
85        let res = isobands(
86            &matrix,
87            &vec![lower_band, lower_band + bandwidth],
88            false,
89            width,
90            height,
91            false,
92        );
93        assert!(res.is_err());
94    }
95
96    #[test]
97    fn isobands_err_threshold_too_short() {
98        let matrix = vec![vec![1., 1.], vec![1., 5.]];
99        let (matrix, width, height) = make_grid_from2d_vec(&matrix);
100
101        let res = isobands(&matrix, &vec![2.], false, width, height, false);
102        assert!(res.is_err());
103    }
104
105    #[test]
106    fn isobands_err_matrix_rows_not_same_length() {
107        let matrix: Vec<Vec<f64>> = vec![vec![1., 1.], vec![1., 5., 5.]];
108        let (matrix, width, height) = make_grid_from2d_vec(&matrix);
109
110        let res = isobands(&matrix, &vec![1., 3.], false, width, height, false);
111        assert!(res.is_err());
112    }
113
114    #[test]
115    fn isobands_minimal() {
116        let matrix = vec![vec![1., 1.], vec![1., 5.]];
117        let (matrix, width, height) = make_grid_from2d_vec(&matrix);
118
119        let lower_band = 1.;
120        let bandwidth = 2.;
121
122        let res = isobands(
123            &matrix,
124            &vec![lower_band, lower_band + bandwidth],
125            false,
126            width,
127            height,
128            false,
129        )
130        .unwrap();
131        assert_eq!(
132            res[0].0,
133            vec![vec![
134                Point::new(0.5, 1.),
135                Point::new(1., 0.5),
136                Point::new(1., 0.),
137                Point::new(0., 0.),
138                Point::new(0., 1.),
139                Point::new(0.5, 1.),
140            ]]
141        );
142        assert_eq!(res[0].1, lower_band);
143        assert_eq!(res[0].2, lower_band + bandwidth);
144    }
145
146    #[test]
147    fn isoband_simple() {
148        let matrix = vec![
149            vec![1., 1., 1., 0.],
150            vec![1., 5., 5., 1.],
151            vec![0., 1., 1., 1.],
152        ];
153        let lower_band = 1.;
154        let bandwidth = 1.;
155        let (matrix, width, height) = make_grid_from2d_vec(&matrix);
156
157        let res = isobands(
158            &matrix,
159            &vec![lower_band, lower_band + bandwidth],
160            false,
161            width,
162            height,
163            false,
164        )
165        .unwrap();
166        assert_eq!(
167            res[0].0,
168            vec![
169                vec![
170                    Point::new(0.25, 1.),
171                    Point::new(1., 0.25),
172                    Point::new(2., 0.25),
173                    Point::new(2.75, 1.),
174                    Point::new(2., 1.75),
175                    Point::new(1., 1.75),
176                    Point::new(0.25, 1.),
177                ],
178                vec![
179                    Point::new(0., 1.),
180                    Point::new(1., 2.),
181                    Point::new(1., 2.),
182                    Point::new(2., 2.),
183                    Point::new(3., 2.),
184                    Point::new(3., 1.),
185                    Point::new(3., 1.),
186                    Point::new(2., 0.),
187                    Point::new(2., 0.),
188                    Point::new(1., 0.),
189                    Point::new(0., 0.),
190                    Point::new(0., 1.),
191                ]
192            ]
193        );
194    }
195
196    #[test]
197    fn isobands_example() {
198        let matrix = vec![
199            vec![18., 13., 10., 9., 10., 13., 18.],
200            vec![13., 8., 5., 4., 5., 8., 13.],
201            vec![10., 5., 2., 1., 2., 5., 10.],
202            vec![9., 4., 1., 12., 1., 4., 9.],
203            vec![10., 5., 2., 1., 2., 5., 10.],
204            vec![13., 8., 5., 4., 5., 8., 13.],
205            vec![18., 13., 10., 9., 10., 13., 18.],
206            vec![18., 13., 10., 9., 10., 13., 18.],
207        ];
208        let (matrix, width, height) = make_grid_from2d_vec(&matrix);
209        let lower_band = 4.5;
210        let bandwidth = 4.5;
211
212        let res = isobands(
213            &matrix,
214            &vec![lower_band, lower_band + bandwidth],
215            false,
216            width,
217            height,
218            false,
219        )
220        .unwrap();
221        assert_eq!(
222            res[0].0,
223            vec![
224                vec![
225                    Point::new(1., 0.8),
226                    Point::new(0.8, 1.),
227                    Point::new(0.2, 2.),
228                    Point::new(0., 3.),
229                    Point::new(0., 3.),
230                    Point::new(0., 3.),
231                    Point::new(0.2, 4.),
232                    Point::new(0.8, 5.),
233                    Point::new(1., 5.2),
234                    Point::new(2., 5.8),
235                    Point::new(3., 6.),
236                    Point::new(3., 7.),
237                    Point::new(3., 7.),
238                    Point::new(3., 7.),
239                    Point::new(3., 6.),
240                    Point::new(4., 5.8),
241                    Point::new(5., 5.2),
242                    Point::new(5.2, 5.),
243                    Point::new(5.8, 4.),
244                    Point::new(6., 3.),
245                    Point::new(6., 3.),
246                    Point::new(6., 3.),
247                    Point::new(5.8, 2.),
248                    Point::new(5.2, 1.),
249                    Point::new(5., 0.8),
250                    Point::new(4., 0.2),
251                    Point::new(3., 0.),
252                    Point::new(3., 0.),
253                    Point::new(3., 0.),
254                    Point::new(2., 0.2),
255                    Point::new(1., 0.8),
256                ],
257                vec![
258                    Point::new(0.9, 3.),
259                    Point::new(1., 2.5),
260                    Point::new(1.1666666666666667, 2.),
261                    Point::new(2., 1.1666666666666667),
262                    Point::new(2.5, 1.),
263                    Point::new(3., 0.9),
264                    Point::new(3.5, 1.),
265                    Point::new(4., 1.1666666666666667),
266                    Point::new(4.833333333333333, 2.),
267                    Point::new(5., 2.5),
268                    Point::new(5.1, 3.),
269                    Point::new(5., 3.5),
270                    Point::new(4.833333333333333, 4.),
271                    Point::new(4., 4.833333333333333),
272                    Point::new(3.5, 5.),
273                    Point::new(3., 5.1),
274                    Point::new(2.5, 5.),
275                    Point::new(2., 4.833333333333333),
276                    Point::new(1.1666666666666667, 4.),
277                    Point::new(1., 3.5),
278                    Point::new(0.9, 3.),
279                ],
280                vec![
281                    Point::new(2.7272727272727275, 3.),
282                    Point::new(3., 2.7272727272727275),
283                    Point::new(3.2727272727272725, 3.),
284                    Point::new(3., 3.2727272727272725),
285                    Point::new(2.7272727272727275, 3.),
286                ],
287                vec![
288                    Point::new(3., 2.3181818181818183),
289                    Point::new(2.3181818181818183, 3.),
290                    Point::new(3., 3.6818181818181817),
291                    Point::new(3.6818181818181817, 3.),
292                    Point::new(3., 2.3181818181818183),
293                ],
294            ]
295        );
296    }
297
298    #[test]
299    fn isobands_original_code_issue_6_3_5() {
300        let matrix = vec![
301            vec![1., 1., 1., 1., 1., 1., 1.],
302            vec![1., 5., 5., 5., 5., 5., 1.],
303            vec![1., 5., 15., 15., 15., 5., 1.],
304            vec![1., 5., 10., 10., 10., 5., 1.],
305            vec![1., 5., 5., 5., 5., 5., 1.],
306            vec![1., 1., 1., 1., 1., 1., 1.],
307        ];
308        let (matrix, width, height) = make_grid_from2d_vec(&matrix);
309
310        let lower_band = 3.;
311        let bandwidth = 2.;
312
313        let res = isobands(
314            &matrix,
315            &vec![lower_band, lower_band + bandwidth],
316            false,
317            width,
318            height,
319            false,
320        )
321        .unwrap();
322        assert_eq!(
323            res[0].0,
324            vec![
325                vec![
326                    Point::new(1.0, 0.5),
327                    Point::new(0.5, 1.0),
328                    Point::new(0.5, 2.0),
329                    Point::new(0.5, 3.0),
330                    Point::new(0.5, 4.0),
331                    Point::new(1.0, 4.5),
332                    Point::new(2.0, 4.5),
333                    Point::new(3.0, 4.5),
334                    Point::new(4.0, 4.5),
335                    Point::new(5.0, 4.5),
336                    Point::new(5.5, 4.0),
337                    Point::new(5.5, 3.0),
338                    Point::new(5.5, 2.0),
339                    Point::new(5.5, 1.0),
340                    Point::new(5.0, 0.5),
341                    Point::new(4.0, 0.5),
342                    Point::new(3.0, 0.5),
343                    Point::new(2.0, 0.5),
344                    Point::new(1.0, 0.5)
345                ],
346                vec![
347                    Point::new(1.0, 2.0),
348                    Point::new(2.0, 1.0),
349                    Point::new(3.0, 1.0),
350                    Point::new(4.0, 1.0),
351                    Point::new(5.0, 2.0),
352                    Point::new(5.0, 3.0),
353                    Point::new(4.0, 4.0),
354                    Point::new(3.0, 4.0),
355                    Point::new(2.0, 4.0),
356                    Point::new(1.0, 3.0),
357                    Point::new(1.0, 2.0)
358                ]
359            ]
360        );
361    }
362
363    #[test]
364    fn isobands_original_code_issue_6_5_7() {
365        let matrix = vec![
366            vec![1., 1., 1., 1., 1., 1., 1.],
367            vec![1., 5., 5., 5., 5., 5., 1.],
368            vec![1., 5., 15., 15., 15., 5., 1.],
369            vec![1., 5., 10., 10., 10., 5., 1.],
370            vec![1., 5., 5., 5., 5., 5., 1.],
371            vec![1., 1., 1., 1., 1., 1., 1.],
372        ];
373        let (matrix, width, height) = make_grid_from2d_vec(&matrix);
374
375        let lower_band = 5.;
376        let bandwidth = 2.;
377
378        let res = isobands(
379            &matrix,
380            &vec![lower_band, lower_band + bandwidth],
381            false,
382            width,
383            height,
384            false,
385        )
386        .unwrap();
387        assert_eq!(
388            res[0].0,
389            vec![
390                vec![
391                    Point::new(1.0, 1.0),
392                    Point::new(1.0, 1.0),
393                    Point::new(1.0, 2.0),
394                    Point::new(1.0, 3.0),
395                    Point::new(1.0, 4.0),
396                    Point::new(1.0, 4.0),
397                    Point::new(2.0, 4.0),
398                    Point::new(3.0, 4.0),
399                    Point::new(4.0, 4.0),
400                    Point::new(5.0, 4.0),
401                    Point::new(5.0, 4.0),
402                    Point::new(5.0, 3.0),
403                    Point::new(5.0, 2.0),
404                    Point::new(5.0, 1.0),
405                    Point::new(5.0, 1.0),
406                    Point::new(4.0, 1.0),
407                    Point::new(3.0, 1.0),
408                    Point::new(2.0, 1.0),
409                    Point::new(1.0, 1.0)
410                ],
411                vec![
412                    Point::new(1.2, 2.0),
413                    Point::new(2.0, 1.2),
414                    Point::new(3.0, 1.2),
415                    Point::new(4.0, 1.2),
416                    Point::new(4.8, 2.0),
417                    Point::new(4.6, 3.0),
418                    Point::new(4.0, 3.6),
419                    Point::new(3.0, 3.6),
420                    Point::new(2.0, 3.6),
421                    Point::new(1.4, 3.0),
422                    Point::new(1.2, 2.0)
423                ]
424            ]
425        );
426    }
427
428    #[test]
429    fn isobands_multiple_bands() {
430        let matrix = vec![
431            vec![1., 1., 1., 1., 1., 1., 1.],
432            vec![1., 5., 5., 5., 5., 5., 1.],
433            vec![1., 5., 15., 15., 15., 5., 1.],
434            vec![1., 5., 10., 10., 10., 5., 1.],
435            vec![1., 5., 5., 5., 5., 5., 1.],
436            vec![1., 1., 1., 1., 1., 1., 1.],
437        ];
438        let (matrix, width, height) = make_grid_from2d_vec(&matrix);
439
440        let intervals = vec![3., 5., 7.];
441
442        let res = isobands(&matrix, &intervals, false, width, height, false).unwrap();
443
444        assert_eq!(res.len(), 2);
445
446        assert_eq!(
447            res[0].0,
448            vec![
449                vec![
450                    Point::new(0.9999999999999749, 1.0),
451                    Point::new(1.0, 0.9999999999999749),
452                    Point::new(2.0, 0.9999999999999749),
453                    Point::new(3.0, 0.9999999999999749),
454                    Point::new(4.0, 0.9999999999999749),
455                    Point::new(5.0, 0.9999999999999749),
456                    Point::new(5.000000000000025, 1.0),
457                    Point::new(5.000000000000025, 2.0),
458                    Point::new(5.000000000000025, 3.0),
459                    Point::new(5.000000000000025, 4.0),
460                    Point::new(5.0, 4.000000000000025),
461                    Point::new(4.0, 4.000000000000025),
462                    Point::new(3.0, 4.000000000000025),
463                    Point::new(2.0, 4.000000000000025),
464                    Point::new(1.0, 4.000000000000025),
465                    Point::new(0.9999999999999749, 4.0),
466                    Point::new(0.9999999999999749, 3.0),
467                    Point::new(0.9999999999999749, 2.0),
468                    Point::new(0.9999999999999749, 1.0),
469                ],
470                vec![
471                    Point::new(1.0, 0.5),
472                    Point::new(0.5, 1.0),
473                    Point::new(0.5, 2.0),
474                    Point::new(0.5, 3.0),
475                    Point::new(0.5, 4.0),
476                    Point::new(1.0, 4.5),
477                    Point::new(2.0, 4.5),
478                    Point::new(3.0, 4.5),
479                    Point::new(4.0, 4.5),
480                    Point::new(5.0, 4.5),
481                    Point::new(5.5, 4.0),
482                    Point::new(5.5, 3.0),
483                    Point::new(5.5, 2.0),
484                    Point::new(5.5, 1.0),
485                    Point::new(5.0, 0.5),
486                    Point::new(4.0, 0.5),
487                    Point::new(3.0, 0.5),
488                    Point::new(2.0, 0.5),
489                    Point::new(1.0, 0.5),
490                ]
491            ]
492        );
493        assert_eq!(res[0].1, 3.0);
494        assert_eq!(res[0].2, 5.0);
495
496        assert_eq!(
497            res[1].0,
498            vec![
499                vec![
500                    Point::new(1.0, 1.0),
501                    Point::new(1.0, 1.0),
502                    Point::new(1.0, 2.0),
503                    Point::new(1.0, 3.0),
504                    Point::new(1.0, 4.0),
505                    Point::new(1.0, 4.0),
506                    Point::new(2.0, 4.0),
507                    Point::new(3.0, 4.0),
508                    Point::new(4.0, 4.0),
509                    Point::new(5.0, 4.0),
510                    Point::new(5.0, 4.0),
511                    Point::new(5.0, 3.0),
512                    Point::new(5.0, 2.0),
513                    Point::new(5.0, 1.0),
514                    Point::new(5.0, 1.0),
515                    Point::new(4.0, 1.0),
516                    Point::new(3.0, 1.0),
517                    Point::new(2.0, 1.0),
518                    Point::new(1.0, 1.0)
519                ],
520                vec![
521                    Point::new(1.2, 2.0),
522                    Point::new(2.0, 1.2),
523                    Point::new(3.0, 1.2),
524                    Point::new(4.0, 1.2),
525                    Point::new(4.8, 2.0),
526                    Point::new(4.6, 3.0),
527                    Point::new(4.0, 3.6),
528                    Point::new(3.0, 3.6),
529                    Point::new(2.0, 3.6),
530                    Point::new(1.4, 3.0),
531                    Point::new(1.2, 2.0)
532                ]
533            ]
534        );
535        assert_eq!(res[1].1, 5.0);
536        assert_eq!(res[1].2, 7.0);
537    }
538
539    #[test]
540    /// Test that isobands returns the same result when using a quadtree or not (simple dataset)
541    fn isobands_simple_same_with_quadtree() {
542        let matrix = vec![
543            vec![1., 1., 1., 0.],
544            vec![1., 5., 5., 1.],
545            vec![0., 1., 1., 1.],
546        ];
547        let (matrix, width, height) = make_grid_from2d_vec(&matrix);
548        let lower_band = 1.;
549        let bandwidth = 1.;
550
551        let res1 = isobands(
552            &matrix,
553            &[lower_band, lower_band + bandwidth],
554            false,
555            width,
556            height,
557            false,
558        )
559        .unwrap();
560        let res2 = isobands(
561            &matrix,
562            &[lower_band, lower_band + bandwidth],
563            true,
564            width,
565            height,
566            false,
567        )
568        .unwrap();
569
570        assert_eq!(res1, res2);
571    }
572
573    #[test]
574    /// Test that isobands returns the same result when using a quadtree or not (volcano dataset)
575    fn isobands_volcano_same_with_quadtree() {
576        let volcano_str = include_str!("../tests/fixtures/volcano.json");
577        let raw_data: serde_json::Value = serde_json::from_str(volcano_str).unwrap();
578        let matrix: Vec<f64> = raw_data["data"]
579            .as_array()
580            .unwrap()
581            .iter()
582            .map(|x| x.as_f64().unwrap())
583            .collect();
584        let h = raw_data["height"].as_u64().unwrap() as usize;
585        let w = raw_data["width"].as_u64().unwrap() as usize;
586
587        let intervals = [
588            90., 95., 100., 105., 110., 115., 120., 125., 130., 135., 140., 145., 150., 155., 160.,
589            165., 170., 175., 180., 185., 190., 195., 200.,
590        ];
591
592        let res1 = isobands(&matrix, &intervals, false, w, h, false).unwrap();
593        let res2 = isobands(&matrix, &intervals, true, w, h, false).unwrap();
594
595        assert_eq!(res1, res2);
596    }
597
598    #[test]
599    /// Test that isobands returns the same result when using a quadtree or not
600    /// (dataset from https://observablehq.com/@mthh/stewarts-potentials-on-the-gpu)
601    fn isobands_pot_pop_same_with_quadtree() {
602        let pot_pop = include_str!("../tests/fixtures/pot_pop_fr.json");
603        let raw_data: serde_json::Value = serde_json::from_str(pot_pop).unwrap();
604        let matrix: Vec<f64> = raw_data["data"]
605            .as_array()
606            .unwrap()
607            .iter()
608            .map(|x| x.as_f64().unwrap())
609            .collect();
610        let h = raw_data["height"].as_u64().unwrap() as usize;
611        let w = raw_data["width"].as_u64().unwrap() as usize;
612
613        let intervals = [
614            0.001, 105483.25, 527416.25, 1054832.5, 2109665., 3164497.5, 4219330., 5274162.5,
615            6328995., 7383827.5, 8438660., 9704459., 10548326.,
616        ];
617
618        let res1 = isobands(&matrix, &intervals, false, w, h, false).unwrap();
619        let res2 = isobands(&matrix, &intervals, true, w, h, false).unwrap();
620
621        assert_eq!(res1, res2);
622    }
623}
contour_isobands/lib.rs

contour_isobands/
lib.rs
