Cubes

Pyspeckit can do a few things with spectral cubes. The most interesting is the spectral line fitting.

Cube objects have a fiteach() method that will fit each spectral line within a cube. It can be made to do this in parallel with the multicore option.

As of version 0.16, pyspeckit cubes can be read from SpectralCube objects:

>>> pcube = pyspeckit.Cube(cube=mySpectralCube)

Otherwise, they can be created from FITS cubes on disk:

>>> pcube = pyspeckit.Cube(filename="mycube.fits")

or from arrays:

>>> mycube = np.random.randn(250,50,50)
>>> myxaxis = np.linspace(-100,100,250)
>>> pcube = pyspeckit.Cube(cube=mycube, xarr=myxaxis, xunit='km/s')

The most interesting features of the Cube object are the fiteach() method, which fits a model spectrum to each element of the cube, and mapplot, which plots up various projections of the cube.

Cube.mapplot will create an interactive plot window. You can click on any pixel shown in that window and pull up a second window showing the spectrum at that pixel. If you’ve fitted the cube, the associated best-fit model will also be shown. This interactive setup can be a bit fragile, though, so please report bugs aggressively so we can weed them out!

The interactive viewer has a few button interactions described here.

Cubes

Tools to deal with spectroscopic data cubes.

Some features in Cubes require additional packages:

• smoothing - requires agpy‘s smooth and parallel_map routines
• pyregion

The ‘grunt work’ is performed by the cubes module

class pyspeckit.cubes.SpectralCube.Cube(filename=None, cube=None, xarr=None, xunit=None, errorcube=None, header=None, x0=0, y0=0, maskmap=None, **kwargs)

Bases: pyspeckit.spectrum.classes.Spectrum

A pyspeckit Cube object. Can be created from a FITS file on disk or from an array or a spectral_cube.SpectralCube object. If an array is used to insantiate the cube, the xarr keyword must be given, specifying the X-axis units

Parameters:

filename : str, optional

The name of a FITS file to open and read from. Must be 3D

cube : np.ndarray, spectral_cube.SpectralCube, or astropy.units.Quantity

The data from which to instantiate a Cube object. If it is an array or an astropy Quantity (which is an array with attached units), the X-axis must be specified. If this is given as a SpectralCube object, the X-axis and units should be handled automatically.

xarr : np.ndarray or astropy.units.Quantity, optional

The X-axis of the spectra from each cube. This actually corresponds to axis 0, or what we normally refer to as the Z-axis of the cube, but it indicates the X-axis in a plot of intensity vs wavelength. The units for this array are specified in the xunit keyword unless a Quantity is given.

xunit : str, optional

The unit of the xarr array if xarr is given as a numpy array

errorcube : np.ndarray, spectral_cube.SpectralCube, or Quantity, optional

A cube with the same shape as the input cube providing the 1-sigma error for each voxel. This can be specified more efficiently as an error map for most use cases, but that approach has not yet been implemented. However, you can pass a 2D error map to fiteach.

header : fits.Header or dict, optional

The header associated with the data. Only needed if the cube is given as an array or a quantity.

x0, y0 : int

The initial spectrum to use. The Cube object can be treated as a pyspeckit.Spectrum object, with all the associated tools (plotter, fitter) using the set_spectrum method to select a pixel from the cube to plot and fit. However, it is generally more sensible to extract individual spectra and treat them separately using the get_spectrum method, so these keywords MAY BE DEPRECATED in the future.

maskmap : np.ndarray, optional

A boolean mask map, where True implies that the data are good. This will be used for both plotting using mapplot and fitting using fiteach.

copy(deep=True)

Create a copy of the spectrum with its own plotter, fitter, etc. Useful for, e.g., comparing smoothed to unsmoothed data

crop(x1, x2, unit=None, **kwargs)

Replace the current spectrum with a subset from x1 to x2 in current units

Fixes CRPIX1 and baseline and model spectra to match cropped data spectrum

cubes = <module 'pyspeckit.cubes.cubes' from '/var/build/user_builds/pyspeckit-keflavich/checkouts/ammonia_tau_example/pyspeckit/cubes/cubes.pyc'>

downsample(dsfactor)

Downsample the spectrum (and all of its subsidiaries) without smoothing

Parameters:

dsfactor : int

Downsampling Factor

fiteach(errspec=None, errmap=None, guesses=(), verbose=True, verbose_level=1, quiet=True, signal_cut=3, usemomentcube=False, blank_value=0, integral=True, direct=False, absorption=False, use_nearest_as_guess=False, start_from_point=(0, 0), multicore=0, continuum_map=None, **fitkwargs)

Fit a spectrum to each valid pixel in the cube

For guesses, priority is use_nearest_as_guess, usemomentcube, guesses, None

Parameters:

use_nearest_as_guess: bool :

Unless the fitted point is the first, it will find the nearest other point with a successful fit and use its best-fit parameters as the guess

start_from_point: tuple(int,int) :

Either start from the center or from a point defined by a tuple. Work outward from that starting point.

guesses: tuple or ndarray[naxis=3] :

Either a tuple/list of guesses with len(guesses) = npars or a cube of guesses with shape [npars, ny, nx]. NOT TRUE, but a good idea in principle: You can also use a dictionary of the form {(y,x): [list of length npars]}, where (y,x) specifies a pixel location. If the dictionary method is used, npars must be specified and it sets the length of the first parameter axis

signal_cut: float :

Minimum signal-to-noise ratio to “cut” on (i.e., if peak in a given spectrum has s/n less than this value, ignore it)

blank_value: float :

Value to replace non-fitted locations with. A good alternative is numpy.nan

verbose: bool :

verbose_level: int :

Controls how much is output. 0,1 - only changes frequency of updates in loop 2 - print out messages when skipping pixels 3 - print out messages when fitting pixels 4 - specfit will be verbose

multicore: int :

if >0, try to use multiprocessing via parallel_map to run on multiple cores

continuum_map: np.ndarray :

Same shape as error map. Subtract this from data before estimating noise.

flux

The data in the spectrum (flux = data, for compatibility with astropy’s Spectrum1D object).

classmethod from_hdu(hdu)

Create a pyspeckit Spectrum object from an HDU

classmethod from_spectrum1d(spec1d)

Tool to load a pyspeckit Spectrum from a specutils object

(this is intended to be temporary; long-term the pyspeckit Spectrum object will inherit from a specutils Spectrum1D object)

get_apspec(aperture, coordsys=None, method='mean', **kwargs)

Extract an aperture using cubes.extract_aperture (defaults to Cube pixel coordinates)

aperture [tuple or list] (x, y, radius)

The aperture to use when extracting the data

coordsys [ ‘celestial’ | ‘galactic’ | None]

the coordinate system the aperture is specified in None indicates pixel coordinates (default)

wunit [str]

arcsec, arcmin, or degree

get_modelcube(update=False)

get_spectrum(x, y)

Very simple: get the spectrum at coordinates x,y

(inherits fitter from self)

Returns a SpectroscopicAxis instance

getlines(linetype='radio', **kwargs)

Access a registered database of spectral lines. Will add an attribute with the name linetype, which then has properties defined by the speclines module (most likely, a table and a “show” function to display the lines)

interpnans(spec)

Interpolate over NAN values, replacing them with values interpolated from their neighbors using linear interpolation.

load_fits(fitsfile)

load_model_fit(fitsfilename, npars, npeaks=1, fittype=None, _temp_fit_loc=(0, 0))

Load a parameter + error cube into the .parcube and .errcube attributes.

Parameters:

fitsfilename : str

The filename containing the parameter cube written with write_fit

npars : int

The number of parameters in the model fit for a single spectrum

npeaks : int

The number of independent peaks fit toward each spectrum

fittype : str, optional

The name of the fittype, e.g. ‘gaussian’ or ‘voigt’, from the pyspeckit fitter registry. This is optional; it should have been written to the FITS header and will be read from there if it is not specified

_temp_fit_loc : tuple (int,int)

The initial spectrum to use to generate components of the class. This should not need to be changed.

load_spectral_cube(cube)

Load the cube from a spectral_cube.SpectralCube object

measure(z=None, d=None, fluxnorm=None, miscline=None, misctol=10.0, ignore=None, derive=True, **kwargs)

Initialize the measurements class - only do this after you have run a fitter otherwise pyspeckit will be angry!

momenteach(verbose=True, verbose_level=1, multicore=0, **kwargs)

Return a cube of the moments of each pixel

Parameters:

multicore: int :

if >0, try to use multiprocessing via parallel_map to run on multiple cores

moments(unit='km/s', **kwargs)

Return the moments of the spectrum. In order to assure that the 1st and 2nd moments are meaningful, a ‘default’ unit is set. If unit is not set, will use current unit.

Documentation imported from the moments module:

Returns the gaussian parameters of a 1D distribution by calculating its moments. Depending on the input parameters, will only output a subset of the above.

Theory, from first principles (in the absence of noise): integral(gaussian) = sqrt(2*pi*sigma^2) * amp sigma = integral / amp / sqrt(2*pi)

In the presence of noise, this gets much more complicated. The noisy approach is inspired by mpfit

Parameters:

height: :

is the background level

amplitude: :

is the maximum (or minimum) of the data after background subtraction

x: :

is the first moment

width_x: :

is the second moment

estimator: function :

A function to estimate the “height” or “background level” of the data, e.g. mean or median. If masked arrays are being used, use the np.ma versions of the numpy functions

negamp: bool or None :

Force the peak negative (True), positive (False), or the sign of the peak will be “autodetected” (negamp=None)

nsigcut: float or None :

If specified, the code will attempt to estimate the noise and only use data above/below n-sigma above the noise. The noise will be estimated from the data unless the noise is specified with noise_estimate

noise_estimate: float or None :

Guess for the noise value. Only matters if nsigcut is specified.

Returns:

(height, amplitude, x, width_x) :

parse_hdf5_header(hdr)

HDF5 reader will create a hdr dictionary from HDF5 dataset attributes if they exist. This routine will convert that dict to a pyfits header instance.

parse_header(hdr, specname=None)

Parse parameters from a .fits header into required spectrum structure parameters

parse_text_header(Table)

Grab relevant parameters from a table header (xaxis type, etc)

This function should only exist for Spectrum objects created from .txt or other atpy table type objects

plot_apspec(aperture, coordsys=None, reset_ylimits=True, wunit='arcsec', method='mean', **kwargs)

Extract an aperture using cubes.extract_aperture (defaults to Cube coordinates)

Parameters:

aperture : list

A list of aperture parameters, e.g.

• For a circular aperture, len(ap)=3: + ap = [xcen,ycen,radius]
• For an elliptical aperture, len(ap)=5: + ap = [xcen,ycen,height,width,PA]

coordsys : None or str

The coordinate system of the aperture (e.g., galactic, fk5, None for pixel)

method : ‘mean’ or ‘sum’

Either average over parellel spectra or sum them.

plot_fit(x, y, silent=False, **kwargs)

If fiteach has been run, plot the best fit at the specified location

Parameters:

x : int

y : int

The x, y coordinates of the pixel (indices 2 and 1 respectively in numpy notation)

plot_spectrum(x, y, plot_fit=False, **kwargs)

Fill the .data array with a real spectrum and plot it

set_apspec(aperture, coordsys=None, method='mean')

Extract an aperture using cubes.extract_aperture (defaults to Cube coordinates)

set_spectrum(x, y)

shape

Return the data shape (a property of the Spectrum)

show_fit_param(parnumber, **kwargs)

If pars have been computed, display them in the mapplot window

Parameters:

parnumber : int

The index of the parameter in the parameter cube

show_moment(momentnumber, **kwargs)

If moments have been computed, display them in the mapplot window

slice(start=None, stop=None, unit='pixel', preserve_fits=False, copy=True)

Slice a cube along the spectral axis (equivalent to “spectral_slab” from the spectral_cube package)

Parameters:

start : numpy.float or int

start of slice

stop : numpy.float or int

stop of slice

unit : str

allowed values are any supported physical unit, ‘pixel’

smooth(smooth, **kwargs)

Smooth the spectrum by factor smooth.

Documentation from the cubes.spectral_smooth module:

stats(statrange=(), interactive=False)

Return some statistical measures in a dictionary (somewhat self-explanatory)

Parameters:

statrange : 2-element tuple

X-range over which to perform measures

interactive : bool

specify range interactively in plotter

unit

units

write(filename, type=None, **kwargs)

Write the spectrum to a file. The available file types are listed in spectrum.writers.writers

type - what type of file to write to? If not specified, will attempt to determine type from suffix

write_cube()

write_fit(fitcubefilename, clobber=False)

Write out a fit cube using the information in the fit’s parinfo to set the header keywords

Parameters:

fitcubefilename: string :

Filename to write to

clobber: bool :

Overwrite file if it exists?

class pyspeckit.cubes.SpectralCube.CubeStack(cubelist, xunit='GHz', x0=0, y0=0, maskmap=None, **kwargs)

Bases: pyspeckit.cubes.SpectralCube.Cube

The Cube equivalent of Spectra: for stitching multiple cubes with the same spatial grid but different frequencies together

Initialize the Cube. Accepts FITS files.

x0,y0 - initial spectrum to use (defaults to lower-left corner)

MapPlot

Make plots of the cube and interactively connect them to spectrum plotting. This is really an interactive component of the package; nothing in here is meant for publication-quality plots, but more for user interactive analysis.

That said, the plotter makes use of APLpy, so it is possible to make publication-quality plots.

author:Adam Ginsburg date:03/17/2011

---

class pyspeckit.cubes.mapplot.MapPlotter(Cube=None, figure=None, doplot=False, **kwargs)

Bases: object

Class to plot a spectrum

See mapplot for use documentation; this docstring is only for initialization.

Create a map figure for future plotting

circle(x1, y1, x2, y2, **kwargs)

Plot the spectrum of a circular aperture

click(event)

Record location of downclick

copy(parent=None)

Create a copy of the map plotter with blank (uninitialized) axis & figure

[ parent ]

A spectroscopic axis instance that is the parent of the specfit instance. This needs to be specified at some point, but defaults to None to prevent overwriting a previous plot.

makeplane(estimator=<function mean at 0x7f1c8aad5b90>)

Create a “plane” view of the cube, either by slicing or projecting it or by showing a slice from the best-fit model parameter cube.

Parameters:

estimator : [ function | ‘max’ | ‘int’ | FITS filename | integer | slice ]

A non-pythonic, non-duck-typed variable. If it’s a function, apply that function along the cube’s spectral axis to obtain an estimate (e.g., mean, min, max, etc.). ‘max’ will do the same thing as passing np.max ‘int’ will attempt to integrate the image (which is why I didn’t duck-type) (integrate means sum and multiply by dx) a .fits filename will be read using pyfits (so you can make your own cover figure) an integer will get the n’th slice in the parcube if it exists If it’s a slice, slice the input data cube along the Z-axis with this slice

mapplot(convention='calabretta', colorbar=True, useaplpy=True, vmin=None, vmax=None, cmap=None, plotkwargs={}, **kwargs)

Plot up a map based on an input data cube.

The map to be plotted is selected using makeplane. The estimator keyword argument is passed to that function.

The plotted map, once shown, is interactive. You can click on it with any of the three mouse buttons.

Button 1 or keyboard ‘1’:

Plot the selected pixel’s spectrum in another window. Mark the clicked pixel with an ‘x’

Button 2 or keyboard ‘o’:

Overplot a second (or third, fourth, fifth...) spectrum in the external plot window

Button 3:

Disconnect the interactive viewer

You can also click-and-drag with button 1 to average over a circular region. This same effect can be achieved by using the ‘c’ key to set the /c/enter of a circle and the ‘r’ key to set its /r/adius (i.e., hover over the center and press ‘c’, then hover some distance away and press ‘r’).

Parameters:

convention : ‘calabretta’ or ‘griesen’

The default projection to assume for Galactic data when plotting with aplpy.

colorbar : bool

Whether to show a colorbar

plotkwargs : dict, optional

A dictionary of keyword arguments to pass to aplpy.show_colorscale or matplotlib.pyplot.imshow

useaplpy : bool

Use aplpy if a FITS header is available

vmin, vmax: float or None :

Override values for the vmin/vmax values. Will be automatically determined if left as None

.. todo: :

Allow mapplot in subfigure

plot_spectrum(event, plot_fit=True)

Connects map cube to Spectrum...

refresh()

cubes.py

From agpy, contains functions to perform various transformations on data cubes and their headers.

pyspeckit.cubes.cubes.aper_world2pix(ap, wcs, coordsys='galactic', wunit='arcsec')

Converts an elliptical aperture (x,y,width,height,PA) from WCS to pixel coordinates given an input wcs (an instance of the pywcs.WCS class). Must be a 2D WCS header.

pyspeckit.cubes.cubes.baseline_cube(cube, polyorder=None, cubemask=None, splineorder=None, numcores=None, sampling=1)

Given a cube, fit a polynomial to each spectrum

Parameters:

cube: np.ndarray :

An ndarray with ndim = 3, and the first dimension is the spectral axis

polyorder: int :

Order of the polynomial to fit and subtract

cubemask: boolean ndarray :

Mask to apply to cube. Values that are True will be ignored when fitting.

numcores : None or int

Number of cores to use for parallelization. If None, will be set to the number of available cores.

pyspeckit.cubes.cubes.blfunc_generator(x=None, polyorder=None, splineorder=None, sampling=1)

Generate a function that will fit a baseline (polynomial or spline) to a data set. Either splineorder or polyorder must be set

Parameters:

x : np.ndarray or None

The X-axis of the fitted array. Will be set to np.arange(len(data)) if not specified

polyorder : None or int

The polynomial order.

splineorder : None or int

sampling : int

The sampling rate to use for the data. Can set to higher numbers to effectively downsample the data before fitting

pyspeckit.cubes.cubes.coords_in_image(fitsfile, lon, lat, system='galactic')

Determine whether the coordinates are inside the image

pyspeckit.cubes.cubes.extract_aperture(cube, ap, r_mask=False, wcs=None, coordsys='galactic', wunit='arcsec', debug=False, method='mean')

Extract an aperture from a data cube. E.g. to acquire a spectrum of an outflow that is extended.

Cube should have shape [z,y,x], e.g. cube = fits.getdata(‘datacube.fits’)

Apertures are specified in PIXEL units with an origin of 0,0 (NOT the 1,1 fits standard!) unless wcs and coordsys are specified

Parameters:

ap : list

For a circular aperture, len(ap)=3:

ap = [xcen,ycen,radius]

For an elliptical aperture, len(ap)=5:

ap = [xcen,ycen,height,width,PA]

wcs : wcs

a pywcs.WCS instance associated with the data cube

coordsys : str

the coordinate system the aperture is specified in. Options are ‘celestial’ and ‘galactic’. Default is ‘galactic’

wunit : str

units of width/height. default ‘arcsec’, options ‘arcmin’ and ‘degree’

method : str

‘mean’ or ‘sum’ (average over spectra, or sum them) or ‘error’ for sqrt(sum-of-squares / n)

pyspeckit.cubes.cubes.flatten_header(header, delete=False)

Attempt to turn an N-dimensional fits header into a 2-dimensional header Turns all CRPIX[>2] etc. into new keywords with suffix ‘A’

header must be a fits.Header instance

pyspeckit.cubes.cubes.getspec(lon, lat, rad, cube, header, r_fits=True, inherit=True, wunit='arcsec')

Given a longitude, latitude, aperture radius (arcsec), and a cube file, return a .fits file or a spectrum.

Parameters:

lon: float :

lat: float :

longitude and latitude center of a circular aperture in WCS coordinates must be in coordinate system of the file

rad: float :

radius (default degrees) of aperture

pyspeckit.cubes.cubes.getspec_reg(cubefilename, region, **kwargs)

Aperture extraction from a cube using a pyregion circle region

The region must be in the same coordinate system as the cube header

Warning

The second argument of getspec_reg requires a pyregion region list, and therefore this code depends on pyregion.

pyspeckit.cubes.cubes.integ(file, vrange, xcen=None, xwidth=None, ycen=None, ywidth=None, **kwargs)

wrapper of subimage_integ that defaults to using the full image

pyspeckit.cubes.cubes.plane_smooth(cube, cubedim=0, parallel=True, numcores=None, **kwargs)

parallel-map the smooth function

Parameters:

parallel: bool :

defaults True. Set to false if you want serial (for debug purposes?)

numcores: int :

pass to parallel_map (None = use all available)

pyspeckit.cubes.cubes.speccen_header(header, lon=None, lat=None, proj='TAN', system='celestial', spectral_axis=3, celestial_axes=[1, 2])

Turn a cube header into a spectrum header, retaining RA/Dec vals where possible (speccen is like flatten; spec-ify would be better but, specify? nah)

Assumes 3rd axis is velocity

pyspeckit.cubes.cubes.spectral_smooth(cube, smooth_factor, downsample=True, parallel=True, numcores=None, **kwargs)

Smooth the cube along the spectral direction

pyspeckit.cubes.cubes.subcube(cube, xcen, xwidth, ycen, ywidth, header=None, dvmult=False, return_HDU=False, units='pixels', widthunits='pixels')

Crops a data cube

All units assumed to be pixel units

cube has dimensions (velocity, y, x)

xwidth and ywidth are “radius” values, i.e. half the length that will be extracted

if dvmult is set, multiple the average by DV (this is useful if you set average=sum and dvmul=True to get an integrated value)

pyspeckit.cubes.cubes.subimage_integ(cube, xcen, xwidth, ycen, ywidth, vrange, header=None, average=<function mean at 0x7f1c8aad5b90>, dvmult=False, return_HDU=False, units='pixels', zunits=None)

Returns a sub-image from a data cube integrated over the specified velocity range

All units assumed to be pixel units

cube has dimensions (velocity, y, x)

xwidth and ywidth are “radius” values, i.e. half the length that will be extracted

if dvmult is set, multiply the average by DV (this is useful if you set average=sum and dvmul=True to get an integrated value)

[edit this page on github]

[edit this page on bitbucket]