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_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 |
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_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 |
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]
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]
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()
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 [ ]: