Position-Velocity Diagrams from Disks¶

Authors¶

Adam Ginsburg, Eric Koch

Learning Goals¶

  • Extract a position-velocity diagram from a spectral cube of a protoplanetary disk using pvextractor
  • Extract a position-velocity diagram from a spectral cube of a protoplanetary disk using reproject via spectral-cube using a region mask

Keywords¶

cube, pv-diagram

Summary¶

In this tutorial, we will extract position-velocity (PV) diagrams from a cube of a disk and plot them.

Requirements¶

!pip install --upgrade spectral-cube git+https://github.com/radio-astro-tools/pvextractor.git@61e118aaf28e2d746deeccf06af8fdd7f405b815 radio-beam regions reproject

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
from spectral_cube import SpectralCube
from astropy import wcs
import os

/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 [3]:
import pylab as pl
# set so that these display properly on black backgrounds
pl.rcParams['figure.facecolor']='w'
In [4]:
import radio_beam
In [5]:
from astropy import units as u

We download a 2GB cube from the MAPS survey:

In [6]:
filename = 'HD_163296_CO_220GHz.0.15arcsec.JvMcorr.image.pbcor.fits'
if not os.path.exists(filename):
    try:
        filename = download_file('ftp://ftp.cv.nrao.edu/NRAO-staff/rloomis/MAPS/HD_163296/images/CO/0.15arcsec/HD_163296_CO_220GHz.0.15arcsec.JvMcorr.image.pbcor.fits', cache=True, timeout=10)
    except:
        import ftplib
        ftp = ftplib.FTP('ftp.cv.nrao.edu')
        ftp.login()
        ftp.cwd('NRAO-staff/rloomis/MAPS/HD_163296/images/CO/0.15arcsec')
        with open('HD_163296_CO_220GHz.0.15arcsec.JvMcorr.image.pbcor.fits', 'wb') as fp:
            ftp.retrbinary('RETR HD_163296_CO_220GHz.0.15arcsec.JvMcorr.image.pbcor.fits', fp.write)
        ftp.quit()
In [7]:
# import ftplib
# ftp = ftplib.FTP('ftp.cv.nrao.edu')
# ftp.login()
# ftp.cwd('NRAO-staff/rloomis/MAPS/HD_163296/images/CO/0.15arcsec')
# with open('HD_163296_CO_220GHz.0.15arcsec.JvMcorr.image.pbcor.fits', 'wb') as fp:
#     ftp.retrbinary('RETR HD_163296_CO_220GHz.0.15arcsec.JvMcorr.image.pbcor.fits', fp.write)
# ftp.quit()

We load the cube using the dask backend, which allows for some parallelization:

In [8]:
cube = SpectralCube.read(filename, use_dask=True)
In [9]:
cube
Out[9]:
DaskSpectralCube with shape=(127, 2048, 2048) and unit=Jy / beam and chunk size (127, 514, 514):
 n_x:   2048  type_x: RA---SIN  unit_x: deg    range:   269.082528 deg:  269.094790 deg
 n_y:   2048  type_y: DEC--SIN  unit_y: deg    range:   -21.961977 deg:  -21.950605 deg
 n_s:    127  type_s: FREQ      unit_s: Hz     range: 230523958206.200 Hz:230543336804.442 Hz
In [10]:
mx = cube.max(axis=0)

A quick look at the image cube shows that there is a disk rotated about 45 degrees in the center of the frame:

In [11]:
pl.imshow(mx.value, origin='lower')
Out[11]:
<matplotlib.image.AxesImage at 0x7f76f215c2f0>
No description has been provided for this image

We can draw an ellipse around the disk to downselect only it:

In [12]:
import regions
In [13]:
center = regions.PixCoord(1024, 1024)
ellipse = regions.EllipsePixelRegion(center, width=550, height=400, angle=45*u.deg)
In [14]:
ax = pl.gca()
ax.imshow(mx.value, origin='lower')
ellipse.plot(ax=ax, facecolor='none', edgecolor='red', lw=2)
Out[14]:
<matplotlib.patches.Ellipse at 0x7f76f100c2f0>
No description has been provided for this image

We make a cutout by creating a subcube using the ellipse region as a mask:

In [15]:
cutout = cube.subcube_from_regions([ellipse])
cutout
Out[15]:
DaskSpectralCube with shape=(127, 481, 481) and unit=Jy / beam and chunk size (127, 244, 244):
 n_x:    481  type_x: RA---SIN  unit_x: deg    range:   269.087219 deg:  269.090094 deg
 n_y:    481  type_y: DEC--SIN  unit_y: deg    range:   -21.957622 deg:  -21.954955 deg
 n_s:    127  type_s: FREQ      unit_s: Hz     range: 230523958206.200 Hz:230543336804.442 Hz

Then we want to extract a position-velocity diagram across the disk.

We specify a width of 200 pixels (we could go to ~400) so we average across the short axis of the disk:

In [16]:
import pvextractor
In [17]:
path = pvextractor.Path([(0,0), (481,481)], width=200)

We show the path overlaid on our cutout disk:

In [18]:
ax = pl.subplot(111, projection=cutout.wcs.celestial)
ax.imshow(cutout.max(axis=0).value, origin='lower')
path.show_on_axis(ax, spacing=5, alpha=0.7, linewidth=0.25)
Out[18]:
<matplotlib.collections.PatchCollection at 0x7f76eddf06e0>
No description has been provided for this image

Then, we extract the PV diagram. We choose spacing=5 to average over 5 pixels. This averaging isn't necessary, but does make the operation a little faster and increases the signal-to-noise ratio per spatial bin.

In [19]:
pv = pvextractor.extract_pv_slice(cutout.with_spectral_unit(u.km/u.s, velocity_convention='radio'), path, spacing=5)

























And plot the resulting diagram:


















In [20]:






ax = pl.subplot(111, projection=wcs.WCS(pv.header))
im = ax.imshow(pv.data)
cb = pl.colorbar(mappable=im)
cb.set_label("Brightness Temperature [K]")
ax.set_aspect(1)
























No description has been provided for this image























Second approach¶
We can also reproject the whole cube by rotating 45 degrees.


This requires making our own new header, which is a bit tedious, but effective.


















In [21]:






header = cutout.wcs.to_header()
header['NAXIS'] = 3
header['NAXIS1'] = 600
header['NAXIS2'] = 400
header['NAXIS3'] = cutout.shape[0]
angle = 45*u.deg
header['CD1_1'] = np.cos(angle).value * np.abs(cube.wcs.wcs.cdelt[0])
header['CD2_1'] = -np.sin(angle).value * np.abs(cube.wcs.wcs.cdelt[0])
header['CD1_2'] = np.sin(angle).value * np.abs(cube.wcs.wcs.cdelt[1])
header['CD2_2'] = np.cos(angle).value * np.abs(cube.wcs.wcs.cdelt[1])
header['CD3_3'] = cube.wcs.wcs.cdelt[2]
header['CRPIX1'] = 300
header['CRPIX2'] = 200

























We then reproject the whole cube, which takes a minute or two:


















In [22]:






reproj = cutout.reproject(header)





















In [23]:






reproj





















Out[23]:




DaskSpectralCube with shape=(127, 400, 600) and unit=Jy / beam and chunk size (127, 363, 363):
 n_x:    600  type_x: RA---SIN  unit_x: deg    range:   269.086547 deg:  269.090774 deg
 n_y:    400  type_y: DEC--SIN  unit_y: deg    range:   -21.958249 deg:  -21.954328 deg
 n_s:    127  type_s: FREQ      unit_s: Hz     range: 230523958206.200 Hz:230543336804.442 Hz




















In [24]:






rmax = reproj.max(axis=0)





















In [25]:






rmax.quicklook()

























WARNING: AstropyDeprecationWarning: CoordinateHelper.ticks should not be accessed directly and is deprecated [astropy.visualization.wcsaxes.coordinate_helpers]













INFO: Auto-setting vmin to -2.397e-03 [aplpy.core]
INFO: Auto-setting vmax to  6.239e-02 [aplpy.core]












No description has been provided for this image























Then, the position-velocity diagram is easy: we just take the average along the short axis:


















In [26]:






pv2 = reproj.with_spectral_unit(u.km/u.s, velocity_convention='radio').mean(axis=1)

























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]





















In [27]:






ax = pl.subplot(111, projection=wcs.WCS(pv2.header))
im = ax.imshow(pv2.data)
cb = pl.colorbar(mappable=im)
cb.set_label("Brightness Temperature [K]")
ax.set_aspect(4)
























No description has been provided for this image



















In [ ]: