Cannot export data cube into yt, jy/beam is not recognized

[paulocortes]

Trying to export an ALMA data cube into yt, got the following error:

In [9]: cube = SpectralCube.read('hcopTemp_NGC6334I.I.fits')

In [10]: print(cube)
VaryingResolutionSpectralCube with shape=(51, 1500, 1500) and unit=Jy / beam:
 n_x:   1500  type_x: RA---SIN  unit_x: deg    range:   260.215825 deg:  260.228658 deg
 n_y:   1500  type_y: DEC--SIN  unit_y: deg    range:   -35.788425 deg:  -35.778015 deg
 n_s:     51  type_s: FREQ      unit_s: Hz     range: 260244487520.382 Hz:260266190479.550 Hz

In [11]: ytcube = cube.to_yt(spectral_factor=0.5)
    ...: 
yt : [WARNING  ] 2022-01-12 12:03:20,875 Cannot find time
yt : [INFO     ] 2022-01-12 12:03:20,876 Detected these axes: RA---SIN DEC--SIN FREQ 
yt : [WARNING  ] 2022-01-12 12:03:20,881 No length conversion provided. Assuming 1 = 1 cm.
yt : [WARNING  ] 2022-01-12 12:03:20,882 Assuming 1.0 = 1.0 s
yt : [WARNING  ] 2022-01-12 12:03:20,882 Assuming 1.0 = 1.0 g
yt : [INFO     ] 2022-01-12 12:03:20,932 Parameters: current_time              = 0.0
yt : [INFO     ] 2022-01-12 12:03:20,932 Parameters: domain_dimensions         = [1500 1500   51]
yt : [INFO     ] 2022-01-12 12:03:20,932 Parameters: domain_left_edge          = [0.5 0.5 0.5]
yt : [INFO     ] 2022-01-12 12:03:20,932 Parameters: domain_right_edge         = [1500.5 1500.5   26. ]
yt : [INFO     ] 2022-01-12 12:03:20,932 Parameters: cosmological_simulation   = 0
---------------------------------------------------------------------------
UnitParseError                            Traceback (most recent call last)
~/anaconda3/lib/python3.7/site-packages/spectral_cube/spectral_cube.py in to_yt(self, spectral_factor, nprocs, **kwargs)
   2245             try:
-> 2246                 ds.quan(1.0,units)
   2247             except UnitParseError:

~/anaconda3/lib/python3.7/site-packages/unyt/array.py in __new__(cls, input_scalar, units, registry, dtype, bypass_validation, input_units, name)
   2061             bypass_validation=bypass_validation,
-> 2062             name=name,
   2063         )

~/anaconda3/lib/python3.7/site-packages/unyt/array.py in __new__(cls, input_array, units, registry, dtype, bypass_validation, input_units, name)
    572             # it's a str.
--> 573             units = Unit(input_units, registry=registry)
    574 

~/anaconda3/lib/python3.7/site-packages/unyt/unit_object.py in __new__(cls, unit_expr, base_value, base_offset, dimensions, registry, latex_repr)
    264             # lookup the unit symbols
--> 265             unit_data = _get_unit_data_from_expr(unit_expr, registry.lut)
    266             base_value = unit_data[0]

~/anaconda3/lib/python3.7/site-packages/unyt/unit_object.py in _get_unit_data_from_expr(unit_expr, unit_symbol_lut)
    986         for expr in unit_expr.args:
--> 987             unit_data = _get_unit_data_from_expr(expr, unit_symbol_lut)
    988             base_value *= unit_data[0]

~/anaconda3/lib/python3.7/site-packages/unyt/unit_object.py in _get_unit_data_from_expr(unit_expr, unit_symbol_lut)
    974     if isinstance(unit_expr, Pow):
--> 975         unit_data = _get_unit_data_from_expr(unit_expr.args[0], unit_symbol_lut)
    976         power = unit_expr.args[1]

~/anaconda3/lib/python3.7/site-packages/unyt/unit_object.py in _get_unit_data_from_expr(unit_expr, unit_symbol_lut)
    971     if isinstance(unit_expr, Symbol):
--> 972         return _lookup_unit_symbol(unit_expr.name, unit_symbol_lut)
    973 

~/anaconda3/lib/python3.7/site-packages/unyt/unit_registry.py in _lookup_unit_symbol(symbol_str, unit_symbol_lut)
    330     raise UnitParseError(
--> 331         "Could not find unit symbol '%s' in the provided " "symbols." % symbol_str
    332     )

UnitParseError: Could not find unit symbol 'beam' in the provided symbols.

During handling of the above exception, another exception occurred:

RuntimeError                              Traceback (most recent call last)
<ipython-input-11-ec8e7e445f61> in <module>
----> 1 ytcube = cube.to_yt(spectral_factor=0.5)

~/anaconda3/lib/python3.7/site-packages/spectral_cube/spectral_cube.py in to_yt(self, spectral_factor, nprocs, **kwargs)
   2247             except UnitParseError:
   2248                 raise RuntimeError("The unit %s was not parsed by yt. " % units+
-> 2249                                    "Check to make sure it is correct.")
   2250 
   2251         else:

RuntimeError: The unit Jy / beam was not parsed by yt. Check to make sure it is correct.

[keflavich]

Darn! Didn't realize yt can't handle beam units. We might need to propose a fix upstream for that, or drop units when sending to yt.

[paulocortes]

Hi Adam, I did try to convert to to Kelvin, but failed when trying to visualize.

` In [16]: from radio_beam import Beam

In [17]: new_beam = Beam(0.4 * u.arcsec)

In [18]: new_cube = cube.with_beam(new_beam)

In [19]: kcube=new_cube.to(u.K)

In [21]: ytcube = kcube.to_yt(spectral_factor=0.5) yt : [WARNING ] 2022-01-12 14:36:46,551 Cannot find time yt : [INFO ] 2022-01-12 14:36:46,552 Detected these axes: RA---SIN DEC--SIN FREQ yt : [WARNING ] 2022-01-12 14:36:46,558 No length conversion provided. Assuming 1 = 1 cm. yt : [WARNING ] 2022-01-12 14:36:46,558 Assuming 1.0 = 1.0 s yt : [WARNING ] 2022-01-12 14:36:46,558 Assuming 1.0 = 1.0 g yt : [INFO ] 2022-01-12 14:36:46,602 Parameters: current_time = 0.0 yt : [INFO ] 2022-01-12 14:36:46,602 Parameters: domain_dimensions = [1500 1500 51] yt : [INFO ] 2022-01-12 14:36:46,603 Parameters: domain_left_edge = [0.5 0.5 0.5] yt : [INFO ] 2022-01-12 14:36:46,603 Parameters: domain_right_edge = [1500.5 1500.5 26. ] yt : [INFO ] 2022-01-12 14:36:46,603 Parameters: cosmological_simulation = 0`

`In [22]: images = ytcube.quick_render_movie('outdir') /home/paulo/anaconda3/lib/python3.7/site-packages/spectral_cube/spectral_cube.py:532: RuntimeWarning: All-NaN slice encountered result = function(np.dstack((result, plane)), axis=2, **kwargs) yt : [INFO ] 2022-01-12 14:38:28,858 Adding field flux to the list of fields.

AttributeError Traceback (most recent call last)

in ----> 1 images = ytcube.quick_render_movie('outdir') ~/anaconda3/lib/python3.7/site-packages/spectral_cube/ytcube.py in quick_render_movie(self, outdir, size, nframes, camera_angle, north_vector, rot_vector, colormap, cmap_range, transfer_function, start_index, image_prefix, output_filename, log_scale, rescale) 140 141 center = self.dataset.domain_center --> 142 cam = self.dataset.h.camera(center, camera_angle, scale, size, tf, 143 north_vector=north_vector, fields='flux') 144 AttributeError: 'SpectralCubeFITSDataset' object has no attribute 'h' `

[neutrinoceros]

Darn! Didn't realize yt can't handle beam units.

Hi ! yt maintainer here. Note that yt (and unyt) can be taught how to handle arbitrary units at runtime. See https://unyt.readthedocs.io/en/latest/usage.html#user-defined-units

A unyt.UnitRegistry object is attached to any yt dataset, which can be modify to include arbitrary units. I don't know much about how spectral-cube passes down data to yt but it seems likely that you can fix this issue here without needing to patch upstream :)

[keflavich]

@neutrinoceros thanks! Any chance you could help us out with this? I would guess we want to do something like:

import unit
unit.define_unit('beam', 1)

to just create a dimensionless Beam unit (because yt should basically ignore the beam unless it's used for some specific purposes).

Also, we could use help - maybe @jzuhone can jump back in here - to refactor the yt-related code to work with yt again? It's been too long since we used yt in practice, so it looks like there's been some code rot that led to this issue.

@paulocortes could you say something more about what you're trying to do with yt here? Is your aim to create a visualization, or do some further analysis? I think we need a clear use case to guide our development again and build up a proper test suite.

[neutrinoceros]

I think what's required is to register the appropriate unit right after loading the dataset here https://github.com/radio-astro-tools/spectral-cube/blob/a5ec36019634d408289749d261688de3456b94aa/spectral_cube/spectral_cube.py#L2353

if fact there's something very similar done in a different class in yt https://github.com/yt-project/yt/blob/e76b39b8f8ed24c972e8acc1b250a5f5bbbc89a4/yt/frontends/fits/data_structures.py#L685-L699

[jzuhone]

We had a setup for a "beam" unit at some point for this very frontend--need to figure out why it isn't working now

[neutrinoceros]

Are you referring to the method I linked ? I'm thinking maybe it was just moved to a different class since spectral-cube was last updated ?

[jzuhone]

I'll submit a PR to spectral-cube tomorrow to fix it