pixell cannot load CompImageHDU

[mhasself]

• pixell version: 0.8.0+17
• Python version: 3.8
• Operating System: Linux

Description

pixell can't load images that are saved into compressed HDUs of FITS files, even though this is supported by astropy. (The example file below comes from SPT folks so it's a real thing CMB people use.)

This is because enmap.ndmap_proxy_fits uses an API specific to uncompressed HDUs (for reasons of I/O and RAM efficiency in the case where only a small submap should ultimately be loaded from disk). A new proxy class will be needed to access CompImageHDU objects (or plausibly any HDU types that do not expose the raw API).

What I did

Unzip the attached FITS file and try to load it, e.g. like:

m = enmap.read_fits('compressed_zeros_TQU_zea.fits', hdu=1)
Traceback (most recent call last):
  File "/home/mhasse/.local/lib/python3.8/site-packages/pixell/enmap.py", line 1958, in read_fits
    proxy = ndmap_proxy_fits(hdu, wcs, fname=fname, threshold=sel_threshold)
  File "/home/mhasse/.local/lib/python3.8/site-packages/pixell/enmap.py", line 2067, in __init__
    dtype    = fix_endian(hdu.section[(slice(0,1),)*hdu.header["NAXIS"]]).dtype
AttributeError: 'CompImageHDU' object has no attribute 'section'

Here's the demo file: compressed_zeros_TQU_zea.fits.zip

[msyriac]

+1 because I was about to raise this issue a few days ago when working with Matt's files, but never actually got around to it.

[msyriac]

This issue has resurfaced. Here is the output from attempting to load the public DR5 cluster survey mask on LAMBDA:

In [2]: a = enmap.read_map('DR5_cluster-search-area-mask_v1.0.fits')
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-2-6d2a772965be> in <module>
----> 1 a = enmap.read_map('DR5_cluster-search-area-mask_v1.0.fits')

~/repos/pixell/pixell/enmap.py in read_map(fname, fmt, sel, box, pixbox, geometry, wrap, mode, sel_threshold, wcs, hdu, delayed, verbose)
   2082                 else: fmt = "fits"
   2083         if fmt == "fits":
-> 2084                 res = read_fits(fname, sel=sel, box=box, pixbox=pixbox, geometry=geometry, wrap=wrap, mode=mode, sel_threshold=sel_threshold, wcs=wcs, hdu=hdu, delayed=delayed, verbose=verbose)
   2085         elif fmt == "hdf":
   2086                 res = read_hdf(fname, sel=sel, box=box, pixbox=pixbox, geometry=geometry, wrap=wrap, mode=mode, sel_threshold=sel_threshold, wcs=wcs, delayed=delayed)

~/repos/pixell/pixell/enmap.py in read_fits(fname, hdu, sel, box, pixbox, geometry, wrap, mode, sel_threshold, wcs, delayed, verbose)
   2172         ndim = len(hdu.shape)
   2173         if hdu.header["NAXIS"] < 2:
-> 2174                 raise ValueError("%s is not an enmap (only %d axes)" % (fname, hdu.header["NAXIS"]))
   2175         if wcs is None:
   2176                 with warnings.catch_warnings():

ValueError: DR5_cluster-search-area-mask_v1.0.fits is not an enmap (only 0 axes)

The file is here: https://lambda.gsfc.nasa.gov/product/act/actpol_dr5_szcluster_catalog_get.cfm

[msyriac]

Pinging @mattyowl and @mhasself for help here

[mhasself]

Seems to work if you point it to second hdu:

enmap.read_fits('DR5_cluster-search-area-mask_v1.0.fits', hdu=1)

Not sure if a compressed image (as this is) needs to live at hdu>0 in order to be a binary table?

[mhasself]

Indeed a compressed image cannot be a PrimaryHDU (see here). A reasonable behavior for read_fits, when hdu=None, might be for it to skip "empty" (for some definition of empty) HDUs until it finds a good one?

[mattyowl]

Sorry I missed this (but have sorted out my GitHub mail filtering) - I agree with what @mhasself said above (spinning through hdus to find the first non-empty one, if hdu is not given).