# Astronomy reference stores — SDSS DR17 spectra

Public, anonymously readable **virtual-reference** datasets: kilobyte-sized
[Icechunk](https://icechunk.io) repositories of byte-range pointers into original
**FITS** files. No pixels were copied into a new layout, nothing was resampled, and
nothing was resharded into a sample format. Opening one gives you an ordinary
zarr-v3 array; reading it fetches ranged GETs from the FITS files themselves.

Scan once, use many: building these required a full pass over the FITS headers with
[VirtualiZarr](https://github.com/zarr-developers/VirtualiZarr) / `kerchunk.fits`.
That pass is done. You need only `icechunk` to read them.

Everything here is in the public domain — see
[`sdss/dr17/PROVENANCE.md`](sdss/dr17/PROVENANCE.md) for sources, the SDSS licence
statement, and citation guidance.

## What is here

| store | shape | chunks | samples/chunk | pixels fetched from |
| --- | --- | --- | --- | --- |
| `sdss_dr17_p0266_refs` | `(640, 3864)` f4 | `(64, 3864)` | 64 | `data.sdss.org` |
| `sdss_dr17_6plate_refs` | `(3840, 3841)` f4 | `(1, 3841)` | 1 | `data.sdss.org` |
| `sdss_dr17_6plate_mirror_refs` | `(3840, 3841)` f4 | `(1, 3841)` | 1 | this bucket (mirror) |
| `sdss/dr17/…` | — | — | — | the six mirrored `spPlate` FITS files (341 MiB) |

All three expose one variable, **`flux`**, as `(fiber, wavelength)` `float32`. Axis 0
is one **fiber** — a single galaxy/star/quasar spectrum. Axis 1 is a **common
log-wavelength grid**, so spectra are directly comparable row to row with no
resampling.

Fluxes are in SDSS units of 10⁻¹⁷ erg s⁻¹ cm⁻² Å⁻¹. Bad pixels are `NaN`.

**Which to use.** `sdss_dr17_6plate_mirror_refs` unless you have a reason otherwise:
same geometry as `6plate`, but the references point at the FITS mirror in this
bucket, so it is one provider and, from Google Cloud, in-region.
`data.sdss.org` rate-limits new connections and this layout needs thousands of
ranged reads.

### The two layouts are a real trade-off

`spPlate` frames hold ~640 fibers on one wavelength grid, contiguous fiber-major. So a
block of fibers is a contiguous byte range and can become a chunk by arithmetic alone.

- **`p0266_refs` — one plate, 64 fibers per chunk.** Full wavelength width. One read
  decodes 64 spectra, so per-chunk cost is amortized across many samples.
- **`6plate_refs` — six plates, one fiber per chunk.** Different plates cover slightly
  different wavelength ranges, so joining them means cropping each to the shared
  window — which breaks fiber contiguity, leaving one fiber per contiguous range.
  You get archive-scale concatenation (3840 spectra on one flat axis, extending to the
  ~2800-plate archive) and you pay one GET per spectrum.

The crop is exact, not approximate: every SDSS plate uses the same `dloglam` step and
the plates' `COEFF0` start wavelengths differ by whole bins, so the windows align onto
a single grid with no interpolation.

## Read one

Only `icechunk` is required (`pip install icechunk zarr`). No credentials, no
requester-pays, no build step.

```python
import icechunk
import zarr

BUCKET = "insitubatch-bench-insitubatch"
PREFIX = "astronomy/sdss_dr17_6plate_mirror_refs"
# where this store's virtual chunks live -- the mirror, for this store
PIXELS = f"https://storage.googleapis.com/{BUCKET}/astronomy/sdss/dr17/"

config = icechunk.RepositoryConfig.default()
config.set_virtual_chunk_container(
    icechunk.VirtualChunkContainer(PIXELS, icechunk.http_store())
)
repo = icechunk.Repository.open(
    icechunk.gcs_storage(bucket=BUCKET, prefix=PREFIX, anonymous=True),
    config=config,
    authorize_virtual_chunk_access=icechunk.containers_credentials({PIXELS: None}),
)

flux = zarr.open_array(repo.readonly_session("main").store, path="flux", mode="r")
print(flux.shape, flux.dtype, flux.chunks)   # (3840, 3841) float32 (1, 3841)

spectrum = flux[0]           # one fiber
block = flux[64:128]         # 64 fibers

# each store records the prefix its virtual chunks resolve against:
print(repo.config.virtual_chunk_containers)
```

For the two `data.sdss.org`-backed stores, set `PIXELS = "https://data.sdss.org/"` instead.

### Wavelengths

The grid is uniform in log₁₀(wavelength): bin `i` is `10**(COEFF0 + i*COEFF1)` Å, the
SDSS convention, taken from the `spPlate` primary headers. Each store carries its own
solution as array attributes — `flux.attrs["COEFF0"]` and `flux.attrs["COEFF1"]` — so you
do not need to hard-code the table below.

| store | `COEFF0` | `COEFF1` | bins | range |
| --- | --- | --- | --- | --- |
| `sdss_dr17_p0266_refs` | `3.5785` | `0.0001` | 3864 | 3789 – 9221 Å |
| `sdss_dr17_6plate_refs`, `…_mirror_refs` | `3.5797` | `0.0001` | 3841 | 3799 – 9198 Å |

The six-plate stores use `3.5797` because it is the *shared* window start — the
largest of the six plates' `COEFF0` values, which is where the common grid begins.

```python
import numpy as np

COEFF0, COEFF1, N_BINS = 3.5797, 0.0001, 3841       # or read them off flux.attrs
wavelength_angstrom = 10.0 ** (COEFF0 + COEFF1 * np.arange(N_BINS))
wavelength_angstrom[[0, -1]]                        # array([3799.3, 9198.1])
```

## Read it as training batches

One fiber per chunk means one ranged GET per spectrum, so throughput depends entirely
on read concurrency. Any loader works; below is
[insitubatch](https://github.com/emfdavid/insitubatch), which streams shuffled,
split-aware batches straight from the store with no local copy:

```bash
pip install "insitubatch[torch]" icechunk
```

Complete and standalone — it does not continue from the block above:

```python
import icechunk
from insitubatch import InSituDataset, open_geometries, split_by_chunk

BUCKET = "insitubatch-bench-insitubatch"
PREFIX = "astronomy/sdss_dr17_6plate_mirror_refs"
PIXELS = f"https://storage.googleapis.com/{BUCKET}/astronomy/sdss/dr17/"

config = icechunk.RepositoryConfig.default()
config.set_virtual_chunk_container(
    icechunk.VirtualChunkContainer(PIXELS, icechunk.http_store())
)
repo = icechunk.Repository.open(
    icechunk.gcs_storage(bucket=BUCKET, prefix=PREFIX, anonymous=True),
    config=config,
    authorize_virtual_chunk_access=icechunk.containers_credentials({PIXELS: None}),
)
store = repo.readonly_session("main").store

geoms = open_geometries(store, variables=["flux"], sample_axis=0)
ds = InSituDataset(
    store,
    split_by_chunk(geoms["flux"], fractions=(0.8, 0.1, 0.1)),
    geometries=geoms,
    batch_size=64,
    max_inflight=32,     # concurrent ranged GETs -- the knob that matters here
)
for batch in ds.train:
    x = batch.arrays["flux"]          # (64, 3841) float32
```

A full pass over all 3840 spectra takes roughly **7 s** from a VM in `us-central1` (59 MB in
3840 ranged reads — latency-bound, so concurrency around 32 is the knee; higher does not
help). From outside Google Cloud, expect slower.

If you would rather not install a loader, plain zarr slicing (`flux[0:640]`) is a fine way to
pull a subset — but raise zarr's own concurrency first by wrapping the reads in
`with zarr.config.set({"async.concurrency": 32}):`. It defaults to 10, which costs about a
third of the throughput on a store chunked one fiber at a time (17.9 s vs 12.2 s for the full
pass; insitubatch does it in 7.3 s, medians of five interleaved runs).

Treat that last comparison as parity rather than a ranking: the only thing tuned on the zarr
side is concurrency, and slice width, access pattern and thread-pool settings could all narrow
it. At one fiber per chunk there is nothing to amortise within a chunk, so no reader has a
structural advantage — what a batch loader adds here is shuffling, splits and transforms, not
raw throughput.

Note the examples above pass a bare `icechunk.http_store()`. Handing it a tuned connection
pool is easy to get wrong: capping `pool_max_idle_per_host` below your read concurrency makes
a sustained reader reconnect constantly, which measured 12% slower than the default.

### The other layout is much faster per spectrum

One fiber per chunk buys archive scale and costs one request per spectrum. For comparison,
plate 0266 built both ways from the same bytes — same 640 spectra, byte-identical output — is
**0.11 s at 64 fibers/chunk against 1.26 s at one fiber/chunk**, an 11.5× difference from the
chunk geometry alone. Use `sdss_dr17_p0266_refs` if throughput per spectrum matters more to
you than corpus size.

## Rebuilding, or extending to more plates

These stores were built with `examples/sdss/data.py` in
[insitubatch](https://github.com/emfdavid/insitubatch), which reads `spPlate` FITS
with VirtualiZarr and writes the Icechunk references. Point it at more plate URLs to
index a larger slice of the archive — the reference store stays kilobytes.

Correct handling of FITS big-endian data through the virtual chain needs
**VirtualiZarr ≥ 2.7.2** and **zarr ≥ 3.3**.