This refactors the STAC metadata -> xarray coordinates logic to support multi-dimensional coordinates. For example, imagine a field likerescaling (orhref, for that matter), that's different for every asset of every item. Previously, we'd drop this field. Now, we can store it as a 2D array indexing the dimensionstime, band.

This builds off of@Berhinj's work in#222, which provided a lot of inspiration for this design.

This is a significant rewrite of the metadata logic. The basic idea now is, for each field in STAC metadata:

1. Make a 2D NumPy array (for thetime, band dimensions)
2. Iterate through all the STAC metadata, and write values into the array
3. De-duplicate the array: if all values along a dimension are identical, drop that dimension. For example, thetype field will be an MxN array ofimage/tiff; application=geotiff over and over; this is collapsed down into a 0D scalar.

Unfortunately, it does come at a slight performance cost:stackstac.stack on 10,000Landast-8 items is 13.6% slower, going from 13.2s to 15s. (Of course, there's no performance difference with actually computing the results, which is where the bulk of time is spent anyway.)

There are a number of benefits, however:

• Robust handling of JSON-y subfields. For example, theroles field, which contains variable-length lists of tags per band like["data"],["cloud", "cloud-shadow"], is now stored. Same with complex things likeclassification:bitfields. More interestingly, we could now easily retaingeometry,bbox, etc. xrefStore per-item bounding boxes as a coordinate #6
• Generally, should be more robust to new fields/extensions adopted in STAC, or even inserted by users (they may not be formatted ideally, but almost certainly won't be dropped)
• raster:bands subfields are now in coordinates, and supporting extensions likeeo:bands,raster:bands, etc. is now generalized and easy to extend in the future
• Finally have tests for coordinates logic!
• rescale logic could be simplified into a one-liner after the stack is built (stack * stack.coords['raster:bands_scale'] + stack.coords['raster:bands_offset']), rather than plumbing scale/offset factors through many layers. (TBD if this makes a nasty Dask graph though.)

There's maybe even a future where the logic to generate the Dask array takes these coordinates as input, rather than raw STAC dicts or pystac items? Perhaps if we can someday store the geometries as a GeoSeries, or something like that.

Closes#216