radio-astro-tools / spectral-cube

Library for reading and analyzing astrophysical spectral data cubes
http://spectral-cube.rtfd.org
BSD 3-Clause "New" or "Revised" License
97 stars 65 forks source link

Undefined variable 'key' prevents cubes from being opened #785

Open marcinglowacki opened 2 years ago

marcinglowacki commented 2 years ago

Attempting to load in a new spectral-line cube with SpectralCube.read (previous cubes worked fine), but received the following error for this one:

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-6-5d60b59ea817> in <module>
----> 1 SpectralCube.read('/idia/projects/laduma/DR1.0/laduma_dr1.0_image.1380~1420MHz.fits')

/anaconda3/lib/python3.6/site-packages/spectral_cube/io/core.py in __call__(self, filename, *args, **kwargs)
    116         kwargs['target_cls'] = BaseSpectralCube
    117         try:
--> 118             return registry.read(BaseSpectralCube, filename, *args, **kwargs)
    119         except IsADirectoryError:  # See note above StringWrapper
    120             return registry.read(BaseSpectralCube, StringWrapper(filename), *args, **kwargs)

/anaconda3/lib/python3.6/site-packages/astropy/io/registry.py in read(cls, format, cache, *args, **kwargs)
    518 
    519         reader = get_reader(format, cls)
--> 520         data = reader(*args, **kwargs)
    521 
    522         if not isinstance(data, cls):

/anaconda3/lib/python3.6/site-packages/spectral_cube/io/fits.py in load_fits_cube(input, hdu, meta, target_cls, use_dask, **kwargs)
    186         VRSC = VaryingResolutionSpectralCube
    187 
--> 188     data, header, beam_table, beam_units = read_data_fits(input, hdu=hdu, **kwargs)
    189 
    190     if data is None:

/anaconda3/lib/python3.6/site-packages/spectral_cube/io/fits.py in read_data_fits(input, hdu, mode, **kwargs)
    160 
    161         with fits_open(input, mode=mode, **kwargs) as hdulist:
--> 162             return read_data_fits(hdulist, hdu=hdu)
    163 
    164     return array_hdu.data, array_hdu.header, beam_table, beam_units

/anaconda3/lib/python3.6/site-packages/spectral_cube/io/fits.py in read_data_fits(input, hdu, mode, **kwargs)
    104                     for i in range(1, 4):
    105                         if not f"TUNIT{i}" in hdu_item.header:
--> 106                             raise BeamUnitsError(f"Missing beam units keyword {key}{i}"
    107                                                     " in the header.")
    108 

NameError: name 'key' is not defined

This variable does not seem to be defined in the code; perhaps one easy fix is changing line 104 to for i, key in zip(range(1, 4), ["BMAJ", "BMIN", "BPA"])

keflavich commented 2 years ago

It looks like that should actually be modified to read 'TUNIT{i}'. We'd welcome a pull request if you'd like to make that change yourself.

That won't help with your issue, though - it looks like your table has an incorrectly specified header. It has 'BPA' in the column names, but does not have TUNIT specified for all of the columns. Could you tell us a little more about the origin of these data? What telescope, what reduction software produced them? We may need to add a special case for these.

keflavich commented 2 years ago

This is affecting some of my old cubes too:

cube = SpectralCube.read('merge/fullcubes/full_SgrB2_TETC7m_r2_spw0_lines.fits')
Traceback (most recent call last):
  File "<ipython-input-5-d7489d58e7fb>", line 1, in <module>
    cube = SpectralCube.read('merge/fullcubes/full_SgrB2_TETC7m_r2_spw0_lines.fits')
  File "/blue/adamginsburg/adamginsburg/repos/spectral-cube/spectral_cube/io/core.py", line 118, in __call__
    return registry.read(BaseSpectralCube, filename, *args, **kwargs)
  File "/orange/adamginsburg/miniconda3/envs/python39/lib/python3.9/site-packages/astropy/io/registry.py", line 527, in read
    data = reader(*args, **kwargs)
  File "/blue/adamginsburg/adamginsburg/repos/spectral-cube/spectral_cube/io/fits.py", line 188, in load_fits_cube
    data, header, beam_table, beam_units = read_data_fits(input, hdu=hdu, **kwargs)
  File "/blue/adamginsburg/adamginsburg/repos/spectral-cube/spectral_cube/io/fits.py", line 162, in read_data_fits
    return read_data_fits(hdulist, hdu=hdu)
  File "/blue/adamginsburg/adamginsburg/repos/spectral-cube/spectral_cube/io/fits.py", line 106, in read_data_fits
    raise BeamUnitsError(f"Missing beam units keyword {key}{i}"
NameError: name 'key' is not defined
keflavich commented 2 years ago

This error is masking a deeper one: old versions of CASA wrote Beam tables with no units, e.g.:

In [15]: fits.open('merge/fullcubes/full_SgrB2_TETC7m_r2_spw0_lines.fits')[1].header
Out[15]:
XTENSION= 'BINTABLE'           / binary table extension
BITPIX  =                    8 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                   20 / length of dimension 1
NAXIS2  =                 7680 / length of dimension 2
PCOUNT  =                    0 / number of group parameters
GCOUNT  =                    1 / number of groups
TFIELDS =                    5 / number of table fields
TTYPE1  = 'BMAJ    '
TFORM1  = 'E       '
TTYPE2  = 'BMIN    '
TFORM2  = 'E       '
TTYPE3  = 'BPA     '
TFORM3  = 'E       '
TTYPE4  = 'CHAN    '
TFORM4  = 'J       '
TTYPE5  = 'POL     '
TFORM5  = 'J       '

We should assume default units of deg, right? (@e-koch ?)

e-koch commented 2 years ago

I'm unsure about the default units. I don't have a cube where the units weren't defined, but the ones I do have are in arcsec.

marcinglowacki commented 2 years ago

Late reply: thank you for the response!

This particular cube came from image-plane combination via a custom script, and thanks to your point we've now fixed the header and gotten it working. Other cubes we've been working with related to this (MeerKAT, reduced via CASA/CASA-based pipelines, one of the 3. versions) have worked fine with Spectral-cube.

keflavich commented 2 years ago

Thanks @marcinglowacki , that's helpful.

@e-koch agreed, it's arcsecs by default

keflavich commented 2 years ago

My bad beam table appears to have been produced by ORIGIN = 'CASA 4.7.2-REL (r39762)'

e-koch commented 2 years ago

We should add a test of reading in a VRSC written by spectral-cube to ensure the units are written properly to the beams table

Similar to adding a units check to this 1D spectrum example, but with a read/write of the whole cube https://github.com/radio-astro-tools/spectral-cube/blob/master/spectral_cube/tests/test_io.py#L116