viewof projName = Inputs.radio(
  new Map([
    ["Orthographic", "ortho"],
    ["Equirectangular", "equirectangular"],
    ["Hill", "hill"],
    ["Equal-Earth", "equalearth"],
    ["AzEqualArea", "azequalarea"]
  ]),
  {
    label: "Choose projection",
    value: "hill"
  }
)

d3 = require("d3@7", "d3-geo-projection@4", "d3-inertia@0.4")

netcdf = import("https://cdn.skypack.dev/netcdfjs@2.0.2?min").then((d) => d.NetCDFReader)

reader = new netcdf(await fetch("https://raw.githubusercontent.com/cfranken/s2_data/main/files/" + source.src).then(function(response) {
    return response.arrayBuffer();
  }))

sources = [
  {src: "S2_200_2000_median.nc", variable: "cloud_free_fraction02", label: "Median Ratio"},
]

source = sources[dataSource]

isize = reader.header.dimensions[0].size
jsize = reader.header.dimensions[1].size
zsize = reader.header.dimensions[2].size

isize

jsize

zsize

topojson     = require('topojson')

values = {
  var values = new Float32Array(reader.getDataVariable(source.variable));
  return values;
}

polygons = contours(cylinder(topdown(values))).map(invert).map(function(p) {
  return p.type == "Polygon" ? simplify(p, 0.01) : p;
})

function simplify(polygon, s) {
    var topology = topojson.topology({polygon: polygon});
    topology = topojson.simplify(topojson.presimplify(topology), s);
    topology = topojson.filter(topology, r => topology.arcs[r[0]].length > 2);
    var p = topojson.feature(topology, topology.objects.polygon);
    p.value = polygon.value;
    return p;
  }

function invert(d) {
    var shared = {};
    var p = {
      type: "Polygon",
      coordinates: d3.merge(d.coordinates.map(function(polygon) {
        return polygon.map(function(ring) {
          return ring.map(function(point) {
            return [point[0] / (isize + 1) * 360 - 180, 90 - point[1] / jsize * 180];
          }).reverse();
        });
      }))
    };

    // Add contact points with the antimeridian.
    var epsilon = 1e-5;
    p.coordinates = p.coordinates
    .map(function(ring) {
      ring.pop();
      var r = [];
      for (var i = 0, l = ring.length; i < l; i++) {
        var p = ring[i];
        if (p[0] === 180 && ring[(i+l-1)%l][0] < 180)
            r.push([p[0], p[1] + epsilon]);
        if (p[0] === -180 && ring[(i+l-1)%l][0] > -180)
            r.push([p[0], p[1] - epsilon]);
        r.push(p);
        if (p[0] === 180 && ring[(i+1)%l][0] < 180)
            r.push([p[0], p[1] - epsilon]);
        if (p[0] === -180 && ring[(i+1)%l][0] > -180)
            r.push([p[0], p[1] + epsilon]);
      }
      r.push(r[0]);
      return r;
    });

    // Record the y-intersections with the antimeridian.
    p.coordinates.forEach(function(ring) {
      ring.forEach(function(p) {
        if (p[0] === -180 || p[0] === 180) {
          shared[p[1]] |= p[0] === -180 ? 1 : 2;
        }
      });
    });

    // Offset any unshared antimeridian points to prevent their stitching.
    p.coordinates.forEach(function(ring) {
      ring.forEach(function(p) {
        if ((p[0] === -180 || p[0] === 180) && shared[p[1]] !== 3) {
          p[0] = p[0] === -180 ? -179.9995 : 179.9995;
        }
      });
    });

    p = d3.geoStitch(p);

    // After stitching remove the shared points -- they are off
    // see https://github.com/d3/d3-contour/issues/25
    p.coordinates = p.coordinates
    .map(function(ring) {
      ring.pop();
      ring = ring
      .map(function(p) {
        var p0 = (p[0] + 180) * (isize + 1) / isize - 180;
        return [p0, p[1]];
      });
      ring = ring.filter(function(p) {
        return Math.abs(p[0]) < 179.9999;
      });
      ring.push(ring[0]);
      return ring;
    })
    .filter(function(ring) {
      return ring.length > 0;
    });

    // If the MultiPolygon is empty, treat it as the Sphere.
    return p.coordinates.length
        ? {type: "Polygon", coordinates: p.coordinates, value: d.value}
        : {type: "Sphere", value: d.value};
    
  }

function topdown(values) {
  var v = new Float32Array(values);
  for (var j = 0; j < jsize; j++) {
    for (var i = 0; i < isize; i++) {
      v[ (jsize-1-j) * isize + i ] = values[ j * isize + i ];
    }
  }
  return v;
}

// torus: copy left column to the right
function cylinder(values) {
  var v = new Float32Array(jsize * (isize+1));
  for (var j = 0; j < jsize; j++) {
    for (var i = 0; i < isize + 1; i++) {
      v[ j * (isize+1) + i ] = values[ j * isize + i%isize ];
     }
 }
 return v;
}

contours = d3.contours()
      .smooth(true)
      .size([isize + 1, jsize])

polygons