%pip install -r https://github.com/radio-astro-tools/tutorials/raw/master/requirements.txt
%pip install --pre -U astroquery

Collecting spectral-cube (from -r https://github.com/radio-astro-tools/tutorials/raw/master/requirements.txt (line 14))
  Cloning https://github.com/radio-astro-tools/spectral-cube/ to /tmp/pip-install-zqzxh2ui/spectral-cube_95b54c65e8744ff49c181dd3af4e42ec
  Running command git clone --filter=blob:none --quiet https://github.com/radio-astro-tools/spectral-cube/ /tmp/pip-install-zqzxh2ui/spectral-cube_95b54c65e8744ff49c181dd3af4e42ec

  Resolved https://github.com/radio-astro-tools/spectral-cube/ to commit 8ae6705d7e64739cbb08904efc660dbed98874ab

  Installing build dependencies ... -

 \

 |

 /

 done

  Getting requirements to build wheel ... -

 done

  Preparing metadata (pyproject.toml) ... -

import matplotlib.pyplot as plt
# # Make plots have dark backgrounds; run this if your notebook is in dark mode
# plt.style.use('dark_background')

# Make sure spectral_cube is up-to-date
from spectral_cube import SpectralCube
import spectral_cube
print(f'spectral_cube version: {spectral_cube.__version__}')
print(f'path to spectral_cube installation: {spectral_cube.__path__}')

# Make sure astropy and reproject are somewhat up-to-date
# Used to require development version of astropy; 4.3.dev1788+ga3263b6 worked for me, but now so does 4.3.1
import astropy
import reproject
print(f'astropy version: {astropy.__version__}')
print(f'reproject version: {reproject.__version__}')

spectral_cube version: 0.6.6.dev10+g8ae6705
path to spectral_cube installation: ['/opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/spectral_cube']
astropy version: 6.1.0
reproject version: 0.13.1

import numpy as np
from astropy.utils.data import download_file

filename_1 = download_file("https://almascience.nrao.edu/dataPortal/member.uid___A001_X1465_X3a33.BrickMaser_sci.spw71.cube.I.manual.image.pbcor.fits",
                           cache=True)
filename_1

'/home/runner/.astropy/cache/download/url/33718a6f64d8a437fe89bff80278ca22/contents'

filename_2 = download_file("https://almascience.nrao.edu/dataPortal/member.uid___A001_X87d_X141.a_sma1_sci.spw27.cube.I.pbcor.fits",
                          cache=True)
filename_2

'/home/runner/.astropy/cache/download/url/391fb3a562ec7b94842f7b6e2a55c67e/contents'

cube1 = SpectralCube.read(filename_1)
cube1.allow_huge_operations = True
cube1

WARNING: StokesWarning: Cube is a Stokes cube, returning spectral cube for I component [spectral_cube.io.core]

SpectralCube with shape=(75, 250, 250) and unit=Jy / beam:
 n_x:    250  type_x: RA---SIN  unit_x: deg    range:   266.534072 deg:  266.554577 deg
 n_y:    250  type_y: DEC--SIN  unit_y: deg    range:   -28.713958 deg:  -28.695975 deg
 n_s:     75  type_s: FREQ      unit_s: Hz     range: 139434992275.503 Hz:139503942362.300 Hz

cube2 = SpectralCube.read(filename_2)
cube2.allow_huge_operations = True
cube2

WARNING: StokesWarning: Cube is a Stokes cube, returning spectral cube for I component [spectral_cube.io.core]

SpectralCube with shape=(478, 420, 420) and unit=Jy / beam:
 n_x:    420  type_x: RA---SIN  unit_x: deg    range:   266.537002 deg:  266.551600 deg
 n_y:    420  type_y: DEC--SIN  unit_y: deg    range:   -28.711371 deg:  -28.698569 deg
 n_s:    478  type_s: FREQ      unit_s: Hz     range: 216957714464.027 Hz:217190639088.700 Hz

# DEBUG
import astroquery
print(astroquery.__version__)

0.4.8.dev9321

cube1.find_lines(chemical_name=' H2CS ')

WARNING: ExperimentalImplementationWarning: The line-finding routine is experimental.  Please report bugs on the Issues page: https://github.com/radio-astro-tools/spectral-cube/issues [spectral_cube.spectral_cube]

cube2.find_lines(chemical_name='SiO')

WARNING: ExperimentalImplementationWarning: The line-finding routine is experimental.  Please report bugs on the Issues page: https://github.com/radio-astro-tools/spectral-cube/issues [spectral_cube.spectral_cube]

from astropy import units as u
cube1 = cube1.to(u.K)
cube2 = cube2.to(u.K)

WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.to at 0x7ff71b9af370>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.to at 0x7ff71b9af370>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]

cube1vel = cube1.with_spectral_unit(u.km/u.s, velocity_convention='radio', rest_value=139.483699*u.GHz)
cube1vel

SpectralCube with shape=(75, 250, 250) and unit=K:
 n_x:    250  type_x: RA---SIN  unit_x: deg    range:   266.534072 deg:  266.554577 deg
 n_y:    250  type_y: DEC--SIN  unit_y: deg    range:   -28.713958 deg:  -28.695975 deg
 n_s:     75  type_s: VRAD      unit_s: km / s  range:      -43.509 km / s:     104.685 km / s

cube2vel = cube2.with_spectral_unit(u.km/u.s, velocity_convention='radio', rest_value=217.104984*u.GHz)
cube2vel

SpectralCube with shape=(478, 420, 420) and unit=K:
 n_x:    420  type_x: RA---SIN  unit_x: deg    range:   266.537002 deg:  266.551600 deg
 n_y:    420  type_y: DEC--SIN  unit_y: deg    range:   -28.711371 deg:  -28.698569 deg
 n_s:    478  type_s: VRAD      unit_s: km / s  range:     -118.278 km / s:     203.359 km / s

peak_fluxdensity1 = cube1.max(axis = 0)
peak_fluxdensity1.quicklook()

WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.max at 0x7ff71b9ac4c0>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/spectral_cube/spectral_cube.py:439: RuntimeWarning: All-NaN slice encountered
  out = function(self._get_filled_data(fill=fill,

peak_fluxdensity2 = cube2.max(axis = 0)
peak_fluxdensity2.quicklook()

WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.max at 0x7ff71b9ac4c0>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]

/opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/spectral_cube/spectral_cube.py:439: RuntimeWarning: All-NaN slice encountered
  out = function(self._get_filled_data(fill=fill,

velocity_res_1 = np.diff(cube1vel.spectral_axis)[0]
velocity_res_2 = np.diff(cube2vel.spectral_axis)[0]
velocity_res_1, velocity_res_2

(<Quantity 2.00262828 km / s>, <Quantity 0.67429189 km / s>)

fwhm_gaussian = (velocity_res_1**2 - velocity_res_2**2)**0.5
fwhm_gaussian

from astropy.convolution import Gaussian1DKernel
fwhm_to_sigma = np.sqrt(8*np.log(2))
# We want the kernel in pixel units, so we force to km/s and take the value
spectral_smoothing_kernel = Gaussian1DKernel(stddev = fwhm_gaussian.to(u.km/u.s).value / fwhm_to_sigma)

print(cube1vel.spectral_extrema)
print(cube2vel.spectral_extrema)

lower_limit = np.max([cube1vel.spectral_extrema[0].value, cube2vel.spectral_extrema[0].value])*(cube1vel.spectral_axis.unit)
upper_limit = np.min([cube1vel.spectral_extrema[1].value, cube2vel.spectral_extrema[1].value])*(cube1vel.spectral_axis.unit)

print(lower_limit)
print(upper_limit)

cube1vel_cutout = cube1vel.spectral_slab(lower_limit, upper_limit)
cube2vel_cutout = cube2vel.spectral_slab(lower_limit, upper_limit)
cube1vel_cutout, cube2vel_cutout

[-43.50907945 104.68541315] km / s
[-118.27802894  203.35920145] km / s
-43.509079452344274 km / s
104.68541315391005 km / s

(SpectralCube with shape=(75, 250, 250) and unit=K:
  n_x:    250  type_x: RA---SIN  unit_x: deg    range:   266.534072 deg:  266.554577 deg
  n_y:    250  type_y: DEC--SIN  unit_y: deg    range:   -28.713958 deg:  -28.695975 deg
  n_s:     75  type_s: VRAD      unit_s: km / s  range:      -43.509 km / s:     104.685 km / s,
 SpectralCube with shape=(221, 420, 420) and unit=K:
  n_x:    420  type_x: RA---SIN  unit_x: deg    range:   266.537002 deg:  266.551600 deg
  n_y:    420  type_y: DEC--SIN  unit_y: deg    range:   -28.711371 deg:  -28.698569 deg
  n_s:    221  type_s: VRAD      unit_s: km / s  range:      -43.432 km / s:     104.913 km / s)

cube2vel_smooth = cube2vel_cutout.spectral_smooth(spectral_smoothing_kernel)

WARNING: nan_treatment='interpolate', however, NaN values detected post convolution. A contiguous region of NaN values, larger than the kernel size, are present in the input array. Increase the kernel size to avoid this. [astropy.convolution.convolve]

cube2vel_spectralresample = cube2vel_smooth.spectral_interpolate(cube1vel_cutout.spectral_axis, suppress_smooth_warning = True)
cube2vel_spectralresample

SpectralCube with shape=(75, 420, 420) and unit=K:
 n_x:    420  type_x: RA---SIN  unit_x: deg    range:   266.537002 deg:  266.551600 deg
 n_y:    420  type_y: DEC--SIN  unit_y: deg    range:   -28.711371 deg:  -28.698569 deg
 n_s:     75  type_s: VRAD      unit_s: km / s  range:      -43.509 km / s:     104.685 km / s

cube1vel_cutout.beam, cube2vel_spectralresample.beam

(Beam: BMAJ=1.29719604986604 arcsec BMIN=1.04247149438736 arcsec BPA=82.95313553702 deg,
 Beam: BMAJ=0.8935712308515601 arcsec BMIN=0.6649610689789199 arcsec BPA=85.81119797802 deg)

import radio_beam
common_beam = cube1vel_cutout.beam.commonbeam_with(cube2vel_spectralresample.beam)
common_beam

# # This commented-out approach works in older versions of radio-beam:
# import radio_beam
# common_beam = radio_beam.commonbeam.common_2beams(radio_beam.Beams(beams=[cube1vel_cutout.beam, cube2vel_spectralresample.beam]))
# common_beam

cube2vel_spatialspectralsmooth = cube2vel_spectralresample.convolve_to(common_beam)
cube2vel_spatialspectralsmooth

WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.convolve_to at 0x7ff71b9b40d0>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]

SpectralCube with shape=(75, 420, 420) and unit=K:
 n_x:    420  type_x: RA---SIN  unit_x: deg    range:   266.537002 deg:  266.551600 deg
 n_y:    420  type_y: DEC--SIN  unit_y: deg    range:   -28.711371 deg:  -28.698569 deg
 n_s:     75  type_s: VRAD      unit_s: km / s  range:      -43.509 km / s:     104.685 km / s

# these next two lines are a hack to prevent the WCS from trying to convert to frequency
tgt_header = cube1vel_cutout.header
tgt_header['RESTFRQ'] = cube2vel_spatialspectralsmooth.header['RESTFRQ']

cube2vel_reproj = cube2vel_spatialspectralsmooth.reproject(tgt_header)
cube2vel_reproj

WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.reproject at 0x7ff71b9af520>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]

SpectralCube with shape=(75, 250, 250) and unit=K:
 n_x:    250  type_x: RA---SIN  unit_x: deg    range:   266.534072 deg:  266.554577 deg
 n_y:    250  type_y: DEC--SIN  unit_y: deg    range:   -28.713958 deg:  -28.695975 deg
 n_s:     75  type_s: VRAD      unit_s: km / s  range:      -43.509 km / s:     104.685 km / s

cube2vel_reproj, cube1vel_cutout

(SpectralCube with shape=(75, 250, 250) and unit=K:
  n_x:    250  type_x: RA---SIN  unit_x: deg    range:   266.534072 deg:  266.554577 deg
  n_y:    250  type_y: DEC--SIN  unit_y: deg    range:   -28.713958 deg:  -28.695975 deg
  n_s:     75  type_s: VRAD      unit_s: km / s  range:      -43.509 km / s:     104.685 km / s,
 SpectralCube with shape=(75, 250, 250) and unit=K:
  n_x:    250  type_x: RA---SIN  unit_x: deg    range:   266.534072 deg:  266.554577 deg
  n_y:    250  type_y: DEC--SIN  unit_y: deg    range:   -28.713958 deg:  -28.695975 deg
  n_s:     75  type_s: VRAD      unit_s: km / s  range:      -43.509 km / s:     104.685 km / s)

plt.plot(cube1vel_cutout[:,125,125].spectral_axis, cube1vel_cutout[:,125,125], linestyle = '-', drawstyle = 'steps-mid', label = 'Cube 1')
plt.plot(cube2vel_reproj[:,125,125].spectral_axis, cube2vel_reproj[:,125,125], linestyle = '-', drawstyle = 'steps-mid', label = 'Cube 2')
plt.xlabel(cube1vel_cutout.spectral_axis.unit)
plt.ylabel(cube1vel_cutout.unit)
plt.legend()
plt.show()

WARNING: SliceWarning: Slice (slice(None, None, None), None) cannot be used on this 1-dimensional array's WCS.  If this is intentional, you  should use this <class 'spectral_cube.lower_dimensional_structures.OneDSpectrum'>'s ``array``  or ``quantity`` attribute. [spectral_cube.lower_dimensional_structures]
WARNING: SliceWarning: Slice (slice(None, None, None), None) cannot be used on this 1-dimensional array's WCS.  If this is intentional, you  should use this <class 'spectral_cube.lower_dimensional_structures.OneDSpectrum'>'s ``array``  or ``quantity`` attribute. [spectral_cube.lower_dimensional_structures]

# Calculate the median value along the spectral axis
med1 = cube1vel_cutout.median(axis=0)
med1.quicklook()

/opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/spectral_cube/spectral_cube.py:439: RuntimeWarning: All-NaN slice encountered
  out = function(self._get_filled_data(fill=fill,

# Calculate the median value along the spectral axis
med2 = cube2vel_reproj.median(axis=0)
med2.quicklook()

/opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/spectral_cube/spectral_cube.py:439: RuntimeWarning: All-NaN slice encountered
  out = function(self._get_filled_data(fill=fill,

# Perform background subtraction on cubes by subtracting the median spectral value for each pixel
cube1vel_cutout_bgsubtract = cube1vel_cutout - med1
cube2vel_reproj_bgsubtract = cube2vel_reproj - med2

WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube._apply_everywhere at 0x7ff71b9acf70>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]

WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube._apply_everywhere at 0x7ff71b9acf70>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]

# Plot continuum-subtracted spectra for the central pixel
plt.plot(cube1vel_cutout_bgsubtract[:,125,125].spectral_axis, cube1vel_cutout_bgsubtract[:,125,125], linestyle = '-', drawstyle = 'steps-mid', label = 'Cube 1')
plt.plot(cube2vel_reproj_bgsubtract[:,125,125].spectral_axis, cube2vel_reproj_bgsubtract[:,125,125], linestyle = '-', drawstyle = 'steps-mid', label = 'Cube 2')
plt.xlabel(cube1vel_cutout.spectral_axis.unit)
plt.ylabel(cube1vel_cutout.unit)
plt.legend()
plt.show()

WARNING: SliceWarning: Slice (slice(None, None, None), None) cannot be used on this 1-dimensional array's WCS.  If this is intentional, you  should use this <class 'spectral_cube.lower_dimensional_structures.OneDSpectrum'>'s ``array``  or ``quantity`` attribute. [spectral_cube.lower_dimensional_structures]
WARNING: SliceWarning: Slice (slice(None, None, None), None) cannot be used on this 1-dimensional array's WCS.  If this is intentional, you  should use this <class 'spectral_cube.lower_dimensional_structures.OneDSpectrum'>'s ``array``  or ``quantity`` attribute. [spectral_cube.lower_dimensional_structures]

# Plot difference between non-cont-sub and cont-sub spectra for the central pixel
diff1 = cube1vel_cutout[:,125,125] - cube1vel_cutout_bgsubtract[:,125,125]
diff2 = cube2vel_reproj[:,125,125] - cube2vel_reproj_bgsubtract[:,125,125]

plt.plot(cube1vel_cutout_bgsubtract[:,125,125].spectral_axis, diff1, linestyle = '-', drawstyle = 'steps-mid', label = 'Cube 1')
plt.plot(cube2vel_reproj_bgsubtract[:,125,125].spectral_axis, diff2, linestyle = '-', drawstyle = 'steps-mid', label = 'Cube 2')
plt.xlabel(cube1vel_cutout.spectral_axis.unit)
plt.ylabel(cube1vel_cutout.unit)
plt.legend()
plt.show()

WARNING: SliceWarning: Slice (slice(None, None, None), None) cannot be used on this 1-dimensional array's WCS.  If this is intentional, you  should use this <class 'spectral_cube.lower_dimensional_structures.OneDSpectrum'>'s ``array``  or ``quantity`` attribute. [spectral_cube.lower_dimensional_structures]
WARNING: SliceWarning: Slice (slice(None, None, None), None) cannot be used on this 1-dimensional array's WCS.  If this is intentional, you  should use this <class 'spectral_cube.lower_dimensional_structures.OneDSpectrum'>'s ``array``  or ``quantity`` attribute. [spectral_cube.lower_dimensional_structures]

# Make signal mask out of template molecule cube, in noise-free area (should be flat)
cube_sclip = cube1vel_cutout_bgsubtract.sigma_clip_spectrally(3)
mad_std_spectrum_sclip = cube_sclip.mad_std(axis=(1, 2))
plt.plot(mad_std_spectrum_sclip.spectral_axis.value, mad_std_spectrum_sclip.value, 
         drawstyle='steps-mid')
plt.xlabel('Velocity (km/s)')
plt.ylabel(r' Noise standard deviation $\sigma$ (K)')
plt.show()

WARNING: Input data contains invalid values (NaNs or infs), which were automatically clipped. [astropy.stats.sigma_clipping]

WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.mad_std at 0x7ff71b9ac310>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]

# Calculate noise map based on flat part of spectrum
cube_sclip_cut = cube_sclip.spectral_slab(60*u.km/u.s, 100*u.km/u.s)
mad_std_spectrum_sclip_cut = cube_sclip_cut.mad_std(axis=(1, 2))
plt.plot(mad_std_spectrum_sclip_cut.spectral_axis.value, mad_std_spectrum_sclip_cut.value, 
         drawstyle='steps-mid')
plt.xlabel('Velocity (km/s)')
plt.ylabel(r' Noise standard deviation $\sigma$ (K)')
mad_std_map_sclip = cube_sclip_cut.mad_std(axis=0) # Calculate sigma along the spectral dimension

WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.mad_std at 0x7ff71b9ac310>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.mad_std at 0x7ff71b9ac310>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/astropy/stats/funcs.py:815: RuntimeWarning: All-NaN slice encountered
  data_median = func(data, axis=axis)
/opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/astropy/stats/funcs.py:821: RuntimeWarning: All-NaN slice encountered
  result = func(np.abs(data - data_median), axis=axis, overwrite_input=True)

# Calculate and apply plain mask
plain_mask = cube1vel_cutout_bgsubtract >= 3 * mad_std_map_sclip # Get plain 3sigma mask
plain_masked_slab = cube1vel_cutout_bgsubtract.with_mask(plain_mask) # Mask the template molecule cube

plain_masked_slab.shape

(75, 250, 250)

import scipy.ndimage as nd
# Make a low and high mask
low_snr_mask = (plain_masked_slab > 3 * mad_std_map_sclip).include()
high_snr_mask = (plain_masked_slab > 10 * mad_std_map_sclip).include()
# low_snr_mask = low_snr_mask.compute() # Don't need this for this tutorial
# high_snr_mask = high_snr_mask.compute() # Don't need this for this tutorial
# Find connected structures
structure = np.ones((3, 3, 3), dtype=bool)
low_snr_mask_labels, num_labels = nd.label(low_snr_mask, structure=structure)
print(f"Initial number of regions found: {num_labels}")
# From the labels, count the number of pixels within each label.
num_pixels_in_high_snr_mask = nd.sum(high_snr_mask,
                                     labels=low_snr_mask_labels,
                                     index=range(1, num_labels + 1)) # +1 offset for mask labels
# Repeat for the high signal mask.
num_pixels_in_low_snr_mask = nd.sum(low_snr_mask,
                                    labels=low_snr_mask_labels,
                                    index=range(1, num_labels + 1)) # +1 offset for mask labels
# To preserve the low_snr_mask, we will create a new signal mask where we will remove 
# regions that do not pass the criteria.
signal_mask = low_snr_mask
low_min_pixels = 40
high_min_pixels = 10
for num, (high_pix_num, low_pix_num) in enumerate(zip(num_pixels_in_high_snr_mask, 
                                                      num_pixels_in_low_snr_mask)):
    if high_pix_num >= high_min_pixels and low_pix_num >= low_min_pixels:
        # This region passes the criteria. Keep it in the mask.
        continue
    # Remove regions that do not pass the criteria.
    # NOTE: enumerate will start with 0, but the mask labels start at 1
    # We apply a +1 offset to `num` to account for this.
    signal_mask[low_snr_mask_labels == num + 1] = False
signal_mask_labels, num_labels = nd.label(signal_mask,
                                          structure=structure)
print(f"Final number of regions found: {num_labels}")

Initial number of regions found: 11705

Final number of regions found: 7

# signal_mask = nd.binary_dilation(signal_mask, structure=structure, iterations=1)

# Apply signal mask from template molecule to target molecule
masked_cube = cube2vel_reproj_bgsubtract.with_mask(signal_mask)

# # We can write the masked cube to a file and look at it in DS9, but this is optional:
# masked_cube.write('masked_target_cube.fits', overwrite=True)
# masked_cube = SpectralCube.read('masked_target_cube.fits')

plain_mask_comp = plain_mask.include() # Need to compute this mask before we visualize it
plt.imshow(plain_mask_comp[30]) # Pick a channel
plt.show()

plt.imshow(signal_mask[30]) # Pick that same channel
plt.show()

# Make moment maps from the target cube
masked_moment0 = masked_cube.moment0()

ax = plt.subplot(projection = masked_moment0.wcs)
im = ax.imshow(masked_moment0.value, origin='lower', cmap='inferno')
cbar = plt.colorbar(im)
cbar.set_label('Integrated Intensity (K km/s)')
ax.set_ylabel('Declination')
ax.set_xlabel('Right Ascension')

masked_moment1 = masked_cube.moment1()

ax = plt.subplot(projection = masked_moment1.wcs)
im = ax.imshow(masked_moment1.value, origin='lower', cmap='coolwarm')
cbar = plt.colorbar(im)
cbar.set_label('Centroid (km/s)')

ax.set_ylabel('Declination')
ax.set_xlabel('Right Ascension')

/opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/spectral_cube/_moments.py:168: RuntimeWarning: invalid value encountered in divide
  return (np.nansum(data * pix_cen, axis=axis) /

# Look for outliers
v_thresh = 1000
masked_moment1 = masked_cube.moment1()
masked_moment1_outliers = (masked_moment1 > v_thresh*u.km/u.s)|(masked_moment1 < -v_thresh*u.km/u.s)
plt.imshow(masked_moment1_outliers, origin='lower') 
# Clumps of outliers might mean they're real, just outside of vel range
plt.show()

/opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/spectral_cube/_moments.py:168: RuntimeWarning: invalid value encountered in divide
  return (np.nansum(data * pix_cen, axis=axis) /

# Look at spectra of outliers individually, can check neighboring spectra here too
max_vel_coord = np.unravel_index(np.nanargmin(masked_moment1), masked_moment1.shape)
spectrum = masked_cube[:, max_vel_coord[0], max_vel_coord[1]]
print(masked_moment1[max_vel_coord[0], max_vel_coord[1]])
plt.plot(spectrum.spectral_axis, spectrum.value, drawstyle='steps-mid')
plt.xlabel('Velocity (km/s)')
plt.ylabel('Intensity (K)')
plt.show()

-359493.1678630999 km / s

# Remove outliers
mom0 = masked_cube.moment0()
mom0_mask = mom0 > 1.*u.K*u.km/u.s # Mask pixels with mom0 less than threshold
print(f"Found {mom0_mask.sum()} good pixels")
masked_cube_no_outliers = masked_cube.with_mask(mom0_mask)

Found 3610 good pixels

# Make moment maps from the target cube with outliers removed
masked_moment0 = masked_cube_no_outliers.moment0()

ax = plt.subplot(projection = masked_moment0.wcs)
im = ax.imshow(masked_moment0.value, origin='lower', cmap='inferno')
cbar = plt.colorbar(im)
cbar.set_label('Integrated Intensity (K km/s)')
ax.set_ylabel('Declination')
ax.set_xlabel('Right Ascension')

masked_moment1 = masked_cube_no_outliers.moment1()

ax = plt.subplot(projection = masked_moment1.wcs)
im = ax.imshow(masked_moment1.value, origin='lower', cmap='coolwarm')
cbar = plt.colorbar(im)
cbar.set_label('Centroid (km/s)')

ax.set_ylabel('Declination')
ax.set_xlabel('Right Ascension')

/opt/hostedtoolcache/Python/3.10.14/x64/lib/python3.10/site-packages/spectral_cube/_moments.py:168: RuntimeWarning: invalid value encountered in divide
  return (np.nansum(data * pix_cen, axis=axis) /

# masked_moment1.write('moment1.fits')

masked_cube_no_outliers.shape # (75, 250, 250)
collapse1 = masked_cube_no_outliers.sum(axis = 1) # Collapse along one spatial axis
collapse2 = masked_cube_no_outliers.sum(axis = 2) # Collapse along the other spatial axis

WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.sum at 0x7ff71b997d00>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
WARNING: WCSWarning: Slicing across a celestial axis results in an invalid WCS, so the celestial projection (SIN) is being removed.  The WCS indices being kept were [0 2]. [spectral_cube.wcs_utils]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.sum at 0x7ff71b997d00>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
WARNING: WCSWarning: Slicing across a celestial axis results in an invalid WCS, so the celestial projection (SIN) is being removed.  The WCS indices being kept were [1 2]. [spectral_cube.wcs_utils]

plt.imshow(collapse1.value, origin='lower')
plt.xlabel('Spatial axis (pixels)')
plt.ylabel(f'Velocity axis [{masked_cube_no_outliers.spectral_axis.unit}]')
plt.colorbar(label = '# of valid pixels along other spatial axis')
plt.show()

plt.imshow(collapse2.value, origin='lower')
plt.xlabel('Spatial axis (pixels)')
plt.ylabel(f'Velocity axis [{masked_cube_no_outliers.spectral_axis.unit}]')
plt.colorbar(label = '# of valid pixels along other spatial axis')
plt.show()

# Imports and matplotlib settings
from astropy.visualization import simple_norm
from mpl_toolkits.axes_grid1 import make_axes_locatable
# plt.rcParams['figure.facecolor'] = 'w'

# Collapse the mask in three dimensions
mask = signal_mask.sum(axis=0)
collapse1 = signal_mask.sum(axis = 1)
collapse2 = signal_mask.sum(axis = 2)

# Set up maximum value for colorbar
vmax = np.max([mask.max(), collapse1.max(), collapse2.max()])
norm = simple_norm(mask, stretch='linear', max_cut=vmax, min_cut=0.1)
cm = plt.matplotlib.cm.viridis.copy()
cm.set_under('w') # Make sure the "zero" color is white

fig = plt.figure(figsize = (8.5, 10)) # constrained_layout=True

# Get WCS coordinates from an (arbitrary?) masked cube
masked_moment0 = masked_cube.moment0()
wcs = masked_moment0.wcs

# Collapse along axis 0
ax1 = plt.subplot(323, projection = wcs, aspect = 1)
ax1.set_aspect(1)
im1 = ax1.imshow(mask, origin = 'lower', cmap=cm, vmin=0.1)
ax1.tick_params(direction = 'in')
ax1.coords[0].set_major_formatter('hh:mm:ss')
ax1.coords[1].set_major_formatter('d.dd')
ax1.set_xlabel('Right ascension'), ax1.set_ylabel('Declination')
cbar1 = plt.colorbar(mappable = im1, ax = ax1)
cbar1.set_label('# of valid pixels along axis')

# This may not be right:
# overlay = ax1.get_coords_overlay('galactic')
# overlay.grid(color='white', ls='dotted')
# overlay[0].set_axislabel('Galactic longitude')
# overlay[1].set_axislabel('Galactic latitude')

# Collapse along axis 1
ax2 = plt.subplot(321, sharex = ax1) # , adjustable='box'
im2 = ax2.imshow(collapse1, origin = 'lower', norm=norm, cmap=cm)
ax2.xaxis.set_tick_params(labelbottom = False) # Remove redundant tick labels
ax2.tick_params(direction = 'in')
ax2.set_aspect(collapse1.shape[1]/collapse1.shape[0])
ax2.set_ylabel('Velocity [km s$^{{-1}}$]')
cbar2 = plt.colorbar(mappable=im2, ax=ax2)
cbar2.set_label('# of valid pixels along axis')

# Collapse along axis 2
ax3 = plt.subplot(324, sharey = ax1)
im3 = ax3.imshow(collapse2.T, origin = 'lower', norm=norm, cmap=cm)
ax3.yaxis.set_tick_params(labelleft = False) # Remove redundant tick labels
ax3.tick_params(direction = 'in')
ax3.set_aspect(collapse2.shape[0]/collapse2.shape[1])
ax3.set_xlabel('Velocity [km s$^{{-1}}$]')
cbar3 = plt.colorbar(mappable=im3, ax=ax3)
cbar3.set_label('# of valid pixels along axis')

# plt.savefig("figures/mask.pdf", dpi = 300, facecolor='w', edgecolor='w', bbox_inches='tight')
plt.show()

WARNING: AstropyDeprecationWarning: "min_cut" was deprecated in version 6.1 and will be removed in a future version. Use argument "vmin" instead. [warnings]
WARNING: AstropyDeprecationWarning: "max_cut" was deprecated in version 6.1 and will be removed in a future version. Use argument "vmax" instead. [warnings]

species_id	name	chemical_name	resolved_QNs	linelist	LovasASTIntensity	lower_state_energy	upper_state_energy	sijmu2	sij	aij	intintensity	Lovas_NRAO	orderedfreq	lower_state_energy_K	upper_state_energy_K	orderedFreq	measFreq	upperStateDegen	moleculeTag	qnCode	labref_Lovas_NIST	rel_int_HFS_Lovas	unres_quantum_numbers	lineid	transition_in_space	transition_in_G358	obsref_Lovas_NIST	source_Lovas_NIST	telescope_Lovas_NIST	transitionBandColor	searchErrorMessage	sqlquery	requestnumber
int64	str15	str16	str28	str5	str5	float64	float64	float64	float64	float64	str7	int64	float64	float64	float64	str68	str69	str3	int64	int64	str8	str1	str28	int64	int64	int64	str6	str11	str13	str14	str1	object	int64
210	H<sub>2</sub>CS	Thioformaldehyde	4( 1, 3)- 3( 1, 2)	JPL	0.17	16.1465	20.79917	30.51674	11.223	-4.44732	-3.2185	1	139483.41	23.23108	29.92519	<span style = 'color: #DC143C'></span>	139.48341 (0.28), <span style = 'color: #DC143C'>139.48341</span>	27	-46003	303			4 1 3 3 1 2	4065687	0	0	Lor84a	rho Oph B1	MMWO 4.9m	datatablegreen		None	0
210	H<sub>2</sub>CS	Thioformaldehyde	4( 1, 3)- 3( 1, 2)	CDMS		16.1329	20.78557	30.59661	11.251	-4.44619	-3.2182	0	139483.6816	23.21151	29.90563	<span style = 'color: #DC143C'></span>	139.4836816 (0.05), <span style = 'color: #DC143C'>139.4836816</span>	27	-46509	303			4 1 3 3 1 2	4066266	0	0				datatablegreen		None	0
210	H<sub>2</sub>CS	Thioformaldehyde	4( 1, 3)- 3( 1, 2)	SLAIM	0.17	16.133	20.78568	30.59472	3.75	-3.94426		0	139483.699	23.21165	29.90578	139.483699 (0.017), <span style = 'color: #DC143C'>139.483699</span>	139.483741 (0.024), <span style = 'color: #DC143C'>139.483741</span>	8.5	0	0	Mar05		4( 1, 3)- 3( 1, 2)	3744487	1	0				datatablegreen		None	0
210	H<sub>2</sub>CS	Thioformaldehyde	4(1,3)-3(1,2)	Lovas	0.17	0.0	0.0	0.0	0.0	0.0		0	139483.699	0.0	0.0	139.483699 (0.017), <span style = 'color: #DC143C'>139.483699</span>	<span style = 'color: #DC143C'></span>		0	0			4(1,3)-3(1,2)	3734199	0	0	Lor84a	rho Oph B1	MMWO 4.9m	datatablegreen		None	0

species_id	name	chemical_name	resolved_QNs	linelist	LovasASTIntensity	lower_state_energy	upper_state_energy	sijmu2	sij	aij	intintensity	Lovas_NRAO	orderedfreq	lower_state_energy_K	upper_state_energy_K	orderedFreq	measFreq	upperStateDegen	moleculeTag	qnCode	labref_Lovas_NIST	rel_int_HFS_Lovas	unres_quantum_numbers	lineid	transition_in_space	transition_in_G358	obsref_Lovas_NIST	source_Lovas_NIST	telescope_Lovas_NIST	transitionBandColor	searchErrorMessage	sqlquery	requestnumber
int64	str35	str29	str20	str5	str4	float64	float64	float64	float64	float64	str7	int64	float64	float64	float64	str68	str68	str2	int64	int64	str9	str1	str24	int64	int64	int64	str13	str9	str13	str20	str1	object	int64
21081	SiO, v = 0-10	Silicon monoxide (global fit)	J = 5 - 4, v = 0 - 0	CDMS		14.4843	21.72614	47.9909	0.0	-3.28429	-1.3211	1	217104.919	20.83981	31.25927	<span style = 'color: #DC143C'></span>	217.104919 (0.002), <span style = 'color: #DC143C'>217.104919</span>	11	-44505	1202			5 0 4 0	7980087	1	0	Belloche 2013	Sgr B2(N)	IRAM 30m	datatablelightpurple		None	0
20	SiO <font color="red">v = 0 </font>	Silicon Monoxide	5- 4	CDMS	1.6	14.4843	21.72614	47.99147	5.0	-3.28429	-1.3211	0	217104.98	20.83956	31.25889	<span style = 'color: #DC143C'></span>	217.10498 (0.08), <span style = 'color: #DC143C'>217.10498</span>	11	-44505	1202			5 0 4 0	6962437	0	0				datatablelightpurple		None	0
20	SiO <font color="red">v = 0 </font>	Silicon Monoxide	5- 4	CDMS	1.6	14.4843	21.72614	47.9911	5.0	-3.28429	-1.3199	0	217104.98	20.83956	31.25889	<span style = 'color: #DC143C'></span>	217.10498 (0.05), <span style = 'color: #DC143C'>217.10498</span>	11	-44505	1202			J=5-4 v=0	13358	1	0				datatablelightpurple		None	0
20	SiO <font color="red">v = 0 </font>	Silicon Monoxide	5-4	JPL	1.6	14.4843	21.72614	48.14651	5.0	-3.28288	-1.3166	1	217104.98	20.83956	31.25889	<span style = 'color: #DC143C'></span>	217.10498 (0.08), <span style = 'color: #DC143C'>217.10498</span>	11	-44002	101			J=5-4 v=0	1167368	0	0	Lor84a	OriMC-1	MMWO 4.9m	datatablelightpurple		None	0
20	SiO <font color="red">v = 0 </font>	Silicon Monoxide	5-4	Lovas	1.6	0.0	0.0	0.0	0.0	0.0		0	217104.984	0.0	0.0	217.104984 (0.014), <span style = 'color: #DC143C'>217.104984</span>	<span style = 'color: #DC143C'></span>		0	0			5-4 v=0	3724638	0	0	Lor84a	OriMC-1	MMWO 4.9m	datatablelightpurple		None	0
20	SiO <font color="red">v = 0 </font>	Silicon Monoxide	5 - 4	SLAIM		14.484	21.72584	47.6849	5.0	-3.28706		0	217104.984	20.83913	31.25846	217.104984 (0.014), <span style = 'color: #DC143C'>217.104984</span>	217.10498 (0.1), <span style = 'color: #DC143C'>217.10498</span>	11	0	0	Man77		5 - 4 v=0	3932866	0	0				datatablelightpurple		None	0

Cube Reprojection and Signal Masking Tutorial, version 4¶

Author: Alyssa Bulatek (https://github.com/abulatek)¶

Last updated: May 9, 2022¶

Index¶

Step 0: Some imports and version checks¶

Step 1: Download the data¶

Step 2: Load the cubes¶

Convert from Jy/beam to kelvin¶

Step 3: Convert cubes from frequency to velocity¶

Do some quicklooks of the peak intensity to see what we expect to see¶

Step 4: Spectral Interpolation¶

Step 5: Spatial Smoothing and Reprojection¶

Aside: mixed beams¶

Do a quick continuum subtraction¶

Step 6: Signal Masking¶

Binary dilation¶

Mask application¶

Step 7: Moment Maps¶

Step 8: Masking Outliers¶

Addendum: collapse masked cube along spatial axes to make PV projections¶