%pip install -r https://github.com/radio-astro-tools/tutorials/raw/master/requirements.txt
%pip install --pre -U astroquery

Collecting spectral-cube (from -r https://github.com/radio-astro-tools/tutorials/raw/master/requirements.txt (line 14))
  Cloning https://github.com/radio-astro-tools/spectral-cube/ to /tmp/pip-install-omjlwbt8/spectral-cube_ab44d9ea2ac7486ca8cbcf273ce4696b
  Running command git clone --filter=blob:none --quiet https://github.com/radio-astro-tools/spectral-cube/ /tmp/pip-install-omjlwbt8/spectral-cube_ab44d9ea2ac7486ca8cbcf273ce4696b

  Resolved https://github.com/radio-astro-tools/spectral-cube/ to commit 8ae6705d7e64739cbb08904efc660dbed98874ab

  Installing build dependencies ... -

 \

 |

 /

 done

  Getting requirements to build wheel ... -

 done

  Preparing metadata (pyproject.toml) ... -

import pylab as pl
import numpy as np
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

cube = SpectralCube.read('http://www.astropy.org/astropy-data/l1448/l1448_13co.fits')
cube

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

pl.imshow(cube[25].value, origin='lower')

<matplotlib.image.AxesImage at 0x7fca50600340>

path = Path([(20,20), (40,40), (60,20)])

ax = pl.subplot(111, projection=cube.wcs.celestial)
ax.imshow(cube[25].value)
path.show_on_axis(ax, spacing=1, color='r')

<matplotlib.lines.Line2D at 0x7fca946abc10>

pvdiagram = extract_pv_slice(cube=cube, path=path, spacing=1)
pvdiagram

<astropy.io.fits.hdu.image.PrimaryHDU at 0x7fca4ffc41f0>

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

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

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

pl.figure(figsize=(12,6))
ax = pl.subplot(121, projection=cube.wcs.celestial)
ax.imshow(mx)
path.show_on_axis(ax, spacing=1, color='r')

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

from astropy.coordinates import SkyCoord

skypath = Path(SkyCoord([3.4, 3.43, 3.42]*u.h, [30.5, 30.75, 30.5]*u.deg, frame='fk5'))

ax = pl.subplot(111, projection=cube.wcs.celestial)
ax.imshow(cube[25].value)
skypath.show_on_axis(ax, spacing=1, color='r')

<matplotlib.lines.Line2D at 0x7fca4cc94640>

pvdiagram2 = extract_pv_slice(cube=cube, path=skypath)
pvdiagram2

<astropy.io.fits.hdu.image.PrimaryHDU at 0x7fca9470a290>

pl.figure(figsize=(12,6))
ax = pl.subplot(121, projection=cube.wcs.celestial)
ax.imshow(mx)
skypath.show_on_axis(ax, spacing=1, color='r')

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

pl.figure(figsize=(12,6))
ax = pl.subplot(121, projection=cube.wcs.celestial)
ax.imshow(mx)
skypath.show_on_axis(ax, spacing=1, color='r')

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

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)

pvdiagram3 = extract_pv_slice(cube=cube, path=skypath2)

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

skypath2.show_on_axis(ax, spacing=1,
                      edgecolor='r', 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]")

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

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

Extracting and Plotting Position-Velocity Diagrams¶

Authors¶

Learning Goals¶

Keywords¶

Summary¶

Header material¶

PV Extraction from Pixel Coordinates¶

PV Extraction from Sky Coordinates¶

PV Extraction with Spatial Averaging¶

Saving¶