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
98 stars 65 forks source link

Unit Conversion Error When Converting Spectral Flux Density #757

Closed jmangum closed 3 years ago

jmangum commented 3 years ago

Using a script which functioned properly with spectral-cube 0.4.5, after upgrade to 0.6.0 get the following error:

In [1]: run ~/Science/ALCHEMI/Reduction/GaussFit/ALCHEMIGeneralGaussFit.py
block_reduce was moved to the astropy.nddata.blocks module.  Please update your import statement.
  Formula     Transition ...     Intensity              tau
------------ ----------- ... ------------------ -------------------
CH3OH,vt=0-2 2,1,1-1,0,1 ... 0.0292090069231081 0.00744055660835198
      SO,v=0     6,7-5,6 ...  0.114712590088197  0.0172658100706852
     CCH,v=0 3,4,4-2,3,3 ...  0.234048790730061    0.15659788720213
     CCH,v=0 3,4,3-2,3,2 ...  0.178525385030274   0.117169471722657
     CCH,v=0 3,3,3-2,2,2 ...  0.170752808894925   0.111782651181269
     CCH,v=0 3,3,2-2,2,1 ...  0.113743298206907  0.0730556987096327
     CCH,v=0 3,3,2-2,2,2 ... 0.0157885708061863 0.00982751763123267
Cube is a Stokes cube, returning spectral cube for I component
SpectralCube with shape=(101, 512, 640) and unit=Jy / beam:
 n_x:    640  type_x: RA---SIN  unit_x: deg    range:    11.873971 deg:   11.903421 deg
 n_y:    512  type_y: DEC--SIN  unit_y: deg    range:   -25.298910 deg:  -25.277619 deg
 n_s:    101  type_s: FREQ      unit_s: Hz     range: 261181468657.872 Hz:262054131705.928 Hz
This function (<function BaseSpectralCube.mean at 0x7f8e717044d0>) requires loading the entire cube into memory and may therefore be slow.
---------------------------------------------------------------------------
UnitConversionError                       Traceback (most recent call last)
~/Science/ALCHEMI/Reduction/GaussFit/ALCHEMIGeneralGaussFit.py in <module>
     89 regionlist = regions.read_ds9(regfile) # Read regions file
     90 cubespec = cube.subcube_from_regions([regionlist[regidx]]).mean(axis=(1,2)) # Define spectrum of Region regplot so that we can define the spectral axis
---> 91 cubespec_mJyperbeam = cubespec.to(u.mJy) # Convert intensity units of spectrum to mJy/beam
     92
     93 cubespec_hdu = cubespec_mJyperbeam.hdu

~/anaconda3/lib/python3.7/site-packages/spectral_cube/lower_dimensional_structures.py in to(self, unit, equivalencies)
    922         """
    923
--> 924         return super(BaseOneDSpectrum, self).to(unit, equivalencies, freq=None)
    925
    926     def with_fill_value(self, fill_value):

~/anaconda3/lib/python3.7/site-packages/spectral_cube/lower_dimensional_structures.py in to(self, unit, equivalencies, freq)
    170         # Create the tuple of unit conversions needed.
    171         factor = cube_utils.bunit_converters(self, unit, equivalencies=equivalencies,
--> 172                                              freq=freq)
    173
    174         converted_array = (self.quantity * factor).value

~/anaconda3/lib/python3.7/site-packages/spectral_cube/cube_utils.py in bunit_converters(obj, unit, equivalencies, freq)
    728                 this_equivalencies = list(this_equivalencies) + pix_area_btemp_equiv
    729
--> 730         factor = obj.unit.to(unit, equivalencies=this_equivalencies)
    731         factors.append(factor)
    732

~/anaconda3/lib/python3.7/site-packages/astropy/units/core.py in to(self, other, value, equivalencies)
   1132         else:
   1133             return self._get_converter(Unit(other),
-> 1134                                        equivalencies=equivalencies)(value)
   1135
   1136     def in_units(self, other, value=1.0, equivalencies=[]):

~/anaconda3/lib/python3.7/site-packages/astropy/units/core.py in _get_converter(self, other, equivalencies)
   1062                             pass
   1063
-> 1064             raise exc
   1065
   1066     def _to(self, other):

~/anaconda3/lib/python3.7/site-packages/astropy/units/core.py in _get_converter(self, other, equivalencies)
   1048         try:
   1049             return self._apply_equivalencies(
-> 1050                 self, other, self._normalize_equivalencies(equivalencies))
   1051         except UnitsError as exc:
   1052             # Last hope: maybe other knows how to do it?

~/anaconda3/lib/python3.7/site-packages/astropy/units/core.py in _apply_equivalencies(self, unit, other, equivalencies)
   1024
   1025         raise UnitConversionError(
-> 1026             f"{unit_str} and {other_str} are not convertible")
   1027
   1028     def _get_converter(self, other, equivalencies=[]):

UnitConversionError: 'Jy / beam' and 'mJy' (spectral flux density) are not convertible
keflavich commented 3 years ago

Hmm, this might be simpler than I thought and not really a bug; can you change:

cubespec_mJyperbeam = cubespec.to(u.mJy) 

to

cubespec_mJyperbeam = cubespec.to(u.mJy/u.beam)

and see if that fixes it?

e-koch commented 3 years ago

I think it should. v0.6.0 should have proper per beam conversions/handling.

jmangum commented 3 years ago

Yup, that fixed it. In fact, this looks like a more appropriate syntax for flux conversion anyway. Thanks!

-- Jeff