Skip to content

Enable the reading of tile compressed FITS files with kerchunk #14575

Description

@Cadair

kerchunk is a project which parses data and generates a zarr array using fsspec. kerchunk currently supports most FITS files, but not any tile compressed data.

This issue is intended to be a tracking issue for all the things that need to be done in all the places they need to be done to make this work.

Changes in Astropy

The main functionality that astropy needs to provide are numcodecs Codec classes for the decompression and dequantization post processing steps that need to be applied to the data on disk after read. The initial version of these Codec classes are implemented in the new _tiled_compression.codecs subpackge. Currently these codecs have not been tested outside of their internal use in io.fits.

  • Raise an error if more bytes are passed to decode than expected. Currently this segfaults.
  • Modify the codecs to accept the buf= argument.
  • Do not assume that the input to decode is already an array, when it's likely a byte buffer.
  • Refactor decompress_hdu_section so that its parts could be used by kerchunk when generating references to each tile in the heap.
  • Modify the codecs to zero-pad tiles which are smaller than the global tile size. see this comment.
  • See if it is possible to have the codecs not have per-tile settings at all. Currently PLIO, RICE and HCOMPRESS all need to know the actual tile shape of each tile (rather than the global tile shape).
  • Move finalize_array into the codec .decode() methods.

Changes not in Astropy

  • The Quantize codec requires per-tile arguments zscale and zzero, which is not supported by numcodecs / zarr.
  • We might need the decompression codecs to know the expected output size (real tile shape) for each tile.
  • Implement building the references for each tile in the FITS file. Things to consider here:
    • Some tiles in a file can be gzipped when others are compressed with other codecs
    • Data can be quantized

Stretch Goals

  • Enable combining more than one FITS tile into a single zarr chunk. This would be possible if the tiles are contiguous in the heap (not guaranteed), but the codecs would need to be modified to accept the bytes for more than one tile into .decode at once.

Activity

  1. Cadair commented on Mar 24, 2023

    @Cadair
    MemberAuthor

    @martindurant I have some questions for you based on the above list:

    1. What are the input and output types that we should support in .decode? You say that the input "is likely a byte buffer", but it can also be an array? what else is it allowed to be?
    2. How do we move forward on solving the per-tile inputs to the Quantize codec?
    3. What should the behaviour of the buf= argument be? Is this documented anywhere?
    4. What should the output of .decode be? Is it ok for it to always be a numpy array with a specific dtype as calculated by the codec? - The reason for this question is currently we do this post-processing after the call to decode but I suspect it should be inside decode?
  2. martindurant commented on Mar 24, 2023

    @martindurant

    Thanks for putting together the summary, I think you've covered all the major points.

    Some answers:

    1. buf is "buffer-like" according to the docstring, which would include numpy arrays and memoryviews. In kerchunk's usage it will be bytes. The output is something you can pass to another decode() function, or to np.frombuffer
    2. do I understand that currently the arguments for each tile are stored in the bintable of the HDU? I see two possibilities, both of which have some preexisting conversation. There may be other things, but this is the most complicate step.
    • each kerchunk reference has additional fields, and at load time these are added in the proposed "context" dict, so that the codec which is called after fetching bytes can use them
    • each kerchunk key is made of two kerchunk entries to be concatenated, one small binary chunk containing the arguments and one bytes range from the heap, so that the codec can extract all it needs from the bytes it receives.
    1. didn't understand the question
    2. zarr will coerce the output of decode to the dtype and shape expected, if not providing an explicit out= agument. So providing buffer or numpy output are both fine.
  3. Cadair commented on Mar 24, 2023

    @Cadair
    MemberAuthor

    didn't understand the question

    Sorry I mistyped, I meant what's the behaviour of the out= argument.

    do I understand that currently the arguments for each tile are stored in the bintable of the HDU?

    Yes, that's correct.

    the proposed "context" dict

    Is this a zarr proposal? and is context proposed to be an argument to decode?

  4. martindurant commented on Mar 24, 2023

    @martindurant

    Is this a zarr proposal? and is context proposed to be an argument to decode?

    This is an idea within zarr, yes. It is not a formal ZEP-level discussion, because it is an implementation detail rather than any change in the spec. The thread closest to a resolution is this one: zarr-developers/zarr-python#1131 , which was originally about defining what class of array a codec makes (numpy versus cupy).

    what's the behaviour of the out= argument

    According to the documentation, this can also be any mutable buffer (memoryview, bytearray...), but in practice I think it's always numpy or numpy compatible, so you can [:] assign into it and expect it to have the correct dtype and shape for the chunk. It provides for zero-copy in the case of decompression that accepts an output pointer.

  5. moved this from Todo to Backlog in Cadair's Todo Liston Jun 8, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions