Extracting and Plotting Position-Velocity Diagrams¶

Authors¶

Adam Ginsburg, Eric Koch

Learning Goals¶

  • Extract a position-velocity diagram from a cube using both pixel and sky coordinates using pvextractor
  • Display the position-velocity diagram with appropriately labeled coordinates
  • Display the extraction path on the plots

Keywords¶

cube, pv-diagram

Summary¶

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

Header material¶

We import tools from several packages up front:

In [1]:
%pip install spectral-cube dask numpy matplotlib astropy pvextractor
%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: pvextractor in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (0.5.dev30+gd182280)
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: scipy>=1.8 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from pvextractor) (1.16.0)
Requirement already satisfied: qtpy>=2.0 in /opt/hostedtoolcache/Python/3.13.5/x64/lib/python3.13/site-packages (from pvextractor) (2.4.3)

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)

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

import pylab as pl
from astropy.visualization import quantity_support
from astropy import units as u
from astropy import wcs

# set so that these display properly on black backgrounds
pl.rcParams['figure.facecolor']='w'

from spectral_cube import SpectralCube

from pvextractor import extract_pv_slice, Path

/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

Retrieve and open a cube from astropy-data:

In [3]:
cube = SpectralCube.read('http://www.astropy.org/astropy-data/l1448/l1448_13co.fits')
cube
Out[3]:
SpectralCube with shape=(53, 105, 105):
 n_x:    105  type_x: RA---SFL  unit_x: deg    range:    50.924417 deg:   51.740103 deg
 n_y:    105  type_y: DEC--SFL  unit_y: deg    range:    30.301945 deg:   30.966389 deg
 n_s:     53  type_s: VOPT      unit_s: m / s  range:     2528.195 m / s:    5982.223 m / s

Show single channel to find where to draw the path. We use pixel units so it's easier to define the path from the pixel coords in the matplotlib viewer.

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

PV Extraction from Pixel Coordinates¶

First we create an extraction path. The entries are pairs of pixel coordinates, (x,y)

In [5]:
path = Path([(20,20), (40,40), (60,20)])

Then we can overplot it on our figure, now with WCS shown. The plotting uses WCSAxes

In [6]:
ax = pl.subplot(111, projection=cube.wcs.celestial)
ax.imshow(cube[25].value)
path.show_on_axis(ax, spacing=1, color='w')
Out[6]:
<matplotlib.lines.Line2D at 0x7fbcbfcaafd0>
No description has been provided for this image

spacing gives the separation between these points in pixels; we finely sampled by picking one-pixel spacing.

We can then extract the pv diagram, specifying the same spacing.

In [7]:
pvdiagram = extract_pv_slice(cube=cube, path=path, spacing=1)
pvdiagram
Out[7]:
<astropy.io.fits.hdu.image.PrimaryHDU at 0x7fbcbfbc7750>

and plot it. pvdiagram is a PrimaryHDU object, so we need to grab the data separately from the header and convert the header to a WCS object:

In [8]:
ax = pl.subplot(111, projection=wcs.WCS(pvdiagram.header))
im = ax.imshow(pvdiagram.data)
cb = pl.colorbar(mappable=im)
# we could specify the colorbar units like this:
# cb.set_label(cube.unit)
# but the 'BUNIT' keyword is not set for these data, so we don't know the unit.  We instead manually specify:
cb.set_label("Brightness Temperature [K]")
No description has been provided for this image

Changing units to the more commonly used km/s and more readable arcminutes can be done with wcsaxes tools:

In [9]:
ww = wcs.WCS(pvdiagram.header)

ax = pl.subplot(111, projection=ww)
im = ax.imshow(pvdiagram.data)
cb = pl.colorbar(mappable=im)
cb.set_label("Brightness Temperature [K]")

ax0 = ax.coords[0]
ax0.set_format_unit(u.arcmin)
ax1 = ax.coords[1]
ax1.set_format_unit(u.km/u.s)

ax.set_ylabel("Velocity [km/s]")
ax.set_xlabel("Offset [arcmin]")
No description has been provided for this image

We can put all this together:

In [10]:
# we will use the peak intensity for future display
# the warning here can be ignored because the cube is small,
# but we don't silence it because it's a legit warning when dealing with big cubes
mx = cube.max(axis=0).value

WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.max at 0x7fbcc0797ce0>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
In [11]:
pl.figure(figsize=(12,6))
ax = pl.subplot(121, projection=cube.wcs.celestial)
ax.imshow(mx)
path.show_on_axis(ax, spacing=1, color='w')

ww = wcs.WCS(pvdiagram.header)

ax = pl.subplot(122, projection=ww)
im = ax.imshow(pvdiagram.data)

ax0 = ax.coords[0]
ax0.set_format_unit(u.arcmin)
ax1 = ax.coords[1]
ax1.set_format_unit(u.km/u.s)

ax.set_ylabel("Velocity [km/s]")
ax.set_xlabel("Offset [arcmin]")
No description has been provided for this image

PV Extraction from Sky Coordinates¶

We can also make paths using celestial coordinates by passing coordinates defined in an~astropy.coordinates.SkyCoord object to ~pvextractor.Path.

In [12]:
from astropy.coordinates import SkyCoord
In [13]:
skypath = Path(SkyCoord([3.4, 3.43, 3.42]*u.h, [30.5, 30.75, 30.5]*u.deg, frame='fk5'))

We can plot again; the coordinates will be automatically determined

In [14]:
ax = pl.subplot(111, projection=cube.wcs.celestial)
ax.imshow(cube[25].value)
skypath.show_on_axis(ax, spacing=1, color='w')
Out[14]:
<matplotlib.lines.Line2D at 0x7fbcbc8a8690>
No description has been provided for this image
In [15]:
pvdiagram2 = extract_pv_slice(cube=cube, path=skypath)
pvdiagram2
Out[15]:
<astropy.io.fits.hdu.image.PrimaryHDU at 0x7fbcbc8aa490>
In [16]:
pl.figure(figsize=(12,6))
ax = pl.subplot(121, projection=cube.wcs.celestial)
ax.imshow(mx)
skypath.show_on_axis(ax, spacing=1, color='w')

ww = wcs.WCS(pvdiagram2.header)

ax = pl.subplot(122, projection=ww)
im = ax.imshow(pvdiagram2.data)

ax0 = ax.coords[0]
ax0.set_format_unit(u.arcmin)
ax1 = ax.coords[1]
ax1.set_format_unit(u.km/u.s)

ax.set_ylabel("Velocity [km/s]")
ax.set_xlabel("Offset [arcmin]")
No description has been provided for this image

We can also change the aspect ratio of the PV diagram:

In [17]:
pl.figure(figsize=(12,6))
ax = pl.subplot(121, projection=cube.wcs.celestial)
ax.imshow(mx)
skypath.show_on_axis(ax, spacing=1, color='w')

ww = wcs.WCS(pvdiagram2.header)

ax = pl.subplot(122, projection=ww)
im = ax.imshow(pvdiagram2.data)
ax.set_aspect(2)

ax0 = ax.coords[0]
ax0.set_format_unit(u.arcmin)
ax1 = ax.coords[1]
ax1.set_format_unit(u.km/u.s)

ax.set_ylabel("Velocity [km/s]")
ax.set_xlabel("Offset [arcmin]")
No description has been provided for this image

PV Extraction with Spatial Averaging¶

~pvextractor.Path allows you to specify a width to average over, which specifies a spatial range around the path to average over:

In [18]:
skypath2 = Path(SkyCoord([3.4, 3.429, 3.42]*u.h,
                         [30.5, 30.75, 30.5]*u.deg,
                         frame='fk5'),
                width=2*u.arcmin)
In [19]:
pvdiagram3 = extract_pv_slice(cube=cube, path=skypath2)

























We can plot this path as a set of patches to show where we averaged.  The default spacing is 1 pixel,so we plot 1-pixel chunks.


















In [20]:






pl.figure(figsize=(12,6))
ax = pl.subplot(121, projection=cube.wcs.celestial)
ax.imshow(mx)

skypath2.show_on_axis(ax, spacing=1,
                      edgecolor='w', linestyle=':',
                      linewidth=0.75)

ww = wcs.WCS(pvdiagram3.header)

ax = pl.subplot(122, projection=ww)
im = ax.imshow(pvdiagram3.data)
ax.set_aspect(2.5)

cb = pl.colorbar(mappable=im)
cb.set_label("Brightness Temperature [K]")

ax0 = ax.coords[0]
ax0.set_format_unit(u.arcmin)
ax1 = ax.coords[1]
ax1.set_format_unit(u.km/u.s)

ax.set_ylabel("Velocity [km/s]")
ax.set_xlabel("Offset [arcmin]")
























No description has been provided for this image























We can also have more widely spaced chunks.


Note that the spacing given to extract_pv_slice affects the shape of the output PV diagram, so we also change the aspect ratio:


















In [21]:






pvdiagram4 = extract_pv_slice(cube=cube, path=skypath2, spacing=5)

pl.figure(figsize=(12,6))
ax = pl.subplot(121, projection=cube.wcs.celestial)
ax.imshow(mx)
skypath2.show_on_axis(ax, spacing=5,
                      edgecolor='w', linestyle=':',
                      linewidth=0.75)

ww = wcs.WCS(pvdiagram4.header)

ax = pl.subplot(122, projection=ww)
im = ax.imshow(pvdiagram4.data)
cb = pl.colorbar(mappable=im)
cb.set_label("Brightness Temperature [K]")
ax.set_aspect(0.5)

ax0 = ax.coords[0]
ax0.set_format_unit(u.arcmin)
ax1 = ax.coords[1]
ax1.set_format_unit(u.km/u.s)

ax.set_ylabel("Velocity [km/s]")
ax.set_xlabel("Offset [arcmin]")




































No description has been provided for this image























Saving¶
Finally, we can save the extracted PV diagram as a FITS file:


















In [22]:






pvdiagram.writeto("saved_pvdiagram.fits", overwrite=True)