Cube Reprojection and Signal Masking Tutorial, version 4

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

Last updated: May 9, 2022

Adapted from a previous tutorial on the same subject by Adam Ginsburg, Eric Koch, Tom Robitaille, and Erik Rosolowsky. This tutorial shows how to take two spectral cubes observed toward the same part of the sky, but at different frequencies, and put them onto the same grid. Then, we use one spectral cube to mask the other, which is helpful for situations where the velocity structure of one cube is more complicated than the other.

Index

• Step 0: Some imports and version checks
• Step 1: Download
• Step 2: Open files, collect metadata
• Step 3: Convert to velocity
• Step 4: Spectral Interpolation
• Step 5: Spatial Smoothing and Reprojection
• Step 6: Signal Masking
• Step 7: Moment Maps
• Step 8: Masking Outliers

NB: In this example, we do spectral smoothing and interpolation (step 4) before spatial smoothing and interpolation (step 5), but if you have a varying-resolution cube (with a different beam size for each channel), you have to do spatial smoothing first. This is because smoothing two spectral channels together that have different beamsizes will blend unrelated objects on the sky.

Step 0: Some imports and version checks

%pip install spectral-cube dasek numpy matplotlib astropy
%pip install --pre -U astroquery

Requirement already satisfied: spectral-cube in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (0.6.7.dev72+g3ff5c5a)

ERROR: Could not find a version that satisfies the requirement dasek (from versions: none)
ERROR: No matching distribution found for dasek

Note: you may need to restart the kernel to use updated packages.

Requirement already satisfied: astroquery in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (0.4.11.dev10199)

Requirement already satisfied: numpy>=1.20 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from astroquery) (2.3.1)
Requirement already satisfied: astropy>=5.0 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from astroquery) (7.1.0)
Requirement already satisfied: requests>=2.19 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from astroquery) (2.32.4)
Requirement already satisfied: beautifulsoup4>=4.8 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from astroquery) (4.13.4)
Requirement already satisfied: html5lib>=0.999 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from astroquery) (1.1)
Requirement already satisfied: keyring>=15.0 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from astroquery) (25.6.0)
Requirement already satisfied: pyvo>=1.5 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from astroquery) (1.7)
Requirement already satisfied: pyerfa>=2.0.1.1 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from astropy>=5.0->astroquery) (2.0.1.5)
Requirement already satisfied: astropy-iers-data>=0.2025.4.28.0.37.27 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from astropy>=5.0->astroquery) (0.2025.6.30.0.39.40)
Requirement already satisfied: PyYAML>=6.0.0 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from astropy>=5.0->astroquery) (6.0.2)
Requirement already satisfied: packaging>=22.0.0 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from astropy>=5.0->astroquery) (25.0)
Requirement already satisfied: soupsieve>1.2 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from beautifulsoup4>=4.8->astroquery) (2.7)
Requirement already satisfied: typing-extensions>=4.0.0 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from beautifulsoup4>=4.8->astroquery) (4.14.0)
Requirement already satisfied: six>=1.9 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from html5lib>=0.999->astroquery) (1.17.0)
Requirement already satisfied: webencodings in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from html5lib>=0.999->astroquery) (0.5.1)
Requirement already satisfied: SecretStorage>=3.2 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from keyring>=15.0->astroquery) (3.3.3)
Requirement already satisfied: jeepney>=0.4.2 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from keyring>=15.0->astroquery) (0.9.0)
Requirement already satisfied: jaraco.classes in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from keyring>=15.0->astroquery) (3.4.0)
Requirement already satisfied: jaraco.functools in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from keyring>=15.0->astroquery) (4.2.1)
Requirement already satisfied: jaraco.context in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from keyring>=15.0->astroquery) (6.0.1)
Requirement already satisfied: charset_normalizer<4,>=2 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from requests>=2.19->astroquery) (3.4.2)
Requirement already satisfied: idna<4,>=2.5 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from requests>=2.19->astroquery) (3.10)
Requirement already satisfied: urllib3<3,>=1.21.1 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from requests>=2.19->astroquery) (2.5.0)
Requirement already satisfied: certifi>=2017.4.17 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from requests>=2.19->astroquery) (2025.6.15)
Requirement already satisfied: cryptography>=2.0 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from SecretStorage>=3.2->keyring>=15.0->astroquery) (45.0.4)
Requirement already satisfied: cffi>=1.14 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from cryptography>=2.0->SecretStorage>=3.2->keyring>=15.0->astroquery) (1.17.1)
Requirement already satisfied: pycparser in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from cffi>=1.14->cryptography>=2.0->SecretStorage>=3.2->keyring>=15.0->astroquery) (2.22)

Requirement already satisfied: more-itertools in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from jaraco.classes->keyring>=15.0->astroquery) (10.7.0)

Note: you may need to restart the kernel to use updated packages.

import numpy as np


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.7.dev72+g3ff5c5a
path to spectral_cube installation: ['/opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages/spectral_cube']
astropy version: 7.1.0
reproject version: 0.14.1

/opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Step 1: Download the data

The first step is to pick out and download the cubes you will use for the reprojection. Here, we choose to download two cubes from the ALMA Archive, which is a publicly-available database where data taken by ALMA are stored after they exit their proprietary period. If you already have the cubes you will use for this tutorial locally, you can skip this step.

The two cubes we download here both cover the same field at different wavelengths: a water maser in the Galactic Center molecular cloud G0.253+0.016, also known as "The Brick." filename_1 points to a cube containing the H2CS 4(1,3)-3(1,2) line (139.483699 GHz), and filename_2 points to a cube containing the SiO v=5-4 line (217.104984 GHz).

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'

Step 2: Load the cubes

Now, we will load the cubes into memory using the package spectral_cube, which we imported in Step 0. As previously mentioned, the two cubes are at different frequencies (139 GHz and 217 GHz). Once we've loaded the cubes in, we'll use the method find_lines to confirm the presence of the line we expect to find in each of the cubes (H2CS in cube1, SiO in cube2).

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.11.dev10199

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]

Table length=4

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	str1	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		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		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

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]

Table length=4

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	str16	str20	str5	str4	float64	float64	float64	float64	float64	str7	int64	float64	float64	float64	str68	str68	str2	int64	int64	str9	str1	str24	int64	int64	int64	str6	str8	str13	str20	str1	object	int64
20	SiO <font color="red">v = 0 </font>	Silicon Monoxide	5- 4	CDMS		14.4843	21.72614	47.99148	5.0	-3.28429	-1.3211	1	217104.919	20.83956	31.25889	<span style = 'color: #DC143C'></span>	217.104919 (0.002), <span style = 'color: #DC143C'>217.104919</span>	11	-44505	1202			5 0 4 0	16279781	0	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	0	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

We can see four listings for the H2CS 4(1,3)-3(1,2) line (139.483699 GHz) in cube1, each from different linelists. We can also see a few different SiO v=5-4 lines (217.104984 GHz), some from the same linelist (CDMS) and some from the other linelists. These printouts are just for verification purposes; we'll use the center frequencies we have written down in Step 1 for each of these transitions instead of the frequencies given in the find_lines function.

Convert from Jy/beam to kelvin

At this point, we'll convert the cubes to be in units of kelvin instead of Jy/beam. This is because we'll be smoothing the cubes spatially in a future step, and the data values will only be preserved if the cubes are in units of kelvin instead of Jy/beam (since the latter units are dependent on the beamsize).

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

Step 3: Convert cubes from frequency to velocity

Now, we will convert these cubes from frequency units to velocity units using the rest frequencies of the transitions of interest. This is the first step to getting the cubes on a commmon velocity grid: we match up their rest frequencies at 0 km/s. We convert these cubes from frequency to velocity using the method with_spectral_unit, the documentation for which is available here in the spectral_cube API.

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

From the shape of the cube, we can see the H2CS cube is narrower in velocity, so we'll use that as the target spectral reprojection (otherwise, we would have empty velocity channels on either side of the smaller cube that would be getting reprojected). However, the SiO cube is the smaller footprint on the sky.

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

The cubes we are working with are three-dimensional: they have two spatial axes (on the plane of the sky) and one spectral axis. It is hard to visualize them in 3D easily, so we can collapse them into two dimensions in a variety of different ways. Here, we will collapse them along their spectral axis (axis = 0; spectral_cube usually makes the spectral axis the first axis as a default). In particular, for any given spatial coordinate in our cube, we will take the value of the brightest pixel along the spectral axis (the peak flux density). This summary image shows us where the bright emission is (spatially) in the cube, and how bright the emission in a certain region is in relation to the brightness of the emission in another region.

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

WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.max at 0x7f6c8aa46020>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages/spectral_cube/spectral_cube.py:440: RuntimeWarning: All-NaN slice encountered
  out = function(self._get_filled_data(fill=fill,
WARNING: AstropyDeprecationWarning: CoordinateHelper.ticks should not be accessed directly and is deprecated [astropy.visualization.wcsaxes.coordinate_helpers]

INFO: Auto-setting vmin to  1.643e-02 [aplpy.core]
INFO: Auto-setting vmax to  1.839e+00 [aplpy.core]

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

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

/opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages/spectral_cube/spectral_cube.py:440: RuntimeWarning: All-NaN slice encountered
  out = function(self._get_filled_data(fill=fill,
WARNING: AstropyDeprecationWarning: CoordinateHelper.ticks should not be accessed directly and is deprecated [astropy.visualization.wcsaxes.coordinate_helpers]

INFO: Auto-setting vmin to  4.007e-02 [aplpy.core]
INFO: Auto-setting vmax to  2.261e+00 [aplpy.core]

Step 4: Spectral Interpolation

We can safely assume that any two cubes from ALMA are not perfectly matched on a spatial and spectral grid. One cube will have a higher spatial resolution than the other, and one cube will have a finer spectral resolution than the other. We must match these two cubes in two different steps: matching them on a spatial grid (2D), and matching them along their spectral axes (1D). We can treat these as roughly independent steps, and could do either first. However, if your cubes have a varying spatial resolution (different beam size at each spectral channel), you will need to perform spatial smoothing first—you don't want to blend unrelated components of your source by smoothing spectrally before all channels are smoothed to the same spatial resolution. Because each of our cubes have a constant spatial resolution, we will perform spectral smoothing first.

We need to match the spectral resolution of both cubes to each other. We'll start by printing out the spectral resolutions of both cubes.

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>)

cube1, our H2CS cube, has broader channels. We want to smooth cube2 to that channel width before doing the spatial reprojection. To do this, we'll smooth the spectral axis of cube2 using a Gaussian with an appropriate FWHM:

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

$1.8856963 \; \mathrm{\frac{km}{s}}$

Here is the page on smoothing in the spectral_cube documentation. Importantly, here, we're using a Gaussian to smooth the spectral channels, and then we say that the spectral resolution of the resulting cubes is equal to the FWHM of that Gaussian. However, spectral resolution is defined differently by the ALMA correlator: Adam thinks it's half of the channel width smoothed by a Hanning function and then sampled by the channel width, but the answer lies somewhere in the ALMA documentation. The takeaway: we are not necessarily using the "correct" definition of resolution, but we could figure out what that is, if we wanted to.

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)

Then, let's cut down to the same spectral coverage. This is not technically necessary to do, but it is more computationally efficient to smooth over a smaller spectral range. We will do this in a general way, determining the smallest coverage and cutting both cubes based on those values, so no matter which cube has a larger spectral coverage, this code should work.

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.50907945233598 km / s
104.68541315391835 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)

We then smooth with the kernel. Note that this is doing 420x420 = 176400 smoothing operations on a length-221 spectrum: it will take several minutes to run.

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]

Now that we've done spectral smoothing, we can resample the spectrum. We include the parameter suppress_smooth_warning = True in this step, because otherwise we'll see a warning that "the data should be smoothed prior to resampling." Since we have already smoothed the data, this warning can be safely ignored.

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

Spectral Interpolate:   0%|          | 0/176400 [00:00<?, ?it/s]

Spectral Interpolate:   7%|▋         | 12344/176400 [00:00<00:01, 123435.70it/s]

Spectral Interpolate:  14%|█▍        | 24688/176400 [00:00<00:02, 56689.52it/s] 

Spectral Interpolate:  18%|█▊        | 32190/176400 [00:00<00:03, 41859.59it/s]

Spectral Interpolate:  21%|██▏       | 37531/176400 [00:00<00:03, 36070.75it/s]

Spectral Interpolate:  24%|██▎       | 41779/176400 [00:01<00:04, 32688.03it/s]

Spectral Interpolate:  26%|██▌       | 45387/176400 [00:01<00:04, 29855.73it/s]

Spectral Interpolate:  28%|██▊       | 48539/176400 [00:01<00:04, 28273.37it/s]

Spectral Interpolate:  29%|██▉       | 51432/176400 [00:01<00:04, 26717.58it/s]

Spectral Interpolate:  31%|███       | 54114/176400 [00:01<00:05, 23265.95it/s]

Spectral Interpolate:  32%|███▏      | 56459/176400 [00:01<00:05, 22993.01it/s]

Spectral Interpolate:  33%|███▎      | 58762/176400 [00:01<00:05, 22394.19it/s]

Spectral Interpolate:  35%|███▍      | 60992/176400 [00:01<00:05, 22185.33it/s]

Spectral Interpolate:  36%|███▌      | 63200/176400 [00:02<00:05, 21488.08it/s]

Spectral Interpolate:  37%|███▋      | 65336/176400 [00:02<00:05, 21052.54it/s]

Spectral Interpolate:  38%|███▊      | 67428/176400 [00:02<00:05, 19993.39it/s]

Spectral Interpolate:  39%|███▉      | 69420/176400 [00:02<00:05, 19273.62it/s]

Spectral Interpolate:  40%|████      | 71340/176400 [00:02<00:05, 19150.08it/s]

Spectral Interpolate:  42%|████▏     | 73391/176400 [00:02<00:05, 19526.88it/s]

Spectral Interpolate:  43%|████▎     | 75459/176400 [00:02<00:05, 19852.21it/s]

Spectral Interpolate:  44%|████▍     | 77539/176400 [00:02<00:04, 20123.88it/s]

Spectral Interpolate:  45%|████▌     | 79554/176400 [00:02<00:04, 19999.44it/s]

Spectral Interpolate:  46%|████▋     | 81640/176400 [00:03<00:04, 20248.01it/s]

Spectral Interpolate:  47%|████▋     | 83737/176400 [00:03<00:04, 20460.83it/s]

Spectral Interpolate:  49%|████▊     | 85826/176400 [00:03<00:04, 20587.97it/s]

Spectral Interpolate:  50%|████▉     | 87887/176400 [00:03<00:04, 20594.43it/s]

Spectral Interpolate:  51%|█████     | 89962/176400 [00:03<00:04, 20638.47it/s]

Spectral Interpolate:  52%|█████▏    | 92064/176400 [00:03<00:04, 20752.24it/s]

Spectral Interpolate:  53%|█████▎    | 94152/176400 [00:03<00:03, 20787.23it/s]

Spectral Interpolate:  55%|█████▍    | 96244/176400 [00:03<00:03, 20824.44it/s]

Spectral Interpolate:  56%|█████▌    | 98327/176400 [00:03<00:03, 20528.41it/s]

Spectral Interpolate:  57%|█████▋    | 100381/176400 [00:03<00:03, 19934.91it/s]

Spectral Interpolate:  58%|█████▊    | 102379/176400 [00:04<00:04, 17699.69it/s]

Spectral Interpolate:  59%|█████▉    | 104445/176400 [00:04<00:03, 18495.99it/s]

Spectral Interpolate:  60%|██████    | 106548/176400 [00:04<00:03, 19199.37it/s]

Spectral Interpolate:  62%|██████▏   | 108668/176400 [00:04<00:03, 19766.96it/s]

Spectral Interpolate:  63%|██████▎   | 110824/176400 [00:04<00:03, 20282.74it/s]

Spectral Interpolate:  64%|██████▍   | 113056/176400 [00:04<00:03, 20875.83it/s]

Spectral Interpolate:  65%|██████▌   | 115176/176400 [00:04<00:02, 20968.56it/s]

Spectral Interpolate:  66%|██████▋   | 117285/176400 [00:04<00:02, 20468.47it/s]

Spectral Interpolate:  68%|██████▊   | 119343/176400 [00:04<00:02, 20266.78it/s]

Spectral Interpolate:  69%|██████▉   | 121377/176400 [00:05<00:02, 19409.58it/s]

Spectral Interpolate:  70%|███████   | 123539/176400 [00:05<00:02, 20039.08it/s]

Spectral Interpolate:  71%|███████▏  | 125781/176400 [00:05<00:02, 20727.58it/s]

Spectral Interpolate:  73%|███████▎  | 128027/176400 [00:05<00:02, 21233.86it/s]

Spectral Interpolate:  74%|███████▍  | 130377/176400 [00:05<00:02, 21900.02it/s]

Spectral Interpolate:  75%|███████▌  | 132726/176400 [00:05<00:01, 22370.63it/s]

Spectral Interpolate:  77%|███████▋  | 135135/176400 [00:05<00:01, 22879.84it/s]

Spectral Interpolate:  78%|███████▊  | 137639/176400 [00:05<00:01, 23520.59it/s]

Spectral Interpolate:  80%|███████▉  | 140298/176400 [00:05<00:01, 24436.79it/s]

Spectral Interpolate:  81%|████████  | 142809/176400 [00:05<00:01, 24637.93it/s]

Spectral Interpolate:  83%|████████▎ | 145564/176400 [00:06<00:01, 25506.63it/s]

Spectral Interpolate:  84%|████████▍ | 148439/176400 [00:06<00:01, 26474.25it/s]

Spectral Interpolate:  86%|████████▌ | 151469/176400 [00:06<00:00, 27619.24it/s]

Spectral Interpolate:  88%|████████▊ | 154807/176400 [00:06<00:00, 29343.45it/s]

Spectral Interpolate:  90%|████████▉ | 158177/176400 [00:06<00:00, 30647.59it/s]

Spectral Interpolate:  92%|█████████▏| 162795/176400 [00:06<00:00, 35303.13it/s]

Spectral Interpolate:  98%|█████████▊| 173722/176400 [00:06<00:00, 57484.19it/s]

Spectral Interpolate: 100%|██████████| 176400/176400 [00:06<00:00, 26600.92it/s]

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

Step 5: Spatial Smoothing and Reprojection

Recall that our two cubes started with different beam sizes. Our goal is to match the beam sizes of the two cubes, and then to reproject the cubes onto the same spatial pixel grid. We can start to do this by printing out the beam sizes of the two cubes as they currently stand.

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)

cube1 again has the larger beam, so we'll smooth cube2 to its resolution.

Aside: mixed beams

If cube1 and cube2 had different sized beams, but neither was clearly larger, we would have to convolve both to a common beam. In this case, it's redundant and we could have just used cube1's beam, but this is the more general approach:

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

Beam: BMAJ=$1.29719604986604^{''}$ BMIN=$1.04247149438736^{''}$ BPA=$82.95313553702^\circ$

# # 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

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

Now, we can reproject cube2 onto the same spatial grid as cube1.

# 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

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

These two cubes are now on an identical grid, and can be directly compared:

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)

These spectra can now be overplotted as they are in the same unit with the same beam. Using the quicklook maps of the peak intensity from step 2, we know that there is some signal in the central region of our pointing, so we can plot the spectrum for the pixel there. We can approximate the pixel coordinate using the peak intensity map with fewer pixels, since we know both cubes got smoothed to that (coarser) spatial resolution. The center is at about (125, 125).

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]

Do a quick continuum subtraction

If we looked at the spectra we calculated in the previous step and saw that the baseline level of the spectra were not at zero, we would want to perform continuum subtraction to subtract off that non-zero baseline (since we're only working with the spectral lines for our purposes). Though our above spectra look like they are likely already continuum-subtracted, we can perform the following steps anyway. First, we'll calculate the median values of each of the cubes along the spectral axes. This gives us the approximate continuum value in each pixel. The median is more resistant to the presence of the spectral lines than the mean would be.

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

/opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages/spectral_cube/spectral_cube.py:440: RuntimeWarning: All-NaN slice encountered
  out = function(self._get_filled_data(fill=fill,
WARNING: AstropyDeprecationWarning: CoordinateHelper.ticks should not be accessed directly and is deprecated [astropy.visualization.wcsaxes.coordinate_helpers]

INFO: Auto-setting vmin to -2.126e-01 [aplpy.core]
INFO: Auto-setting vmax to  1.949e-01 [aplpy.core]

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

/opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages/spectral_cube/spectral_cube.py:440: RuntimeWarning: All-NaN slice encountered
  out = function(self._get_filled_data(fill=fill,
WARNING: AstropyDeprecationWarning: CoordinateHelper.ticks should not be accessed directly and is deprecated [astropy.visualization.wcsaxes.coordinate_helpers]

INFO: Auto-setting vmin to -9.803e-02 [aplpy.core]
INFO: Auto-setting vmax to  1.078e-01 [aplpy.core]

# 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

# 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]

For our spectra, since they are already continuum-subtracted, there is almost no visible change.

Step 6: Signal Masking

Now, we will use one of the molecules, H2CS, to mask the other one, SiO. This can be useful if SiO has complicated multi-component velocity structure, but H2CS has a similar morphology on the sky and we want to use that to grab the matching component of the SiO velocity structure.

First, we'll estimate the noise throughout the cube by using a signal-resistant measurement of the standard deviation (called the median absolute deviation (MAD) standard deviation). We'll sigma-clip the cube first at 3 sigma and plot the spectrum to look for where the noise is least affected by signal.

# 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 0x7f6c8aa45e40>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]

We can see that between -15 and 45 km/s, the noise estimate is affected by signal. However, above 60 km/s, the noise is flatter. We'll use this part of the spectrum to estimate the noise map.

# 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 0x7f6c8aa45e40>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.mad_std at 0x7f6c8aa45e40>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages/astropy/stats/funcs.py:870: RuntimeWarning: All-NaN slice encountered
  data_median = func(data, axis=axis)
/opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages/astropy/stats/funcs.py:879: RuntimeWarning: All-NaN slice encountered
  result = func(np.abs(data - data_median), axis=axis, overwrite_input=True)

Now, we'll create a "plain mask" where we mask the template cube and only keep pixels that are at or above 3 sigma (where sigma is given on a pixel-by-pixel basis by the noise map we just created). This procedure was originally illustrated in this signal masking tutorial.

# 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)

Then, we'll apply that mask to the template molecule cube and create a more sophisticated "signal mask" from that cube, which takes into account the fact that connected regions of signal are more likely to be real than an individual high-sigma pixel.

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

Binary dilation

The process of binary dilation grows the mask by one pixel outwards in every direction. This is meant to smooth the rough edges of the mask that may be generated when we perform these cuts, and it could help us to recover real contiguous signal that we may have accidentally removed in the masking process. However, it can also pick up noise. User beware: compare the number of outliers in your final moment maps (this is done in a future step) with and without binary dilation to determine if you wish to use it in your processing.

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

Mask application

Finally, we'll apply the signal mask we just created using the template molecule cube (which has a mask at each velocity step) to the target molecule cube.

# 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')

Visualizing the masks:

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()

Step 7: Moment Maps

First, we'll create some initial moment maps of the target molecule using the cube we just masked. The negative values in these moment maps are due to us using a mask from one cube to mask signal in another cube. We can pick up noise that way, and since the noise can be negative, the sum along the spectral axis (which is what we are plotting in a moment 0 map) can also be negative.

# 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.13.5/x64/lib/python3.13/site-packages/spectral_cube/_moments.py:182: RuntimeWarning: invalid value encountered in divide
  return (np.nansum(data * pix_cen, axis=axis) /

The color bar on the moment 1 map is skewed due to the presence of at least one outlier pixel with a large negative velocity.

Step 8: Masking Outliers

At the moment, to get a better idea of the velocity structure of the target molecule, we can mask the outliers by masking pixels with a moment 0 value less than some threshold, which can cause moment 1 values to blow up since the moment 0 value appears in the denominator of the moment 1 value.

# 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.13.5/x64/lib/python3.13/site-packages/spectral_cube/_moments.py:182: 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()

-359492.34390212677 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.13.5/x64/lib/python3.13/site-packages/spectral_cube/_moments.py:182: RuntimeWarning: invalid value encountered in divide
  return (np.nansum(data * pix_cen, axis=axis) /

The velocity structure is more visible here with the outliers removed.

# masked_moment1.write('moment1.fits')

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

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 0x7f6c8aa45800>) 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 0x7f6c8aa45800>) 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()

This is another visualization, but it includes outliers.

# 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. [astropy.utils.decorators]
WARNING: AstropyDeprecationWarning: "max_cut" was deprecated in version 6.1 and will be removed in a future version. Use argument "vmax" instead. [astropy.utils.decorators]