Flux conservation when smoothing in units of Jy/beam

[ashleythomasbarnes]

When spatially smoothing a cube in units of "Jy/beam", does spectral-spectral cube conserve flux?

e.g.

import radio_beam
from spectral_cube import SpectralCube

beam5 = radio_beam.Beam(major=5*au.arcsec, minor=5*au.arcsec, pa=0*au.deg)
inputfile='some cube.fits'
cubejy = SpectralCube.read(inputfile)
cubejy = cubejy.with_spectral_unit('km/s')
cubejy

....SpectralCube with shape=(20, 512, 512) and unit=Jy / beam...

cubejy5 = cubejy.convolve_to(beam5)

spectral-spectral seems to treat cubejy5 as though it was in units of K - i.e. not accounting for the 1/beam when smoothing. This then returns an erroneously low flux in the cube.

[keflavich]

Jy/beam is a unit of surface brightness, just like Kelvin. The problem you're observing is that the Jy/beam unit is now different: the beam size has changed.

SpectralCube presently treats this correctly; the resulting cube has the beam you've convolved to. You can check this by looking at original_cube.beam or original_cube.beams and comparing to cubejy5.beam - they should be different beams.

You can also check by converting your results to Kelvin. They should be correct - but they won't be the same in general, since you're averaging regions of different surface brightness.

[keflavich]

Ash is right here, actually - we are keeping the cube in units of the original beam; if we update the beam (and we should), we should also be scaling up the flux unit. This behavior is unique to Jy/beam scaling - for Jy/sr and Kelvin, the behavior is correct. But for Jy/beam, the beam is changing and we haven't accounted for that.

[keflavich]

@e-koch you want to look into tackling this? We need to special-case the per beam unit handling.

[e-koch]

@keflavich Will do!

[e-koch]

This issue goes a bit deeper: we don't have built in conversion handling for other surface brightness units either (e.g. Jy/beam -> Jy/sr). We can pass equivalencies for simple conversions (https://docs.astropy.org/en/stable/api/astropy.units.equivalencies.beam_angular_area.html#astropy.units.equivalencies.beam_angular_area) but checks for changes in the spatial resolution aren't handled. See #683

[e-koch]

List of operations where the unit change needs reflect the beam change:

• convolve_to
• with_beam - raise error if beam is already defined (invalid to just assign another beam); possibly add a warning telling the user the assumptions of assigning a beam when the beam was not defined before
• reproject? - Does a Jy/pix need handling in this case?
• spatial_smooth_median -- This is not well-defined with the beam.... we could not allow median smoothing when in units of Jy/beam or Jy/pix. Or have a keyword option to disable the error if the user really doesn't care about the unit implications
• spatial_smooth -- Similar to the median above because the resulting convolved kernel shape need not be Gaussian

Add if someone notices a missing operation from this list.

[keflavich]

re: reproject - yes, definitely needs jy/pix handling.

re: spatial_smooth_median: this is an interesting one, as it's better characterized as a filter rather than a smooth operation, even though the effect is a smoother image. You're right, though, default should be to say "NO" and force the user to override it.

[keflavich]

@e-koch have you solved some or all of this now?

[e-koch]

Partially fixed in #700 to properly handle more BUNIT conversions. But I haven't fixed the Jy/beam conversions when the beam changes, yet.