Basics of `GoFish`¶

This Notebook will walk through the very basic steps of how to load up data, set up the rotation profile and start extracting spectra.

Get the Data¶

For this we will use the CS (3-2) emission in TW Hya, originally published in Teague et al. (2018). The emission is already clearly detected, but acts as a good starting point to learn how to use GoFish. You can download the data from Harvard Dataverse, making sure to place it in the directory where this notebook is, or automatically with the cell below.

[1]:

import os
if not os.path.exists('TWHya_CS_32.fits'):
    !wget -O TWHya_CS_32.fits -q https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/LO2QZM/KYZFUQ

Load the Data¶

Loading the data is as simple as firing up GoFish and passing it the path of the fits cube to imagecube.

[2]:

import matplotlib.pyplot as plt
from gofish import imagecube
import numpy as np

[3]:

cube = imagecube('TWHya_CS_32.fits', FOV=10.0)

You can use the FOV argument here to trim down your image cube to a field of view which may speed up calculations.

Coordinate System¶

We can deproject from a pixel position, $(x_{\rm pix},\, y_{\rm pix})$ , to a disk-frame coordinate, $(r,\,\phi)$ . These transforms are done through the disk_coords function which takes geometrical properties of the disk and returns the disk coordinates.

2D Disks¶

For TW Hya we know that i=5.0 and PA=152.0 (both in degrees), so we can return the $(r,\,\phi)$ values for each pixel. Note that in GoFish, position angle is defined as the angle between North and the red-shifted disk major axis in an eastwards (anti-clockwise) direciton. We also assume the source is centered in the image such that x0=0.0 and y0=0.0.

[4]:

rvals, tvals, _ = cube.disk_coords(x0=0.0, y0=0.0, inc=5.0, PA=152.0)

[5]:

fig, ax = plt.subplots()
im = ax.imshow(rvals, origin='lower', extent=cube.extent, vmin=0.0)
cb = plt.colorbar(im, pad=0.02)
ax.set_xlabel('Offset (arcsec)')
ax.set_ylabel('Offset (arcsec)')
cb.set_label('Midplane Radius (arcsec)', rotation=270, labelpad=13)

[6]:

fig, ax = plt.subplots()
im = ax.imshow(tvals, origin='lower', extent=cube.extent)
cb = plt.colorbar(im, pad=0.02)
ax.set_xlabel('Offset (arcsec)')
ax.set_ylabel('Offset (arcsec)')
cb.set_label('Polar Angle (radians)', rotation=270, labelpad=13)

The polar angle, $\theta$ , in is radians, and runs from $-\pi$ to $\pi$ in an eastward direction with $\theta = 0$ aligning with the red-shifted major axis of the disk. Note that this is not the same as the position angle, particualrly when the disk is highly inclined.

3D Disks¶

You’ll notice that when calling disk_coords we left space for three returned parameters. The third is $z$ , the height above the midplane. By default, we assume that the disk is a razor-thin 2D disk. However, with the high spatial resolution afford by ALMA, it is now possible to resolve the vertical structure of the disk as well (see Rosenfeld et al. 2013 for a nice demonstration of this for HD 163296).

For this, we assume that the emission surface is described via,

$z(r) = z_0 \times \left( \frac{r - r_{\rm cavity}}{1^{\prime\prime}} \right)^{\psi} \times \exp \left( -\left[ \frac{r - r_{\rm cavity}}{r_{\rm taper}}\right]^{q_{\rm taper}} \right)$

where the second term acts as a ‘correction’ term in order to account for the drop in emission surface expected due to the drop in gas surface density at large radii. r_cavity can be used to describe an inner cavity, for example as in PDS 70.

As an example, we can model a conical disk with $\psi = 1$ (taking $\psi > 1$ will give a flared emission surface). We increase the inclination to make the changes due to the height more noticable. Here we use the plot_surface function to have an idea of what the surface looks like.

[7]:

fig = cube.plot_surface(inc=50.0, PA=152.0, z0=0.3, psi=1.0,
                        r_max=4.5, fill='zvals', return_fig=True)

ax = fig.axes[0]
ax.set_xlabel('Offset (arcsec)')
ax.set_ylabel('Offset (arcsec)')

[7]:

Text(0, 0.5, 'Offset (arcsec)')

Note that in the above we’ve used the fill argument to fill the background color with zvals, the emission height.

We can also assume the disk is tilted in the opposite direction by flipping the sign of the inclination.

[8]:

fig = cube.plot_surface(inc=-50.0, PA=152.0, z0=0.3, psi=1.0,
                        r_max=4.5, fill='zvals', return_fig=True)

ax = fig.axes[0]
ax.set_xlabel('Offset (arcsec)')
ax.set_ylabel('Offset (arcsec)')

[8]:

Text(0, 0.5, 'Offset (arcsec)')

Note: To comply with standard orbital conventions, $i$ should vary only between 0 $^{\circ}$ and 180 $^{\circ}$ . Thus you can include $i > 90^{\circ}$ , resulting in the correct orientation.

[9]:

fig = cube.plot_surface(inc=130.0, PA=152.0, z0=0.3, psi=1.0,
                        r_max=4.5, fill='zvals', return_fig=True)

ax = fig.axes[0]
ax.set_xlabel('Offset (arcsec)')
ax.set_ylabel('Offset (arcsec)')

[9]:

Text(0, 0.5, 'Offset (arcsec)')

Extracting a Spectrum¶

With an idea of how to specify the coordinate systems, we can use this to help us extract spectra from various locations in the disk. This is made easy with the get_spectrum function - simply provide it the coordinates where you want the spectrum and voila!

[10]:

x, y, dy = cube.get_spectrum(coords=(1.5, -1.0), frame='sky', coord_type='cartesian')

fig, ax = plt.subplots()
ax.errorbar(x, y, dy, fmt=' ', capsize=1.25, capthick=1.25, color='k', lw=1.0)
ax.step(x, y, where='mid', color='k', lw=1.0)
ax.set_xlabel('Velocity (m/s)')
ax.set_ylabel('Line Flux (Jy/beam)')
ax.set_xlim(1.84e3, 3.84e3)

[10]:

(1840.0, 3840.0)

In the above we’ve extracted a spectrum at an offset from the image center of $(1.5^{\prime\prime},\, -1.0^{\prime\prime})$ .

By default this is just a single pixel, but we can average over a larger area with the area argument. This specifies the area you want to average over as a fraction of the beam area. For example we can average over an area twice the size of the beam with area=2.0.

[11]:

x, y, dy = cube.get_spectrum(coords=(1.5, -1.0), frame='sky', coord_type='cartesian', area=2.0)

fig, ax = plt.subplots()
ax.errorbar(x, y, dy, fmt=' ', capsize=1.25, capthick=1.25, color='k', lw=1.0)
ax.step(x, y, where='mid', color='k', lw=1.0)
ax.set_xlabel('Velocity (m/s)')
ax.set_ylabel('Line Flux (Jy/beam)')
ax.set_xlim(1.84e3, 3.84e3)

[11]:

(1840.0, 3840.0)

While this does limit some of the noise, remember that you’re averaging over an area where the line profile will probably change, particularly if there’s large gradients in the source properties.

Sometimes it’s more useful to specify the location in the disk frame of reference, such as if you know where there’s an interesting feature. We can do this by swapping frame='disk'. Again we can specify the coordinates in either cartesian, coord_type='cartesian', or cylindrical frames, coord_type='cylindrical', although disk-frame coordintes are more often than not specified in cylindrical form. Remember from above that GoFish assumes that $\phi = 0$ is the red-shifted major axis and $\phi$ will increase in a clockwise direction on the sky.

Here we extract a spectrum at a $2^{\prime\prime}$ offset along the south-west minor axis.

[12]:

x, y, dy = cube.get_spectrum(coords=(2.0, np.pi / 2.0), frame='disk', coord_type='cylindrical')

fig, ax = plt.subplots()
ax.errorbar(x, y, dy, fmt=' ', capsize=1.25, capthick=1.25, color='k', lw=1.0)
ax.step(x, y, where='mid', color='k', lw=1.0)
ax.set_xlabel('Velocity (m/s)')
ax.set_ylabel('Line Flux (Jy/beam)')
ax.set_xlim(1.84e3, 3.84e3)

[12]:

(1840.0, 3840.0)

Extracing a Disk Averaged Spectrum¶

Now that we know how to deproject the data, we can use this to map a Keplerian velocity field onto the disk and calculated the projected line of sight velocity, $v_0$ , given by,

$v_0(r,\,\phi) = \sqrt{\frac{GM_{\rm star}r^2}{(r^2 + z^2)^{3/2}}} \, \cos (\phi) \, \sin (i),$

where we are expected to know $M_{\rm star}$ ( $0.88~M_{\rm sun}$ for TW Hya), and the distance to the source ( $59.5~{\rm pc}$ ).

To apply this we use the averaged_spectrum function, which is provided the geometrical properties of the disk to calculate $v_0$ , before applying these shifts to concentric annuli of the data before combining them, weighted by the area of each annulus. We must also tell it what radial range to consider through the r_min and r_max parameters, both given in arcseconds.

[13]:

x, y, dy = cube.average_spectrum(r_min=0.0, r_max=1.0, inc=5.0,
                                 PA=152., mstar=0.88, dist=59.5)

Plotting the data gives a nice spectrum:

[14]:

fig, ax = plt.subplots()
ax.errorbar(x, y, dy, fmt=' ', capsize=1.25, capthick=1.25, color='k', lw=1.0)
ax.step(x, y, where='mid', color='k', lw=1.0)
ax.set_xlabel('Velocity (m/s)')
ax.set_ylabel('Line Flux (Jy/beam)')
ax.set_xlim(1.84e3, 3.84e3)

[14]:

(1840.0, 3840.0)

Note also that by default dy is the standard deviation of the samples in the velocity bin (calculated by calling for the standard deviation of the bin with scipy.stats.binned_statistics), rather than the uncertainty on the line flux. This is because the binned data is not strictly independent owing to both small spectral and spatial correlations in the data.

When we do the spectral shifting, we actually decorrelate pixels that are close together. This is discussed in Yen et al. (2016), where we essentially gain more independent samples. In GoFish this is implemented by default and uses an numerical approach, unlike the analytical discussed in Yen et al., to account for the beam shape and orientation. This can be ignored (slightly speeding up the calculation) by setting include_spectral_decorrelation=False, as below.

[15]:

x, y, dy = cube.average_spectrum(r_min=0.0, r_max=1.0, inc=5.0,
                                 PA=152., mstar=0.88, dist=59.5,
                                 include_spectral_decorrelation=False)

fig, ax = plt.subplots()
ax.errorbar(x, y, dy, fmt=' ', capsize=1.25, capthick=1.25, color='k', lw=1.0)
ax.step(x, y, where='mid', color='k', lw=1.0)
ax.set_xlabel('Velocity (m/s)')
ax.set_ylabel('Line Flux (Jy/beam)')
ax.set_xlim(1.84e3, 3.84e3)

[15]:

(1840.0, 3840.0)

This results in a slightly larger uncertainty, but this is most noticeable when only considering area which have a large velocity gradient (most of TW Hya is a bad example for this as it is close to face on!).

We can also provide an empirical measure of the uncertainty by measuring the standard deviation in the spectrum wings. You can easily code up something which is more precise and can better define the line-free regions, but a quick work around would be to use the empirical_uncertainty=True option.

[16]:

x, y, dy = cube.average_spectrum(r_min=0.0, r_max=1.0, inc=5.0,
                                 PA=152., mstar=0.88, dist=59.5,
                                 empirical_uncertainty=True)

fig, ax = plt.subplots()
ax.errorbar(x, y, dy, fmt=' ', capsize=1.25, capthick=1.25, color='k', lw=1.0)
ax.step(x, y, where='mid', color='k', lw=1.0)
ax.set_xlabel('Velocity (m/s)')
ax.set_ylabel('Line Flux (Jy/beam)')
ax.set_xlim(1.84e3, 3.84e3)

[16]:

(1840.0, 3840.0)

Because we end up sampling the line profile at a much higher rate thanks to the Doppler shift of the lines around the azimuth (see the discussion in Teague et al. 2018), we can resample at a higher rate. This is done through the resample parameter. By default this is resample=1, so returns the velocity axis attached to the cube.

If the argument is an int, then we super-sample by that factor, while if it is a float, then that is the channel spacing. For example:

[17]:

fig, ax = plt.subplots()

# Integeter
x, y, dy = cube.average_spectrum(r_min=0.5, r_max=1.0, inc=5.0, PA=152.,
                                 mstar=0.88, dist=59.5, dr=0.1, resample=1)
ax.errorbar(x, y, dy, fmt=' ', capsize=1.25, capthick=1.25, color='k', lw=1.0)
ax.step(x, y, where='mid', color='k', lw=1.0, label='resample=1')

# Float
x, y, dy = cube.average_spectrum(r_min=0.5, r_max=1.0, inc=5.0, PA=152.,
                                 mstar=0.88, dist=59.5, dr=0.1, resample=4)
ax.errorbar(x, y, dy, fmt=' ', capsize=1.25, capthick=1.25, color='r', lw=1.0)
ax.step(x, y, where='mid', color='r', lw=1.0, label='resample=4')

ax.legend(loc=1, markerfirst=False)
ax.set_xlabel('Velocity (m/s)')
ax.set_ylabel('Flux Density (Jy/beam)')
ax.set_xlim(2.54e3, 3.14e3)

[17]:

(2540.0, 3140.0)

Similarlaly, we can use down-sample the velocity axis if we want to improve the SNR of the detection.

[18]:

fig, ax = plt.subplots()

# Integeter
x, y, dy = cube.average_spectrum(r_min=0.5, r_max=1.0, inc=5.0, PA=152.,
                                 mstar=0.88, dist=59.5, dr=0.1, resample=50.0)
ax.errorbar(x, y, dy, fmt=' ', capsize=1.25, capthick=1.25, color='k', lw=1.0)
ax.step(x, y, where='mid', color='k', lw=1.0, label='resample=50.0')

# Float
x, y, dy = cube.average_spectrum(r_min=0.5, r_max=1.0, inc=5.0, PA=152.,
                                 mstar=0.88, dist=59.5, dr=0.1, resample=200.0)
ax.errorbar(x, y, dy, fmt=' ', capsize=1.25, capthick=1.25, color='r', lw=1.0)
ax.step(x, y, where='mid', color='r', lw=1.0, label='resample=200.0')

ax.legend(loc=1, markerfirst=False)
ax.set_xlabel('Velocity (m/s)')
ax.set_ylabel('Flux Density (Jy/beam)')
ax.set_xlim(2.34e3, 3.34e3)

[18]:

(2340.0, 3340.0)

Quick Aside: By binning the data to lower spectral resolutions you are actually losing information, even though it may appear to improve things. If you know that your line profile is Gaussian, you’ll get the best results when fitting a Gaussian line when you have the highest spectral sampling. See, for example Lenz & Ayres (1992).

Using the argument unit='K', we can also return this in units of Kelvin, using the Rayleigh-Jeans approximation. And note it is exactly that, an approximation, and is a very poor conversion at sub-mm wavelengths, so use with caution! It does, however, offer a quick way to compare spectra of the same frequency, but imaged at different beam sizes.

Extracted a Disk Integrated Spectrum¶

In a very similar fashion we can extract a disk integrated spectrum, returning a spectrum in units of Jansky. Again we provide it the inner and outer radii to integrate over. Note here that this returns a single-peaked profile. This is because we’ve corrected for the Keplerian rotation of the disk prior to integrating (hence the need for specifying mstar and dist).

[19]:

fig, ax = plt.subplots()
x, y, dy = cube.integrated_spectrum(r_min=0.5, r_max=1.0, inc=5.0, PA=152.,
                                    mstar=0.88, dist=59.5, dr=0.1, resample=50.0)
ax.errorbar(x, y, dy, fmt=' ', capsize=1.25, capthick=1.25, color='k', lw=1.0)
ax.set_title('Integrated Flux = {:.0f} mJy km/s'.format(np.trapz(y, x)), fontsize=12)
ax.step(x, y, where='mid', color='k', lw=1.0)
ax.set_xlabel('Velocity (m/s)')
ax.set_ylabel('Integrated Flux (Jy)')
ax.set_xlim(1.84e3, 3.84e3)

[19]:

(1840.0, 3840.0)

If we were to remove these arguments then the function will default to a standard integration without correcting for the velocity structure of the disk. Here you can see the more typical double peaked profile of the disk, but note that the velocity integrated flux is the same.

[20]:

fig, ax = plt.subplots()
x, y, dy = cube.integrated_spectrum(r_min=0.5, r_max=1.0, inc=5.0, PA=152.)
ax.errorbar(x, y, dy, fmt=' ', capsize=1.25, capthick=1.25, color='k', lw=1.0)
ax.step(x, y, where='mid', color='k', lw=1.0)
ax.set_title('Integrated Flux = {:.0f} mJy km/s'.format(np.trapz(y, x)), fontsize=12)
ax.set_xlabel('Velocity (m/s)')
ax.set_ylabel('Integrated Flux (Jy)')
ax.set_xlim(1.84e3, 3.84e3)

[20]:

(1840.0, 3840.0)

Radial Profiles¶

For well detected sources, you can also imaging wanting to plot radial profiles of their emission. Most typically people would plot the radial profile of the total integrated intensity, showing ${\rm Jy\,beam^{-1}\,km\,s^{-1}}$ as a function of radius, or the peak flux density, ${\rm Jy\,beam^{-1}}$ , or the brightness temperature, ${\rm K}$ , as a function of radius.

Most of the time this can be readily achieved by making moment maps (shameless plug for `bettermoments <https://github.com/richteague/bettermoments>`__ to make such moment maps) and creating the radial profiles. However, even for bright sources, in the outer disk this approach is limited to what can be detected on a per pixel basis. Shifting and stacking the spectra and give a significant improvement in the outer disk.

The easiest way to get a profile is to use the radial_profile function. As with the averaged_spectrum or integrated_spectrum you want to specify the source geometry. In addition, you can (should!) provide radial bin edges or centers (not both) and the unit you want. The allowed unit values are:

'Jy m/s': Integrated spectrum in Janskys.
'K m/s': Integrated spectrum in Kelvin.
'Jy': Peak of the integrated spectrum.
'Jy/beam': Peak of the averaged spectrum.
'K': Peak of the averaged spectrum in Kelvin.

[21]:

# Calculate the peak flux density as a function of radius.
x, y, dy = cube.radial_profile(inc=5., PA=152., mstar=0.88, dist=59.5, unit='Jy/beam')

# Plot
fig, ax = plt.subplots()
ax.errorbar(x, y, dy, fmt=' ', capsize=1.25, capthick=1.25, color='k', lw=1.0)
ax.step(x, y, where='mid', color='k', lw=1.0)
ax.set_xlim(x.min(), x.max())
ax.set_xlabel('Radius (arcsec)')
ax.set_ylabel('Peak Intensity (Jy/beam)')

[21]:

Text(0, 0.5, 'Peak Intensity (Jy/beam)')

For the integrated quantities, these are integrated over the whole spectrum. You might want to apply your own clipping thresholds or integration ranges. If this is the case you can return an array of all the spectra using the radial_specrta function which takes the same values. Note this time unit can only be 'Jy', 'Jy/beam' or 'K'. First we get the deprojected spectra.

[22]:

rvals, velax, spectra, scatter = cube.radial_spectra(inc=5., PA=152., mstar=0.88, dist=59.5, unit='Jy/beam')

Here rvals is the array of bin centers used for the profile while spectra is of shape (M, 3, N) where M is the number of radial samples and N is the number of velocity samples. The second (index 1) axis is split over velocity, flux density and uncertainty.

A nice way to plot all these spectra is as follows which nicely demonstrates how the lines get weaker in the outer disk but also narrower. This is called a ‘teardrop plot’.

[23]:

cube.plot_teardrop(inc=5., PA=152., mstar=0.88, dist=59.5)

[23]:

<Axes: xlabel='Velocity (km/s)', ylabel='Radius (arcsec)'>

[ ]:

Basics of GoFish¶