Cubes

Cubes

Tools to deal with spectroscopic data cubes.

Many features in Cubes require additional packages:

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

The ‘grunt work’ is performed by the cubes module

---

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

Bases: pyspeckit.spectrum.classes.Spectrum

Initialize the Cube. Accepts files in the following formats:

• .fits

Alternatively, you can specify the xarr, cube, and header kwargs. If nothing is specified, a blank :Cube: will be generated.

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

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, units=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

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, **fitkwargs)

Fit a spectrum to each valid pixel in the cube

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

use_nearest_as_guess [ False ]

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 [ ‘center’, (x,y) ]

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

guesses [ tuple, ndarray[naxis=3] ]

Either a tuple/list of guesses with len(guesses) = npars or a cube of guesses with shape [npars, ny, nx]

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

get_apspec(aperture, coordsys=None)

Extract an aperture using cubes.extract_aperture (defaults to Cube 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)

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 their neighbors...

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

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

measure(z=None, d=None, fluxnorm=None, miscline=None, misctol=None, 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

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.

Returns (height, amplitude, x, width_x)

the gaussian parameters of a 1D distribution by calculating its moments. Depending on the input parameters, will only output a subset of the above. “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

If using masked arrays, pass estimator=np.ma.median ‘estimator’ is used to measure the background level (height)

negamp can be used to force the peak negative (True), positive (False), or it will be “autodetected” (negamp=None)

“nsigcut” - try to estimate the noise and only use data above/below nsigma

above the noise. Estimate the noise from the data unless passed as a keyword

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

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

This should be moved to the FITSSpectrum subclass when that is available

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, **kwargs)

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

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

If fiteach has been run, plot the best fit

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

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

set_apspec(aperture, coordsys=None)

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

shape()

Return the data shape

show_fit_param(parnumber, **kwargs)

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

show_moment(momentnumber, **kwargs)

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

slice(start=None, stop=None, units='pixel', preserve_fits=False)

Slicing the spectrum

start [ numpy.float or int ]

start of slice

stop [ numpy.float or int ]

stop of slice

units [ str ]

allowed values are any supported physical unit, ‘pixel’

smooth(smooth, **kwargs)

Smooth the spectrum by factor “smooth”.

sm.smooth doc:

Smooth and downsample the data array

smooth [ float ]

Number of pixels to smooth by

smoothtype [ ‘gaussian’,’hanning’, or ‘boxcar’ ]

type of smoothing kernel to use

downsample [ bool ]

Downsample the data?

downsample_factor [ int ]

Downsample by the smoothing factor, or something else?

convmode [ ‘full’,’valid’,’same’ ]

see numpy.convolve. ‘same’ returns an array of the same length as ‘data’ (assuming data is larger than the kernel)

stats(statrange=(), interactive=False)

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

range - X-range over which to perform measures interactive - specify range interactively in plotter

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

class pyspeckit.cubes.SpectralCube.CubeStack(cubelist, xunits='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 files in the following formats:

• .fits

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

write_fit(fitcubefilename, clobber=False)

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

fitcubefilename [ string ]

Filename to write to

clobber [ bool ]

Overwrite file if it exists?

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

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 0x104594758>)

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

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

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

Plot up a map based on an input data cube.

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

kwargs are passed to aplpy.show_colorscale or matplotlib.pyplot.imshow (depending on whether aplpy is installed)

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

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 = pyfits.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

INPUTS:

wcs - a pywcs.WCS instance associated with the data cube coordsys - the coordinate system the aperture is specified in.

Options are ‘celestial’ and ‘galactic’. Default is ‘galactic’

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

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

ap = [xcen,ycen,radius]

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

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

Optional inputs:

r_mask - return mask in addition to spectrum (for error checking?)

pyspeckit.cubes.cubes.flatten_header(header)

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 pyfits.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.

lon,lat - longitude and latitude center of a circular aperture in WCS coordinates

must be in coordinate system of the file

rad - 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

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

parallel - defaults True. Set to false if you want serial (for debug

purposes?)

numcores - pass to parallel_map (None = use all available)

pyspeckit.cubes.cubes.rotcrop_cube(x1, y1, x2, y2, cubename, outname, xwidth=25, ywidth=25, in_system='galactic', out_system='equatorial', clobber=True)

Crop a data cube and then rotate it with montage

pyspeckit.cubes.cubes.speccen_header(header, lon=None, lat=None)

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, **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 0x104594758>, 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, multiple the average by DV (this is useful if you set average=sum and dvmul=True to get an integrated value)