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
96 stars 62 forks source link

Enable convolution of model cubes #766

Open keflavich opened 2 years ago

keflavich commented 2 years ago

It is sometimes necessary to convolve model cubes to a target resolution, but they don't have beams so convolution fails:

mod.convolve_to(psf.beams[0])
---------------------------------------------------------------------------
NoBeamError                               Traceback (most recent call last)
<ipython-input-84-510cae9ff8ae> in <module>
----> 1 mod.convolve_to(psf.beams[0])

/blue/adamginsburg/adamginsburg/casa/casa-6.4.0-16/lib/py/lib/python3.8/site-packages/spectral_cube/dask_spectral_cube.py in wrapper(self, *args, **kwargs)
     75     def wrapper(self, *args, **kwargs):
     76         save_to_tmp_dir = kwargs.pop('save_to_tmp_dir', False)
---> 77         cube = function(self, *args, **kwargs)
     78         if save_to_tmp_dir and isinstance(cube, DaskSpectralCubeMixin):
     79             if not ZARR_INSTALLED:

/blue/adamginsburg/adamginsburg/casa/casa-6.4.0-16/lib/py/lib/python3.8/site-packages/spectral_cube/dask_spectral_cube.py in convolve_to(self, beam, convolve, **kwargs)
   1381
   1382         # Check if the beams are the same.
-> 1383         if beam == self.beam:
   1384             warnings.warn("The given beam is identical to the current beam. "
   1385                           "Skipping convolution.")

/blue/adamginsburg/adamginsburg/casa/casa-6.4.0-16/lib/py/lib/python3.8/site-packages/spectral_cube/base_class.py in beam(self)
    818     def beam(self):
    819         if self._beam is None:
--> 820             raise NoBeamError("No beam is defined for this SpectralCube or the"
    821                               " beam information could not be parsed from the"
    822                               " header. A `~radio_beam.Beam` object can be"

NoBeamError: No beam is defined for this SpectralCube or the beam information could not be parsed from the header. A `~radio_beam.Beam` object can be added using `cube.with_beam`.

This can be worked around with some hackery that should be formalized:

delta = radio_beam.Beam(0*u.deg)
mod.beam = delta
mod.convolve_to(beam)

However, what if we want to convolve each plane to a different resolution? There's not even a workaround for that right now:

mod.convolve_to(beams) # won't work
mod.spatial_smooth(beams) # also wont, expects a single kernel

I dug a bit deeper and there are no obvious, easy ways to achieve this by writing special functions to pass to apply_function_parallel_spatial. I think we need to special-case this.

keflavich commented 2 years ago

I suggest we should make a special-case radio_beam.DeltaBeam() in radio_beam to make the Beam(0) make more sense to users.

e-koch commented 2 years ago

For the single beam case:

mod_wbeam = mod.with_beam(delta)

should work, I think

e-koch commented 2 years ago

We have with_beams in: https://github.com/radio-astro-tools/spectral-cube/blob/c3a9eceb5ace5cab774fe4fce58a4a40c2b1b024/spectral_cube/base_class.py#L467

But we don't have a SC -> VSRC method.

Agreed on a DeltaBeam to avoid confusion.

keflavich commented 2 years ago

right, it should - but, we should also probably initialize CASA models with Delta beams if we can find a way to uniquely identify them

e-koch commented 2 years ago

agreed. Is there something in the header to identify a clean model?

keflavich commented 2 years ago

the header unit is Jy pix-1 - can that reliably be treated as a Delta beam? I think yes?

keflavich commented 2 years ago

I thought about this more, and I don't think there is so much reason to add this functionality. For what I'm doing, it would be better just to find a commonbeam and smooth to that anyway.

Nonetheless, the ability should be there. That would enable things like per-plane-beam-merging with uvcombine, which we (and CASA) don't support right now. @urvashirau mentioned needing this coincidentally (?) today.

keflavich commented 2 years ago

https://github.com/ALMA-IMF/reduction/commit/9d3846ff08e80a8c1fe1771d4ff87c3c628ead78#diff-91c493e486904d56f624b0744a0257ec2603050401cdb71aaa9d535316ecc8adR118 shows a key use case we should make sure is well-supported

keflavich commented 2 years ago

https://github.com/ALMA-IMF/reduction/blob/master/misc/beam_volume_tools.py