This implementation is feeling pretty good to me now, but I’m beginning to see what you were going after with #1210: currently there’s not an easy way for the Contour mark to leverage the nearest/barycentric/random-walk interpolation methods to convert non-gridded samples into a grid suitable for d3-contour. But clearly we would like to be able to use them e.g. to produce contours over the ca55 dataset.

So I think what I’d like to do is change the rasterize function option to return a dense array of gridded values (in data space, without scaling), and probably rename it back to interpolate as we originally imagined. That’ll make things a little slower for the raster mark since it’s another copy of the data, and it means that the nearest interpolation method won’t be able to use the Canvas2D API to render, but I think the performance difference will probably be small enough that it won’t be noticeable, and it will allow the Contour mark to reuse these implementations.

Let me know if that sounds good or please take a crack at it as you like, @Fil!

The contours for ca55 are very pretty with more blurring (e.g. blur: 3), and the barycentric interpolate doesn't have the arbitrary (random) high frequency component of random-walk. So I'd suggest to change the test to barycentric+blur: 3, to make it a better example. (DONE)

The contours are empty if we don't specify an interpolate strategy; I think this would probably be addressed by d3/d3-contour#61: it would make contours around the samples. But I don't think that's what people will want in general (and certainly not on this dataset), and I'd suggest we either:

• throw an error or give a warning to say "the samples don't cover the whole raster, please select an interpolate strategy (e.g.: none, or barycentric)"
• default to barycentric+blur3

It would be good to backport this to d3-contour at some point, to allow operations on the returned array of GeoJSON?

I think blurring + random-walk is also a good solution to #1196 (comment) (see dd06ea05c commit log).

• Fixed two bugs I found while playing with the contours and rasters
• Checked degenerate cases such as n=1 or n=2 points.

About the shorthand and defaults, I suggest we follow a similar convention as in Plot.density:

• default to black stroke, no fill; use fill: value if we want a colored fill, same for stroke (for the density mark, the special keyword is "density"; here we can use "value" or "threshold").
• the value of a point is specified with the value option, as a channel if we have data, and as a function f(x,y) if we accept a function sampler.

For the record here's what the vapor-barycentric looked like :)

It's in a test I haven't pushed because it's very slow—I'll try to make a lighter one, but for now I've built 11dbb07 and added it to https://observablehq.com/@observablehq/global-temperatures-raster#issue11dbb0; you can see that it fails to create the contours.

Ah no! It works—it's just that the defaults have changed 👍

I think there are still some issues with the concept of a “dense” grid, and in particular whether this is a linear grid in data space or in screen coordinates. For example, say we start here:

Plot.plot({
  marks: [Plot.raster(await vapor(), {width: 360, height: 180, fill: (d) => d})]
})

We can fix the y inversion y reversing the y scale:

Plot.plot({
  y: {reverse: true},
  marks: [Plot.raster(await vapor(), {width: 360, height: 180, fill: (d) => d})]
})

But what if we wanted instead to specify the x domain as [-180, 180] and the y domain as [90, -90]? Then we get something surprising because the defaults for x1, y1, x2, y2 and 0, 0, width, height.

Plot.plot({
  x: {domain: [-180, 180]},
  y: {domain: [90, -90]},
  marks: [Plot.raster(await vapor(), {width: 360, height: 180, fill: (d) => d})]
})

We can fix that by specifying these options explicitly, I suppose.

Plot.plot({
  x: {domain: [-180, 180]},
  y: {domain: [90, -90]},
  marks: [Plot.raster(await vapor(), {x1: -180, y1: -90, y2: 90, x2: 180, width: 360, height: 180, fill: (d) => d})]
})

But notice what happens when we change the y scale type of sqrt: the raster stays the same, even though the axis has changed.

Plot.plot({
  y: {type: "sqrt"},
  marks: [Plot.raster(await vapor(), {width: 360, height: 180, fill: (d) => d})]
})

In other words, when you use the dense grid shorthand, you have to provide a dense grid in scaled (projected) coordinates, which makes it rather difficult to use. I get more what I expect if I specify x and y explicitly (non-optimized).

Plot.plot({
  y: {type: "sqrt"},
  marks: [
    Plot.raster(await vapor(), {
      width: 360,
      height: 180,
      x: (d, i) => (i % 360) - 180 + 0.5,
      y: (d, i) => 90 - ((i / 360) | 0) + 0.5,
      fill: (d) => d,
      interpolate: "nearest"
    })
  ]
})

Also note that if I try to use the interpolate option with a dense grid, it crashes because the x and y channels are missing.

TypeError: Cannot read properties of undefined (reading '362')

Plot.plot({
  y: {type: "sqrt"},
  marks: [
    Plot.raster(await vapor(), {
      width: 360,
      height: 180,
      fill: (d) => d,
      interpolate: "nearest"
    })
  ]
})

I think the challenge with “dense grid in abstract rather than projected coordinates” is that it doesn’t work as well for evaluating continuous functions. When you have a continuous function, we want to compute the grid in projected (scaled) coordinates and then invert it back into abstract (data) coordinates for the best results. Maybe we need to separate the code path for the continuous function case somehow; we could still provide shorthand for the dense grid case but not worry about optimizing it so much.

Okay, I think this is looking really good now! There’s still some optimization that’s likely possible when you specify a dense grid of values, but only in the case where x and y are linear scales. If we figure out how to make that faster in the future that’s gravy but it shouldn’t change the API! From my side I’m ready to merge this and the parent branch into main and start working on documentation.