SoFiA-2
SoFiA-2 copied to clipboard
SoFiA-Admin
Dealing with BUNIT='beam-1 Jy'
SpectralCube.write() saves BUNIT as 'beam-1 Jy' because of Astropy's FITS standard definition. SoFiA2 fails to accurately calculate W20, and W50. https://github.com/radio-astro-tools/spectral-cube/issues/806
Thank you for reporting this problem. Indeed, SoFiA currently only checks for Jy/beam and common misspellings such as JY/BEAM or Jy/Beam, as this is the standard format used by most software pipelines. Allowing any possible way of expressing the same unit, e.g. beam**-1 Jy or 10^-26*W*m^-2*Hz^-1*beam^-1 and so forth, would require writing a complex physical unit parser that can determine whether a specific string decodes a unit of the desired dimension according to the rules set up by the FITS standard. This is a non-trivial task and not something that we would be able to implement in the short term. If you are aware of an open-source implementation of such a FITS unit string parser in C or C++, please let me know so I can check whether we could incorporate that into SoFiA.
For the time being, a simple workaround might be to simply change the value of the BUNIT keyword to Jy/beam (e.g. with Astropy) before loading the FITS file into SoFiA.
You also mentioned that SoFiA fails to accurately calculate w₂₀ and w₅₀ if BUNIT is not Jy/beam. Could you clarify what exactly you mean by that? I would have expected the calculation of line widths to be independent of the flux density unit, so their value shouldn’t change no matter which BUNIT is defined, but if that’s not the case then this could indicate that there is a bug in the line width measurement algorithm in SoFiA.
The problem reported arises because of the unusual realisation of the unit Jy/beam by SpectralCube, so indeed the problem should be fixed there and it is good to see that two issues have been raised about this: https://github.com/radio-astro-tools/spectral-cube/issues/806 and https://github.com/radio-astro-tools/spectral-cube/issues/794. In fact, as far as I am aware all major radio astronomy packages (CASA, AIPS, MIRIAD) and radio observatories use Jy/beam (sometimes JY/BEAM) in their FITS headers.
Indeed it is best for now to change the unit to Jy/beam in the offending file using the fits routines in astropy and wait for the answer to the two issues posted on https://github.com/radio-astro-tools/spectral-cube/
It may be good to link issue 806 to issue 794 and remind the authors that none of the radio observatories/software packages use beam-1 Jy as unit. By the way, the IAU has Jy as recognised unit, but does not stipulate anything about Jy/beam (https://www.iau.org/publications/proceedings_rules/units/).
I am currently manually changing the BUNIT before running SoFiA2. I've raised the issue with astropy (https://github.com/astropy/astropy/issues/13217)
SoFiA2 propagates the unit just fine (All units in file headers remain the same as BUNIT). The W20, and W50 values are given in pix instead of m/s even with parameter.physical = True. My uninformed guess would be a parsing error because of the space after beam-1 (this was the case with SIP)
@spectram: Thank you for clarifying the w20 and w50 issue. If I understand this correctly, then the line widths are still correct, but specified in channels rather than physical frequency axis units such as Hz or km/s. This is indeed something that we are aware of and intend to change such that w20 and w50 are always specified in physical units whenever parameter.wcs = true (see open issue #63). That would decouple the line widths from the BUNIT issue altogether, which is probably closer to what a SoFiA user would expect.
I think I noticed that W20 and W50 values were < 0.1 (order of magnitude 10^-2). Although I think decoupling BUNIT should fix this as well. Shall I close this then?
We discussed this issue at our recent SoFiA developer‘s meeting this week and have decided to implement a check for a few additional strings that correspond to the same unit of Jy/beam. If one of these is found, we would then convert BUNIT to the commonly-used standard of Jy/beam, which should fix the issue at least for SpectralCube users. Other strings could then be added to that checklist in the future as needed. So please keep the issue open for now until we have implemented this change. I’d be happy to close it myself thereafter.
Posted a comment on Astropy suggesting to at least use Jy beam-1 if Jy/beam is unacceptable and referred to the use by radio observatories and their software and the FITS standards document.
It has been added to issue #13217: https://github.com/astropy/astropy/issues/13217