Cube Reprojection Tutorial¶

Authors¶

Adam Ginsburg, Eric Koch

Learning Goals¶

  • reproject a cube spectrally
  • smooth it spectrally
  • reproject it spatially

Keywords¶

cube, reprojection

Summary¶

This tutorial shows how to take two spectral cubes observed toward the same part of the sky, but different frequencies, and put them onto the same grid using spectral-cube.

Index¶

  • Step 1: Download
  • Step 2: Open files, collect metadata
  • Step 3: Convert to velocity
  • Step 4: Spectral Interpolation
  • Step 5: Spatial Smoothing
  • Step 6: Reprojection

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. For more information see the spectral-cube documentation.

Step 1: Download the data¶

(you might not have to do this step, since you may already have data)

In [1]:
%pip install spectral-cube dask 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)
Requirement already satisfied: dask in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (2025.5.1)
Requirement already satisfied: numpy in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (2.3.1)
Requirement already satisfied: matplotlib in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (3.10.3)
Requirement already satisfied: astropy in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (7.1.0)
Requirement already satisfied: casa-formats-io>=0.1 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from spectral-cube) (0.3.0)
Requirement already satisfied: joblib>=1.1 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from spectral-cube) (1.5.1)
Requirement already satisfied: packaging>=19 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from spectral-cube) (25.0)
Requirement already satisfied: radio-beam>=0.3.5 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from spectral-cube) (0.3.10.dev13+gf75024c)
Requirement already satisfied: setuptools>=62.3.3 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from spectral-cube) (80.9.0)
Requirement already satisfied: tqdm>=4.64 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from spectral-cube) (4.67.1)
Requirement already satisfied: click>=8.1 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from dask) (8.2.1)
Requirement already satisfied: cloudpickle>=3.0.0 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from dask) (3.1.1)
Requirement already satisfied: fsspec>=2021.09.0 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from dask) (2025.5.1)
Requirement already satisfied: partd>=1.4.0 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from dask) (1.4.2)
Requirement already satisfied: pyyaml>=5.3.1 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from dask) (6.0.2)
Requirement already satisfied: toolz>=0.10.0 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from dask) (1.0.0)
Requirement already satisfied: contourpy>=1.0.1 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from matplotlib) (1.3.2)
Requirement already satisfied: cycler>=0.10 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from matplotlib) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from matplotlib) (4.58.4)
Requirement already satisfied: kiwisolver>=1.3.1 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from matplotlib) (1.4.8)
Requirement already satisfied: pillow>=8 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from matplotlib) (11.2.1)
Requirement already satisfied: pyparsing>=2.3.1 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from matplotlib) (3.2.3)
Requirement already satisfied: python-dateutil>=2.7 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from matplotlib) (2.9.0.post0)
Requirement already satisfied: pyerfa>=2.0.1.1 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from astropy) (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) (0.2025.6.30.0.39.40)
Requirement already satisfied: locket in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from partd>=1.4.0->dask) (1.0.0)
Requirement already satisfied: six>=1.5 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from python-dateutil>=2.7->matplotlib) (1.17.0)
Requirement already satisfied: scipy>=1.8 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from radio-beam>=0.3.5->spectral-cube) (1.16.0)
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.
In [2]:
import numpy as np
from astropy.utils.data import download_file

We download the data cubes (18 MB and 337 MB, respectively) from a permalink on the ALMA archives.

If you have trouble with these downloads, try changing to a different ALMA server (e.g., almascience.nrao.edu->almascience.eso.org) or increase the timeout. See https://docs.astropy.org/en/stable/api/astropy.utils.data.download_file.html.

In [3]:
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)
In [4]:
filename_2 = download_file("https://almascience.nrao.edu/dataPortal/member.uid___A001_X87d_X141.a_sma1_sci.spw27.cube.I.pbcor.fits",
                          cache=True)

Step 2: Load the cubes¶

In [5]:
from spectral_cube import SpectralCube
/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
In [6]:
cube1 = SpectralCube.read(filename_1)
cube1
WARNING: StokesWarning: Cube is a Stokes cube, returning spectral cube for I component [spectral_cube.io.core]
Out[6]:
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
In [7]:
cube2 = SpectralCube.read(filename_2)
cube2
WARNING: StokesWarning: Cube is a Stokes cube, returning spectral cube for I component [spectral_cube.io.core]
Out[7]:
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

The cubes are at different frequencies - 139 and 89 GHz.

The first cube covers the H2CS 4(1,3)-3(1,2) line at 139.483699 GHz.

The second covers SiO v=5-4 at 217.104984 GHz

We use the find_lines tool to query splatalogue with astroquery over the spectral range covered by the cube. It returns a table of matching lines. Note that some line names will be repeated because Splatalogue includes several different databases and most chemical species are present in all of these.

In [8]:
# DEBUG
import astroquery
print(astroquery.__version__)
0.4.11.dev10199
In [9]:
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]
Out[9]:
Table length=4
species_idnamechemical_nameresolved_QNslinelistLovasASTIntensitylower_state_energyupper_state_energysijmu2sijaijintintensityLovas_NRAOorderedfreqlower_state_energy_Kupper_state_energy_KorderedFreqmeasFrequpperStateDegenmoleculeTagqnCodelabref_Lovas_NISTrel_int_HFS_Lovasunres_quantum_numberslineidtransition_in_spacetransition_in_G358obsref_Lovas_NISTsource_Lovas_NISTtelescope_Lovas_NISTtransitionBandColorsearchErrorMessagesqlqueryrequestnumber
int64str15str16str28str5str5float64float64float64float64float64str1int64float64float64float64str68str69str3int64int64str8str1str28int64int64int64str6str11str13str14str1objectint64
210H<sub>2</sub>CSThioformaldehyde4( 1, 3)- 3( 1, 2)JPL0.1716.146520.7991730.5167411.223-4.447321139483.4123.2310829.92519<span style = 'color: #DC143C'></span>139.48341 (0.28), <span style = 'color: #DC143C'>139.48341</span>27-460033034 1 3 3 1 2406568700Lor84arho Oph B1MMWO 4.9mdatatablegreenNone0
210H<sub>2</sub>CSThioformaldehyde4( 1, 3)- 3( 1, 2)CDMS16.132920.7855730.5966111.251-4.446190139483.681623.2115129.90563<span style = 'color: #DC143C'></span>139.4836816 (0.05), <span style = 'color: #DC143C'>139.4836816</span>27-465093034 1 3 3 1 2406626600datatablegreenNone0
210H<sub>2</sub>CSThioformaldehyde4( 1, 3)- 3( 1, 2)SLAIM0.1716.13320.7856830.594723.75-3.944260139483.69923.2116529.90578139.483699 (0.017), <span style = 'color: #DC143C'>139.483699</span>139.483741 (0.024), <span style = 'color: #DC143C'>139.483741</span>8.500Mar054( 1, 3)- 3( 1, 2)374448710datatablegreenNone0
210H<sub>2</sub>CSThioformaldehyde4(1,3)-3(1,2)Lovas0.170.00.00.00.00.00139483.6990.00.0139.483699 (0.017), <span style = 'color: #DC143C'>139.483699</span><span style = 'color: #DC143C'></span>004(1,3)-3(1,2)373419900Lor84arho Oph B1MMWO 4.9mdatatablegreenNone0
In [10]:
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]
Out[10]:
Table length=4
species_idnamechemical_nameresolved_QNslinelistLovasASTIntensitylower_state_energyupper_state_energysijmu2sijaijintintensityLovas_NRAOorderedfreqlower_state_energy_Kupper_state_energy_KorderedFreqmeasFrequpperStateDegenmoleculeTagqnCodelabref_Lovas_NISTrel_int_HFS_Lovasunres_quantum_numberslineidtransition_in_spacetransition_in_G358obsref_Lovas_NISTsource_Lovas_NISTtelescope_Lovas_NISTtransitionBandColorsearchErrorMessagesqlqueryrequestnumber
int64str35str16str20str5str4float64float64float64float64float64str7int64float64float64float64str68str68str2int64int64str9str1str24int64int64int64str6str8str13str20str1objectint64
20SiO <font color="red">v = 0 </font>Silicon Monoxide5- 4CDMS14.484321.7261447.991485.0-3.28429-1.32111217104.91920.8395631.25889<span style = 'color: #DC143C'></span>217.104919 (0.002), <span style = 'color: #DC143C'>217.104919</span>11-4450512025 0 4 01627978100datatablelightpurpleNone0
20SiO <font color="red">v = 0 </font>Silicon Monoxide5-4JPL1.614.484321.7261448.146515.0-3.28288-1.31660217104.9820.8395631.25889<span style = 'color: #DC143C'></span>217.10498 (0.08), <span style = 'color: #DC143C'>217.10498</span>11-44002101J=5-4 v=0116736800Lor84aOriMC-1MMWO 4.9mdatatablelightpurpleNone0
20SiO <font color="red">v = 0 </font>Silicon Monoxide5-4Lovas1.60.00.00.00.00.00217104.9840.00.0217.104984 (0.014), <span style = 'color: #DC143C'>217.104984</span><span style = 'color: #DC143C'></span>005-4 v=0372463800Lor84aOriMC-1MMWO 4.9mdatatablelightpurpleNone0
20SiO <font color="red">v = 0 </font>Silicon Monoxide5 - 4SLAIM14.48421.7258447.68495.0-3.287060217104.98420.8391331.25846217.104984 (0.014), <span style = 'color: #DC143C'>217.104984</span>217.10498 (0.1), <span style = 'color: #DC143C'>217.10498</span>1100Man775 - 4 v=0393286600datatablelightpurpleNone0

Step 3: Convert cubes from frequency to velocity¶

In [11]:
from astropy import units as u
In [12]:
cube1vel = cube1.with_spectral_unit(u.km/u.s, velocity_convention='radio', rest_value=139.483699*u.GHz)
cube1vel
Out[12]:
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: VRAD      unit_s: km / s  range:      -43.509 km / s:     104.685 km / s
In [13]:
cube2vel = cube2.with_spectral_unit(u.km/u.s, velocity_convention='radio', rest_value=217.104984*u.GHz)
cube2vel
Out[13]:
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: 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. However, the SiO cube is the smaller footprint on the sky.

Create spatial maps of the peak intensity to quickly explore the cubes:¶

One way to quickly explore the structure in the data cubes is to produce a peak intensity map, or the maximum along the spectral axis (axis=0).

In [14]:
mx = cube1.max(axis=0)
mx.quicklook()
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.max at 0x7f9be0c89b20>) 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  3.023e-04 [aplpy.core]
INFO: Auto-setting vmax to  3.936e-02 [aplpy.core]
No description has been provided for this image

We can do the same thing all on one line (for the other cube this time):

In [15]:
cube2.max(axis=0).quicklook()
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.max at 0x7f9be0c89b20>) 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  9.465e-04 [aplpy.core]
INFO: Auto-setting vmax to  5.172e-02 [aplpy.core]
No description has been provided for this image

Step 4. Spectral Interpolation¶

We can do the spatial or spectral step first. In this case, we choose the spectral step first because the H$_2$CS cube is narrower in velocity (cube1vel) and this will reduce the number of channels we need to spatially interpolate over in the next step.

We need to match resolution to the cube with the largest channel width:

In [16]:
velocity_res_1 = np.diff(cube1vel.spectral_axis)[0]
velocity_res_2 = np.diff(cube2vel.spectral_axis)[0]
velocity_res_1, velocity_res_2
Out[16]:
(<Quantity 2.00262828 km / s>, <Quantity 0.67429189 km / s>)

Next, we will reduce cube2vel to have the same spectral range as cube1vel:

In [17]:
cube2vel_cutout = cube2vel.spectral_slab(cube1vel.spectral_axis.min(),
                                         cube1vel.spectral_axis.max())
cube1vel, cube2vel_cutout
Out[17]:
(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: VRAD      unit_s: km / s  range:      -43.509 km / s:     104.685 km / s,
 SpectralCube with shape=(221, 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:    221  type_s: VRAD      unit_s: km / s  range:      -43.432 km / s:     104.913 km / s)

Note that it is important for the to-be-interpolated cube, in this case cube2, to have pixels bounding cube1's spectral axis, but in this case it does not. If the pixel range doesn't overlap perfectly, it may blank out one of the edge pixels. So, to fix this, we add a little buffer:

In [18]:
cube2vel_cutout = cube2vel.spectral_slab(cube1vel.spectral_axis.min() - velocity_res_2,
                                         cube1vel.spectral_axis.max())
cube1vel, cube2vel_cutout
Out[18]:
(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: VRAD      unit_s: km / s  range:      -43.509 km / s:     104.685 km / s,
 SpectralCube with shape=(222, 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:    222  type_s: VRAD      unit_s: km / s  range:      -44.106 km / s:     104.913 km / s)

Our H2CS cube (cube1vel) has broader channels. We need to first smooth cube2vel to the broader channel width before doing the spatial reprojection.

To do this, we will spectrally smooth with a Gaussian with width set such that smoothing cube2vel will result in the same width as cube1vel. We do this by finding the difference in widths when deconvolving the cube1vel channel width from cube2vel. For further information see the documentation on smoothing.

Note that if we did not do this smoothing step, we would under-sample the cube2vel data in the next downsampling step, reducing our signal-to-noise ratio.

We have adopted a width equal to the channel width; the line spread function is actually a Hanning-smoothed tophat. We are making a coarse approximation here.

In [19]:
fwhm_gaussian = (velocity_res_1**2 - velocity_res_2**2)**0.5
fwhm_gaussian
Out[19]:
$1.8856963 \; \mathrm{\frac{km}{s}}$
In [20]:
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)

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

In [21]:
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 spectral axis of cube2vel_smooth to match cube1vel by interpolating cube2vel_smooth onto cube1vel's grid:

In [22]:
cube2vel_spectralresample = cube2vel_smooth.spectral_interpolate(cube1vel.spectral_axis,
                                                                 suppress_smooth_warning=True)
cube2vel_spectralresample
Spectral Interpolate:   0%|          | 0/176400 [00:00<?, ?it/s]
Spectral Interpolate:   7%|▋         | 11976/176400 [00:00<00:01, 119726.72it/s]
Spectral Interpolate:  14%|█▎        | 23949/176400 [00:00<00:02, 56304.43it/s] 
Spectral Interpolate:  18%|█▊        | 31295/176400 [00:00<00:03, 43338.38it/s]
Spectral Interpolate:  21%|██        | 36655/176400 [00:00<00:03, 35819.76it/s]
Spectral Interpolate:  23%|██▎       | 40849/176400 [00:01<00:04, 32904.26it/s]
Spectral Interpolate:  25%|██▌       | 44459/176400 [00:01<00:04, 30230.87it/s]
Spectral Interpolate:  27%|██▋       | 47637/176400 [00:01<00:04, 28730.91it/s]
Spectral Interpolate:  29%|██▊       | 50572/176400 [00:01<00:04, 26391.75it/s]
Spectral Interpolate:  30%|███       | 53225/176400 [00:01<00:04, 24878.34it/s]
Spectral Interpolate:  32%|███▏      | 55697/176400 [00:01<00:05, 22194.84it/s]
Spectral Interpolate:  33%|███▎      | 57908/176400 [00:01<00:05, 19782.66it/s]
Spectral Interpolate:  34%|███▍      | 60032/176400 [00:01<00:05, 20108.82it/s]
Spectral Interpolate:  35%|███▌      | 62064/176400 [00:02<00:05, 20157.48it/s]
Spectral Interpolate:  36%|███▋      | 64315/176400 [00:02<00:05, 20761.35it/s]
Spectral Interpolate:  38%|███▊      | 66408/176400 [00:02<00:05, 20786.27it/s]
Spectral Interpolate:  39%|███▉      | 68543/176400 [00:02<00:05, 20940.64it/s]
Spectral Interpolate:  40%|████      | 70647/176400 [00:02<00:05, 20542.00it/s]
Spectral Interpolate:  41%|████      | 72708/176400 [00:02<00:05, 20484.67it/s]
Spectral Interpolate:  42%|████▏     | 74826/176400 [00:02<00:04, 20680.30it/s]
Spectral Interpolate:  44%|████▎     | 76937/176400 [00:02<00:04, 20805.42it/s]
Spectral Interpolate:  45%|████▍     | 79021/176400 [00:02<00:04, 20656.57it/s]
Spectral Interpolate:  46%|████▌     | 81089/176400 [00:03<00:04, 19748.20it/s]
Spectral Interpolate:  47%|████▋     | 83092/176400 [00:03<00:04, 19828.16it/s]
Spectral Interpolate:  48%|████▊     | 85082/176400 [00:03<00:04, 19559.87it/s]
Spectral Interpolate:  49%|████▉     | 87123/176400 [00:03<00:04, 19805.68it/s]
Spectral Interpolate:  51%|█████     | 89218/176400 [00:03<00:04, 20138.53it/s]
Spectral Interpolate:  52%|█████▏    | 91236/176400 [00:03<00:04, 18503.10it/s]
Spectral Interpolate:  53%|█████▎    | 93277/176400 [00:03<00:04, 19033.60it/s]
Spectral Interpolate:  54%|█████▍    | 95204/176400 [00:03<00:04, 18911.96it/s]
Spectral Interpolate:  55%|█████▌    | 97230/176400 [00:03<00:04, 19298.14it/s]
Spectral Interpolate:  56%|█████▌    | 99173/176400 [00:03<00:03, 19323.14it/s]
Spectral Interpolate:  57%|█████▋    | 101277/176400 [00:04<00:03, 19825.89it/s]
Spectral Interpolate:  59%|█████▊    | 103267/176400 [00:04<00:03, 19101.95it/s]
Spectral Interpolate:  60%|█████▉    | 105188/176400 [00:04<00:04, 17715.54it/s]
Spectral Interpolate:  61%|██████    | 107287/176400 [00:04<00:03, 18620.44it/s]
Spectral Interpolate:  62%|██████▏   | 109280/176400 [00:04<00:03, 18989.78it/s]
Spectral Interpolate:  63%|██████▎   | 111199/176400 [00:04<00:03, 18700.25it/s]
Spectral Interpolate:  64%|██████▍   | 113302/176400 [00:04<00:03, 19369.93it/s]
Spectral Interpolate:  65%|██████▌   | 115417/176400 [00:04<00:03, 19887.10it/s]
Spectral Interpolate:  67%|██████▋   | 117648/176400 [00:04<00:02, 20597.49it/s]
Spectral Interpolate:  68%|██████▊   | 119717/176400 [00:05<00:02, 20290.48it/s]
Spectral Interpolate:  69%|██████▉   | 121930/176400 [00:05<00:02, 20827.78it/s]
Spectral Interpolate:  70%|███████   | 124191/176400 [00:05<00:02, 21352.88it/s]
Spectral Interpolate:  72%|███████▏  | 126545/176400 [00:05<00:02, 21998.71it/s]
Spectral Interpolate:  73%|███████▎  | 128864/176400 [00:05<00:02, 22351.00it/s]
Spectral Interpolate:  74%|███████▍  | 131103/176400 [00:05<00:02, 20251.94it/s]
Spectral Interpolate:  76%|███████▌  | 133457/176400 [00:05<00:02, 21164.63it/s]
Spectral Interpolate:  77%|███████▋  | 135949/176400 [00:05<00:01, 22231.71it/s]
Spectral Interpolate:  78%|███████▊  | 138258/176400 [00:05<00:01, 22479.25it/s]
Spectral Interpolate:  80%|███████▉  | 140840/176400 [00:05<00:01, 23452.63it/s]
Spectral Interpolate:  81%|████████▏ | 143497/176400 [00:06<00:01, 24369.14it/s]
Spectral Interpolate:  83%|████████▎ | 146352/176400 [00:06<00:01, 25604.83it/s]
Spectral Interpolate:  85%|████████▍ | 149312/176400 [00:06<00:01, 26788.58it/s]
Spectral Interpolate:  86%|████████▋ | 152328/176400 [00:06<00:00, 27790.67it/s]
Spectral Interpolate:  88%|████████▊ | 155916/176400 [00:06<00:00, 30205.80it/s]
Spectral Interpolate:  91%|█████████ | 159785/176400 [00:06<00:00, 32739.45it/s]
Spectral Interpolate:  93%|█████████▎| 164444/176400 [00:06<00:00, 36883.94it/s]
Spectral Interpolate: 100%|██████████| 176400/176400 [00:06<00:00, 26185.72it/s]

Out[22]:
SpectralCube with shape=(75, 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:     75  type_s: VRAD      unit_s: km / s  range:      -43.509 km / s:     104.685 km / s

Note that we included the suppress_smooth_warning=True argument. That is to hide this warning:

WARNING: SmoothingWarning: Input grid has too small a spacing. The data should be smoothed prior to resampling. [spectral_cube.spectral_cube]

which will tell you if the operation will under-sample the original data. The smoothing work we did above is specifically to make sure we are properly sampling, so this warning does not apply.

Step 5. Spatial Smoothing¶

Now that we've done spectral smoothing, we also need to follow a similar procedure of smoothing then resampling for the spatial axes.

The beam is the resolution element of our cubes:

In [23]:
cube1vel.beam, cube2vel_spectralresample.beam
Out[23]:
(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 hase 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:

In [24]:
import radio_beam
common_beam = cube1vel.beam.commonbeam_with(cube2vel.beam)

# This works with older versions of radio-beam:
# common_beam = radio_beam.commonbeam.common_2beams(radio_beam.Beams(beams=[cube1vel.beam, cube2vel.beam]))

common_beam
Out[24]:
Beam: BMAJ=$1.29719604986604^{''}$ BMIN=$1.04247149438736^{''}$ BPA=$82.95313553702^\circ$

We then convolve:

In [25]:
# for v<0.6, we convert to Kelvin to ensure the units are preserved:
# cube2vel_spatialspectralsmooth = cube2vel_spectralresample.to(u.K).convolve_to(common_beam)
# in more recent versions, the unit conversion is handled appropriately,
# so unit conversion isn't needed
cube2vel_spatialspectralsmooth = cube2vel_spectralresample.convolve_to(common_beam)
cube2vel_spatialspectralsmooth
Out[25]:
SpectralCube with shape=(75, 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:     75  type_s: VRAD      unit_s: km / s  range:      -43.509 km / s:     104.685 km / s

Step 6. Reprojection¶

Now we can do the spatial resampling as the final step for producing two cubes matched to the same spatial and spectral pixel grid:

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

cube2vel_reproj = cube2vel_spatialspectralsmooth.reproject(tgt_header)
cube2vel_reproj
Out[26]:
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: 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:

In [27]:
cube2vel_reproj, cube1vel
Out[27]:
(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: VRAD      unit_s: km / s  range:      -43.509 km / s:     104.685 km / s,
 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: 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.

In [28]:
cube1vel[:,125,125].quicklook()
cube2vel_reproj[:,125,125].quicklook()
No description has been provided for this image

Dask¶

All of the above can be done using dask as the underlying framework to parallelize the operations.

The dask approach can be made more memory-efficient (avoid using too much RAM) by writing intermediate steps to disk. The non-dask approach used above will generally need to read the whole cube into memory. Depending on the situation, either approach may be faster, but dask may be needed if the cube is larger than memory.

We repeat all the operations above using dask. We use a ProgressBar so you can see how long it takes. We also suppress warnings to make the output look cleaner (we already saw all the important warnings above).

In [29]:
from dask.diagnostics import ProgressBar
import warnings
In [30]:
with warnings.catch_warnings():
    warnings.simplefilter('ignore')
    with ProgressBar():
        cube2dask = SpectralCube.read(filename_2, use_dask=True)
        cube2daskvel = cube2dask.with_spectral_unit(u.km/u.s,
                                            velocity_convention='radio', rest_value=217.104984*u.GHz)
        cube2daskvel_cutout = cube2daskvel.spectral_slab(cube1vel.spectral_axis.min() - velocity_res_2,
                                                 cube1vel.spectral_axis.max())
        cube2daskvel_smooth = cube2daskvel_cutout.spectral_smooth(spectral_smoothing_kernel)
        cube2daskvel_spectralresample = cube2daskvel_smooth.spectral_interpolate(cube1vel.spectral_axis,
                                                                             suppress_smooth_warning=True)
        cube2daskvel_spatialspectralsmooth = cube2daskvel_spectralresample.convolve_to(common_beam)
        cube2daskvel_reproj = cube2daskvel_spatialspectralsmooth.reproject(tgt_header)
cube2daskvel_reproj
[                                        ] | 0% Completed | 246.49 us
[                                        ] | 0% Completed | 105.17 ms
[                                        ] | 0% Completed | 426.34 ms
[                                        ] | 0% Completed | 698.16 ms
[                                        ] | 0% Completed | 823.82 ms
[                                        ] | 0% Completed | 1.63 s
[                                        ] | 0% Completed | 1.91 s
[                                        ] | 0% Completed | 2.18 s
[                                        ] | 0% Completed | 2.28 s
[                                        ] | 0% Completed | 2.38 s
[                                        ] | 0% Completed | 2.48 s
[                                        ] | 0% Completed | 3.29 s
[                                        ] | 0% Completed | 3.39 s
[                                        ] | 0% Completed | 3.49 s
[                                        ] | 0% Completed | 3.59 s
[                                        ] | 0% Completed | 153.31 s
[                                        ] | 0% Completed | 153.41 s
[                                        ] | 0% Completed | 153.52 s
[                                        ] | 0% Completed | 153.76 s
[                                        ] | 0% Completed | 154.01 s
[                                        ] | 0% Completed | 154.85 s
[                                        ] | 0% Completed | 154.95 s
[                                        ] | 0% Completed | 155.06 s
[########################################] | 100% Completed | 155.16 s

Out[30]:
DaskSpectralCube with shape=(75, 250, 250) and unit=Jy / beam and chunk size (75, 250, 250):
 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
In [ ]: