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
98 stars 65 forks source link

Beam dimensions incorrectly assumed to be in arcsec for VaryingResolutionSpectralCube #737

Closed drtobybrown closed 3 years ago

drtobybrown commented 3 years ago

When reading a cube with varying spectral resolution where BMAJ, BMIN header keyword values are in degrees (which I'd thought was standard) spectral-cube assumes they're in arcsec.

I believe the assumption is at line 3654 in spectral_cube.py

# CASA beam tables are in arcsec, and that's what we support
beams = Beams(major=u.Quantity(beam_data_table['BMAJ'], u.arcsec),
              minor=u.Quantity(beam_data_table['BMIN'], u.arcsec),
              pa=u.Quantity(beam_data_table['BPA'], u.deg),
              meta=[{key: row[key] for key in beam_data_table.names
                     if key not in ('BMAJ','BPA', 'BMIN')}
                    for row in beam_data_table],
             )

I'm not familiar enough with the code to find where the beam_data_table (or beam_table) are parsed from the headers. Beam tables in arcsec seems ok, but the header values are not being converted from degrees.

If it helps, cube I'm using is here (updated link)

keflavich commented 3 years ago

We'll have to look more closely at that cube to see what the problem is. That __init__ code is used for FITS reading, but not for CASA .image reading. I've also verified that CASA generally uses arcseconds, not degrees, in the BMIN/BMAJ columns. Note that the BMIN/BMAJ columns are not equivalent to the BMIN/BMAJ header keywords, which are in degrees by default.

for CASA images, https://github.com/radio-astro-tools/spectral-cube/blob/master/spectral_cube/io/casa_image.py#L101-L126

    if 'major' in beam_:
        beam = Beam(major=u.Quantity(beam_['major']['value'], unit=beam_['major']['unit']),
                    minor=u.Quantity(beam_['minor']['value'], unit=beam_['minor']['unit']),
                    pa=u.Quantity(beam_['positionangle']['value'], unit=beam_['positionangle']['unit']),
                   )
    elif 'beams' in beam_:
        bdict = beam_['beams']
        if beam_['nStokes'] > 1:
            raise NotImplementedError()
        nbeams = len(bdict)
        assert nbeams == beam_['nChannels']
        stokesidx = '*0'

        majors = [u.Quantity(bdict['*{0}'.format(ii)][stokesidx]['major']['value'],
                             bdict['*{0}'.format(ii)][stokesidx]['major']['unit']) for ii in range(nbeams)]
        minors = [u.Quantity(bdict['*{0}'.format(ii)][stokesidx]['minor']['value'],
                             bdict['*{0}'.format(ii)][stokesidx]['minor']['unit']) for ii in range(nbeams)]
        pas = [u.Quantity(bdict['*{0}'.format(ii)][stokesidx]['positionangle']['value'],
                          bdict['*{0}'.format(ii)][stokesidx]['positionangle']['unit']) for ii in range(nbeams)]

        beams = Beams(major=u.Quantity(majors),
                      minor=u.Quantity(minors),
                      pa=u.Quantity(pas))
    else:
        warnings.warn("No beam information found in CASA image.",
                      BeamWarning)
drtobybrown commented 3 years ago

We've been having some issues with header in this cube as it was made with AIPS so maybe that's it. I've updated the link in the original post.. should work now

keflavich commented 3 years ago

Got it. We've never seen one of these before, but this is easy enough to deal with: the header actually says what the units are, and we should be using those units instead of hard-coding them

keflavich commented 3 years ago

@drtobybrown see #738 for the fix. Test it out if you can.

drtobybrown commented 3 years ago

that works for me, thanks