Extract all bands from a Contour Plot [UNLISTED]

Contour plots generated by d3-contour overlay shapes (Multipolygons) of increasing value. To obtain the shape that is between two values a < b, we need to substract the shape B = contour(b) from the shape A = contour(a).

by definition, all polygons from B are included in 1 (or several) polygons of A, and there is no intersection (otherwise we would have value x simultaneously equal to a and to b)
each outer ring from B (where x > b) must be transformed into a hole in its encompassing polygon(s) of A
holes from B must be added as polygons to A.

The substraction algorithm for this case is then much simpler than a generic Multipolygon intersection/substraction:

For each polygon in multipolygon B, take a point from its outer ring. Then find the polygon of A that contains this point, and add that ring as a hole (which implies a reversal of its coordinates). Plus, add holes of B as polygons to A.

But I'm not sure this covers all cases, that is, what happens when a hole of A is inside a hole of B?

A simpler method works OK if one wants to extract 1 band (/fil/94d4df7520ffc8434c40fe9b82ebd536), but it is useless if you need to extract several contiguous bands, as it does not give an exact boundary (there are gaps and overlaps, cf /fil/3f866e2e90c3e019bfe3fd9e0d43fe14).

Forked from mbostock's block: Contour Plot

forked from Fil's block: Extract a band from a Contour Plot [UNLISTED]

<!DOCTYPE html>
<svg width="960" height="673" stroke="#fff" stroke-width="0.5"></svg>
<script src="https://d3js.org/d3.v4.min.js"></script>
<script src="https://d3js.org/d3-hsv.v0.1.min.js"></script>
<script src="https://d3js.org/d3-contour.v1.min.js"></script>
<script>

var svg = d3.select("svg"),
    width = +svg.attr("width"),
    height = +svg.attr("height");

var limit = 155;
  
var i0 = d3.interpolateHsvLong(d3.hsv(120, 1, 0.65), d3.hsv(60, 1, 0.90)),
    i1 = d3.interpolateHsvLong(d3.hsv(60, 1, 0.90), d3.hsv(0, 0, 0.95)),
    interpolateTerrain = function(t) { return t < 0.5 ? i0(t * 2) : i1((t - 0.5) * 2); },
    color = d3.scaleSequential(interpolateTerrain).domain([90, 190]);

d3.json("volcano.json", function(error, volcano) {
  if (error) throw error;

var contours = d3.contours()
        .size([volcano.width, volcano.height])
        .thresholds(d3.range(90, 195, 5))
      (volcano.values);
  svg.selectAll("path")
    .data(contours)
    .enter().append("path")
      .attr("d", d3.geoPath(d3.geoIdentity().scale(width / volcano.width)))
      .attr("stroke", function(d) { return color(d.value); })
  .attr('stroke-width', 3)
  .attr('fill', 'none');

svg.selectAll('circle')
  .data(volcano.values)
  .enter()
  .append('circle')
  .attr('r', 2)
  .attr('cx', (d,i) => {
    return (i % volcano.width) * width / volcano.width;
  })
  .attr('cy', (d,i) => {
    return (i - (i % volcano.width)) / volcano.width * (width / volcano.width);
  })
  .attr('fill', color)
  .append('title').text(d => d);
  
  contours = contours.slice(1)
  contours.forEach((K,i) => {
     var holes = [];
  K.coordinates = K.coordinates.map(d => {
    var k = contours[i+1];
    if (k){
      
    k.coordinates.forEach(e => {
      if (d3.polygonContains(d[0], e[0][0])) {
        // e's outer ring becomes a hole in d
        d.push(e[0].slice().reverse());
        // e's holes become new polygons
        e.slice(1).forEach(
          hole => {  
            holes.push([hole.slice().reverse()]);
          }
        );
      }
    });
    }
    return d;
  });
    holes.forEach(hole => K.coordinates.push(hole))
  ;

if (i%3)
    svg.append('path')
      .attr("fill", color(90 + 5 * i))
      .attr('opacity', 1)
        .datum(K)
       .attr("d", d3.geoPath(d3.geoIdentity().scale(width / volcano.width)))
  });

});

</script>

https://d3js.org/d3.v4.min.js

https://d3js.org/d3-hsv.v0.1.min.js

https://d3js.org/d3-contour.v1.min.js