%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-68r3afwo/spectral-cube_47c7c3cfc112439b96acad3565e0d647
  Running command git clone --filter=blob:none --quiet https://github.com/radio-astro-tools/spectral-cube/ /tmp/pip-install-68r3afwo/spectral-cube_47c7c3cfc112439b96acad3565e0d647

  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 numpy as np
from astropy.utils.data import download_file
from spectral_cube import SpectralCube
from astropy import wcs
import os

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

import radio_beam

from astropy import units as u

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()

# 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()

cube = SpectralCube.read(filename, use_dask=True)

cube

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

mx = cube.max(axis=0)

pl.imshow(mx.value, origin='lower')

<matplotlib.image.AxesImage at 0x7f3e9ce5b310>

import regions

center = regions.PixCoord(1024, 1024)
ellipse = regions.EllipsePixelRegion(center, width=550, height=400, angle=45*u.deg)

ax = pl.gca()
ax.imshow(mx.value, origin='lower')
ellipse.plot(ax=ax, facecolor='none', edgecolor='red', lw=2)

<matplotlib.patches.Ellipse at 0x7f3deb9382b0>

cutout = cube.subcube_from_regions([ellipse])
cutout

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

import pvextractor

path = pvextractor.Path([(0,0), (481,481)], width=200)

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)

<matplotlib.collections.PatchCollection at 0x7f3deb9ec0d0>

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

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)

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

reproj = cutout.reproject(header)

WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.reproject at 0x7f3e71029090>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]

reproj

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

rmax = reproj.max(axis=0)

rmax.quicklook()

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]

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)

Position-Velocity Diagrams from Disks¶

Authors¶

Learning Goals¶

Keywords¶

Summary¶

Requirements¶

Second approach¶