BUNIT='beam-1 Jy' conundrum.

[spectram]

Here from deep within the rabbit hole, I thought I'd raise an issue to keep track of the cascade of other issues raised related to this.

1. Issue raised in Astropy: Regarding FITS standard definition for 'Jy/beam' astropy/astropy#13217
2. Issue raised in SoFiA2: https://github.com/SoFiA-Admin/SoFiA-2/issues/74
3. Issue raised in SpectralCube: Question about writing BUNIT as 'Jy/beam' instead of 'beam-1 Jy' radio-astro-tools/spectral-cube#806

If nothing gets fixed I'd recommend doing something like:
if BUNIT == 'beam-1 Jy:
BUNIT='Jy/beam'
fits.write()
print('please re-run SoFiA2 and come back for accurate NHI, W20, W50)

The NHI calculation can easily be fixed by the addition of BUNIT == 'beam-1 Jy' in the current if statement but W20 and W50 values from the SoFiA2 catalogue file are in 'pix'.

[kmhess]

Note, that there is still something weird with how the w20/w50 is written in the mom1 maps (it's 0!). I don't yet know if it's a problem in SIP or SoFiA, but I'm noting it here so I can look into it again when I actually have time.