mthh/contour-rs

Contour polygon creation in Rust (using marching squares algorithm)

contour-rs

Computes isolines, contour polygons and isobands by applying marching squares to a rectangular array of numeric values.
Outputs ring coordinates or geo-types geometries.

The results can also easily be serialised to GeoJSON.

Note : The core of the algorithm is ported from d3-contour.

Usage

Add this to your Cargo.toml:

[dependencies]
contour = "0.13.1"

and this to your crate root:

extern crate contour;

The API exposes:

• a ContourBuilder struct, which computes isorings coordinates for a Vec of threshold values and transform them either :

  • in Contours (a type containing the threshold value and the geometry as a MultiPolygon), or,
  • in Lines (a type containing the threshold value and the geometry as a MultiLineString).
  • in Bands (a type containing a minimum value, a maximum value and the geometry as a MultiPolygon).

• a contour_rings function, which computes isorings coordinates for a single threshold value (returns a Vec of rings coordinates - this is what is used internally by the ContourBuilder).

ContourBuilder is the recommended way to use this crate, as it is more flexible and easier to use (it enables to specify the origin and the step of the grid, and to smooth the contours, while contour_rings only speak in grid coordinates and doesn't smooth the resulting rings).

Line, Contour and Band can be serialised to GeoJSON using the geojson feature.

Example:

Without defining origin and step:

let c = ContourBuilder::new(10, 10, false); // x dim., y dim., smoothing
let res = c.contours(&vec![
    0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
    0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
    0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
    0., 0., 0., 1., 1., 1., 0., 0., 0., 0.,
    0., 0., 0., 1., 1., 1., 0., 0., 0., 0.,
    0., 0., 0., 1., 1., 1., 0., 0., 0., 0.,
    0., 0., 0., 1., 1., 1., 0., 0., 0., 0.,
    0., 0., 0., 1., 1., 1., 0., 0., 0., 0.,
    0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
    0., 0., 0., 0., 0., 0., 0., 0., 0., 0.
], &[0.5])?; // values, thresholds

With origin and step

let c = ContourBuilder::new(10, 10, true) // x dim., y dim., smoothing
    .x_step(2) // The horizontal coordinate for the origin of the grid.
    .y_step(2) // The vertical coordinate for the origin of the grid.
    .x_origin(100) // The horizontal step for the grid
    .y_origin(200); // The vertical step for the grid

let res = c.contours(&[
0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 1., 1., 0., 1., 0., 0., 0.,
0., 0., 0., 1., 1., 0., 1., 0., 0., 0.,
0., 0., 0., 1., 1., 0., 1., 0., 0., 0.,
0., 0., 0., 1., 1., 0., 1., 0., 0., 0.,
0., 0., 0., 1., 1., 0., 1., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0.
], &[0.5]).unwrap(); // values, thresholds

Using the geojson feature

The geojson feature is not enabled by default, so you need to specify it in your Cargo.toml:

[dependencies]
contour = { version = "0.13.1", features = ["geojson"] }

let c = ContourBuilder::new(10, 10, true) // x dim., y dim., smoothing
    .x_step(2) // The horizontal coordinate for the origin of the grid.
    .y_step(2) // The vertical coordinate for the origin of the grid.
    .x_origin(100) // The horizontal step for the grid
    .y_origin(200); // The vertical step for the grid

let res = c.contours(&[
0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 1., 1., 0., 1., 0., 0., 0.,
0., 0., 0., 1., 1., 0., 1., 0., 0., 0.,
0., 0., 0., 1., 1., 0., 1., 0., 0., 0.,
0., 0., 0., 1., 1., 0., 1., 0., 0., 0.,
0., 0., 0., 1., 1., 0., 1., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0.
], &[0.5]).unwrap(); // values, thresholds
println!("{:?}", res[0].to_geojson()); // prints the GeoJSON representation of the first contour

Output:

Feature {
    bbox: None,
    geometry: Some(Geometry {
        bbox: None,
        value: MultiPolygon([
            [[
                [110.0, 215.0], [110.0, 213.0], [110.0, 211.0], [110.0, 209.0],
                [110.0, 207.0], [109.0, 206.0], [107.0, 206.0], [106.0, 207.0],
                [106.0, 209.0], [106.0, 211.0], [106.0, 213.0], [106.0, 215.0],
                [107.0, 216.0], [109.0, 216.0], [110.0, 215.0]
            ]],
            [[
                [114.0, 215.0], [114.0, 213.0], [114.0, 211.0], [114.0, 209.0],
                [114.0, 207.0], [113.0, 206.0], [112.0, 207.0], [112.0, 209.0],
                [112.0, 211.0], [112.0, 213.0], [112.0, 215.0], [113.0, 216.0],
                [114.0, 215.0]
            ]]
        ]),
        foreign_members: None
    }),
    id: None,
    properties: Some({"threshold": Number(0.5)}),
    foreign_members: None
}

Using the f32 feature

By default, this crate expects f64 values as input and uses f64 values for its computations. If you want to use f32 values instead, you need to specify the f32 feature in your Cargo.toml:

[dependencies]
contour = { version = "0.13.1", features = ["f32"] }

WASM demo

Demo of this crate compiled to WebAssembly and used from JavaScript : wasm_demo_contour.

Difference with the contour-isobands crate (from mthh/contour-isobands-rs repository)

While this crate computes isolines (cf. wikipedia:Marching_squares) from which it derives contour polygons (i.e. polygons that contain all points above the threshold defined for a given isoline) and isobands (i.e. polygons that contain all points between a minimum and a maximum bound), contour-isobands-rs is only dedicated to compute isobands and use a slightly different implementation of the marching squares algorithm for the disambiguation of saddle points (cf. wikipedia:Marching_squares).

License

Licensed under either of

• Apache License, Version 2.0, (LICENSE-APACHE or http://www.apache.org/licenses/LICENSE-2.0)
• MIT license (LICENSE-MIT or http://opensource.org/licenses/MIT)

at your option.

Contribution

Unless you explicitly state otherwise, any contribution intentionally submitted for inclusion in the work by you, as defined in the Apache-2.0 license, shall be dual licensed as above, without any additional terms or conditions.

contour-rs's Issues

Area calculation via shoelace method ignores final vertex (which probably gives wrong result for small polygons)

I think the area function in src/area.rs is not handling the wraparound case correctly:

pub fn area(ring: &[Pt]) -> Float {
    let mut i = 0;
    let n = ring.len() - 1;
    let mut area = ring[n - 1].y * ring[0].x - ring[n - 1].x * ring[0].y;
    while i < n {
        i += 1;
        area += ring[i - 1].y * ring[i].x - ring[i - 1].x * ring[i].y;
    }
    area
}

Say ring has 10 elements; then n == 9 and so the initialization of area is pairing up ring[0] and ring[8] (rather than ring[9]).

A correct implementation could be something like:

pub fn area(poly: &[Pt]) -> Float {
    let n = poly.len();
    let mut area = poly[n - 1].y * poly[0].x - poly[n - 1].x * poly[0].y;
    for i in 1..n {
        area += poly[i - 1].y * poly[i].x - poly[i - 1].x * poly[i].y;
    }
    area
}

(For anyone else who stumbles on this issue: I think normally you'd also want to divide the result by two, but this library is only using area to detect polygon winding and relative size of areas, so it shouldn't be necessary here.)

[FR] Sparse representation of input

I have a dataset where the vast majority of values are nodata values (#16).
The bounding box is substantially larger than what would be required if the data were e.g. tiled or stored in some other compressed format. This prevents me from doing calculations at the precision I'd like to, as I'm running out of memory.
I think it would be useful if contour-rs accepted grids of formats other than Vec. Maybe making the methods generic over some trait that provides fn get_value(x: usize, y: usize) -> Option<V> could work? Such a trait could also help solving #16, since an implementation could return None for such values.

Apparently there are also optimized variants of marching squares for sparse matrices, maybe it would be useful to choose a type of input that would allow to take advantage of them, e.g. https://www.diva-portal.org/smash/record.jsf?pid=diva2%3A1675180&dswid=-8922

Avoid geojson types as internal representation

Every time there is a new release of geojson, the types of this crate must be updated.

This isn't a huge problem - you can just keep bumping the geojson crate (82b2373, 6398263, 520f39d) as you have, but really I don't know that there's a good reason to use the types from the geojson crate as your interchange format.

The downsides are, as mentioned above, that it triggers a downstream update cycle, but also those data structures aren't very efficient ways of representing your data types.

Would you consider a PR which leverages internal types detached from the geojson crate for contours?

If you'd like, it could retain an impl Serialize that serializes to geojson or (my preference) add a to_geojson method that explicitly converts your internal types to geojson types for those people who actually want to serialize geojson.

[FR] Parallelism support

The sister crate over at https://github.com/mthh/contour-isobands-rs supports parallelism via the parallel feature.
Is there a technical reason this crate does not support it as well, or did just no one need it until now 🤔
(Also, in case we do end up adding it here, maybe it would make sense to merge the two crates? Seems they have some shared code, and either algorithm could be exposed via different structs. I'd imagine stuff like #14 would be useful there as well, especially considering the heightened memory usage.)

[BUG] Invalid polygons/closed linestrings when smoothing is enabled

When smoothing is enabled, ContourBuilder::contours and probably ContourBuilder::isobands can end up generating invalid (self-intersecting) polygons.
E.g.:
Without smoothing:

With smoothing:

Overlaid:

I think that this is invalid according to geo-types, and many other geometry standards don't allow this either.
Specifically,

The interior of the polygon must be a connected point-set. That is, any two distinct points in the interior must admit a curve between these two that lies in the interior.

And,

The exterior and interior rings must be valid LinearRings (see LineString).

Further, a closed LineString must not self-intersect. Note that its validity is not enforced, and operations and predicates are undefined on invalid LineStrings.

@mthh You can find the areas the screenshots are from in the dataset I sent you recently at about 557500, 5932680 in its CRS, contours were generated at elevation -5.

Calculate separate isolines

Hi,
First of all, thanks you for providing this implementation.

I'm want to calculate two isolines based on following grid, with a treshold of 1:

    // [
    //   [0.3, 1.5, 0.6],
    //   [0.6, 1.5, 0.6],
    //   [0.3, 1.5, 0.6],
    // ];

We should have two vertical lines, but I end with one circular line. Seems that the code is always trying to obtain a ring. Is there a way to obtain separates lines ?

Here is a working test, that fails:

pub fn isolines(
    x0: f64,
    y0: f64,
    x_step: f64,
    y_step: f64,
    dx: u32,
    dy: u32,
    grid: &[f64],
    thresholds: &[f64],
) -> Result<Vec<Line>, Box<dyn Error>> {
    let res = ContourBuilder::new(dx, dy, true)
        // set origin 
         .x_origin(x0 - x_step / 2.)
        .y_origin(y0 - y_step / 2.)
        .x_step(x_step)
        .y_step(y_step)
        .lines(grid, thresholds)?;

    Ok(res)
}

#[cfg(test)]
mod isolines_test {
    use contour::Line;
    use geo_types::{line_string, MultiLineString};

    use super::isolines;
   
     // transform isolines output into something we can compare with geo_types
    fn as_multi_line(lines: Vec<Line>) -> Vec<MultiLineString> {
        let multi_lines_output: Vec<MultiLineString> =
            lines.iter().map(|l| l.geometry().to_owned()).collect();
        multi_lines_output
    }

 #[test]
    fn create_two_lines_3x3_with_two_vertical_isolines() {
        // Grid data
        let grid: Vec<f64> = vec![0.3, 1.5, 0.6, 0.6, 1.5, 0.6, 0.3, 1.5, 0.6];

        // Thresholds
        let thresholds: Vec<f64> = vec![1.0];

        // Calculate isolines
        let lines = isolines(10., 10., 1., 1., 3, 3, &grid, &thresholds).unwrap();

        let multi_lines_output = as_multi_line(lines);

        // We want two lines
        let expected_lines: Vec<MultiLineString> = vec![
            line_string![
                (x: 10.583333333333334, y: 10.0),
                (x: 10.444444444444445, y: 11.0),
                (x: 10.583333333333334, y: 12.0),
            ]
            .into(),
            line_string![
                (x: 11.555555555555555, y: 12.0),
                (x: 11.555555555555555, y: 11.0),
                (x: 11.555555555555555, y: 10.0),
            ]
            .into(),
        ];

        // Fails: we have only a ring
        assert_eq!(multi_lines_output, expected_lines);
    }
}

Expose an error type

Currently, methods return an error type, which is impossible to access from other crates, and so it is impossible to make any reasonable error handling beyond "something gone wrong".

Error struct is public, but error module is not, and it is not re-exported, so we can't name this type in a code and can't match on its kind,

Overflow/Artifacts with large rasters (>40k*40k)

If I pass in a raster larger than about 40000*40000 (not sure about the exact value, but 50k*50k fails) to ContourBuilder::contours, I get an error similar to this:

thread 'main' panicked at /home/nett/.cargo/registry/src/index.crates.io-6f17d22bba15001f/contour-0.12.1/src/isoringbuilder.rs:137:23:
index out of bounds: the len is 4294836225 but the index is 18446744071562067968
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

Increasing the width of the involved local variables to i64 seems to have fixed this issue, and I'll create a PR containing those changes once I test it against real data.

[FR] Nodata value and incomplete ring support

I want to be able to generate contours that do not include lines on data/nodata borders, 'nodata' being any value either outside the grid or a specifically chosen one inside of it, commonly f32::MAX or f64::MAX, e.g. (nodata marked red):
Without nodata handling:

With nodata handling:

This also prevents a ring forming around the whole dataset, as seen in #7. While I'm aware of the workaround proposed in the issue, it is not sufficiently performant for me.
I started implementing this feature here: https://github.com/SenseLabsDE/contour-rs/tree/nodata, but the current version is relatively inefficient and results in slowdowns even if the feature is unused.