Skip to content

Reimplement the FITS parser natively, without kerchunk - #1090

Draft
TomNicholas wants to merge 1 commit into
zarr-developers:mainfrom
TomNicholas:native-fits-parser
Draft

Reimplement the FITS parser natively, without kerchunk#1090
TomNicholas wants to merge 1 commit into
zarr-developers:mainfrom
TomNicholas:native-fits-parser

Conversation

@TomNicholas

@TomNicholas TomNicholas commented Aug 17, 2026

Copy link
Copy Markdown
Member

Parses FITS with astropy directly and builds the ManifestStore from the HDU offsets it reports, instead of routing through kerchunk's reference dict and translating it back. The fits extra now installs only astropy.

Follows the same pattern as the native netCDF3 parser: read the format's own structure, build ManifestArrays from it, skip the round-trip through kerchunk refs entirely.

Coverage went up

ASCII tables (TABLE HDUs) previously could not be read at all — get_codec_class("numcodecs.FITSAscii") raised KeyError, so they failed at metadata construction. This adds a real zarr v3 FITSAsciiTableCodec. The Hubble file used by the existing test turned out to contain an ASCII table that kerchunk was silently skipping; all six of its columns now match astropy exactly.

Three correctness bugs found along the way:

  • BSCALE was inverted. It was passed straight to FixedScaleOffset, which divides on decode, so every scaled image came back scaled by 1/BSCALE.
  • ASCII column layout was wrong. The previous approach assumed columns tile the row with a one-byte separator. They don't — TBCOL places them, and the real Hubble table has gaps (offsets [0, 16, 32, 48, 64, 92], itemsize 104). Parsing also has to allow the Fortran D exponents the format uses, which numpy cannot parse; without that the real file fails outright.
  • BLANK was ignored, despite custom_parsers.md documenting that FITS parsers should map it. It now becomes a _FillValue, scaled alongside the data when BSCALE/BZERO are present.

Behaviour changes

  • Every data HDU is exposed, not just the first.
  • Axes are named per-HDU (SCI_y, SCI_x) rather than a shared y/x. The shared names only worked because kerchunk emitted a single extension; with several HDUs they collide, since HDUs routinely disagree about the length of a given axis.
  • reader_options is gone. The parser reads through the ObjectStoreRegistry rather than fsspec, so fsspec storage_options no longer apply — configure credentials on the store you register.

Binary tables now raise

Their columns form a big-endian structured dtype, and the zarr v3 struct data type has no way to record a field's byte order: [('flux', '>f4')] serializes to plain float32 and reparses little-endian, and zarr's bytes codec does not apply endian to a structured dtype anyway. So any virtual reference to big-endian record data reads back byte-swapped once the metadata has been through zarr.json.

This is pre-existing and not specific to FITS — it reproduces with a hand-built ManifestArray and no FITS involved. The parser previously accepted these HDUs and returned wrong numbers; it now refuses them, and skip_variables excludes them. Worth fixing properly at the zarr level, but that is a separate change with a much wider blast radius.

Testing

Local fixtures cover images, cubes, scaled images, ASCII tables, BLANK masking, per-HDU axis naming, duplicate EXTNAME disambiguation, and the binary-table refusal, so the parser is no longer only exercised by a network test. The two network tests validate against the real Hubble file, including the previously-skipped ASCII table.

The existing Hubble test passed region="us-west-2", but the bucket is us-east-1. It went unnoticed because kerchunk read via fsspec and the registered store was never actually used; now that the parser reads through it, the region has to be right.

Parse FITS with astropy directly and build the ManifestStore from the HDU
offsets it reports, rather than routing through kerchunk's reference dict.
The fits extra now installs only astropy, and the parser reads through the
ObjectStoreRegistry instead of fsspec, so reader_options is gone.

Expose every data HDU rather than only the first, and name each HDU's axes
after that HDU: separate HDUs routinely disagree about the length of a given
axis, which a shared y/x cannot express.

ASCII tables previously could not be read at all, since nothing was
registered under the codec name kerchunk emitted for them. Add a real zarr
v3 codec, and locate columns by TBCOL rather than assuming they tile the row
with a separator between them, which is not how the format works. Parsing
also has to allow the Fortran D exponents and blank fields the format uses.

Fix BSCALE being passed straight to FixedScaleOffset, which divides rather
than multiplies on decode, so scaled images came back scaled by 1/BSCALE.
Map an integer image's BLANK keyword to a _FillValue attribute, scaled
alongside the data when BSCALE or BZERO is present.

Refuse binary tables. Their columns form a big-endian structured dtype, and
the zarr v3 struct data type cannot record a field's byte order, so the
values would read back byte-swapped once the metadata has been serialized.
Returning nothing is better than returning wrong numbers.

The Hubble test file turned out to contain an ASCII table that was being
silently skipped, and to live in us-east-1 rather than us-west-2 -- the
wrong region went unnoticed while kerchunk did the reading through fsspec
and the registered store was never used.
@TomNicholas TomNicholas added the references generation Reading byte ranges from archival files label Aug 17, 2026
@codecov

codecov Bot commented Aug 17, 2026

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 96.96970% with 4 lines in your changes missing coverage. Please review.
✅ Project coverage is 90.04%. Comparing base (7805559) to head (ce3d802).
⚠️ Report is 1 commits behind head on main.

Files with missing lines Patch % Lines
virtualizarr/codecs.py 95.34% 2 Missing ⚠️
virtualizarr/parsers/fits.py 97.75% 2 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main    #1090      +/-   ##
==========================================
+ Coverage   89.67%   90.04%   +0.37%     
==========================================
  Files          41       41              
  Lines        2682     2803     +121     
==========================================
+ Hits         2405     2524     +119     
- Misses        277      279       +2     
Files with missing lines Coverage Δ
virtualizarr/codecs.py 94.78% <95.34%> (+0.26%) ⬆️
virtualizarr/parsers/fits.py 97.89% <97.75%> (-2.11%) ⬇️

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@emfdavid

emfdavid commented Sep 2, 2026

Copy link
Copy Markdown
Contributor

This looks good -- the native parser is nice, and the BSCALE/TBCOL/BLANK catches are real.

On the binary tables: the byte-order problem is fixed. I reported it as zarr-python#4141 and @d-v-b fixed it in zarr-python#4142, first released in zarr 3.3.0 (2026-07-30). VirtualiZarr#1047/#1048 are the write-side half, released in v2.7.2 (2026-08-05), so convert_to_codec_pipeline already emits endian: big for a big-endian struct. Both landed before this PR was opened.

Main reads real binary tables correctly today. On SDSS DR17 spZbest-0266-51602.fits restricted to its 49 scalar columns, VZ main + zarr 3.3.0 gives 49/49 columns matching astropy, with {'name': 'bytes', 'configuration': {'endian': 'big'}} in the metadata:

codecs: ({'name': 'bytes', 'configuration': {'endian': 'big'}},)
columns matching astropy: 49/49
delivered byteorder: < (native)
  PLATE [266 266 266]  Z [0.02127545 0.21392463 0.65218145]

So the refusal drops a case that works today.

What genuinely doesn't work is array-valued columns: the v3 struct type cannot express a subarray field, so ('>f4', (5,)) degrades to opaque raw_bytes and the codec cannot swap what it can no longer see. Today it surfaces as a bare TypeError: a bytes-like object is required, not 'int' out of Struct.default_scalar(), which is a poor error for a FITS user to land on -- and it is pre-existing, not something this PR introduced. I don't think there is a zarr issue open for it; happy to file one. For now, if you set a fill_value - the issue is bypassed and you get raw bytes.

metadata = create_v3_array_metadata(
    shape=(nrows,), data_type=stored_dtype, chunk_shape=(nrows,),
    fill_value=np.zeros(1, stored_dtype)[0],   # skips Struct.default_scalar()
    ...
)

Suggestion: gate on zarr>=3.3.0 at runtime and parse all the fixed-width columns. The scalars come back correct — 49/49 on the real file. The fixed-length array columns degrade to raw_bytes, which is honest: the bytes are exact, and .view('>i4') recovers them, so worth an attribute or a docstring note rather than a refusal till it can be fixed in zarr. Variable-length columns (TFORM P/Q) are the one case I'd still refuse.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

references generation Reading byte ranges from archival files

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants