Qube

This is the base class for the full package. It defines some basic functions, such as reading in fits files from different instruments. It also gives access to some common operations you might wish to perform, such as getting a spectrum, estimating the noise and creating moment maps.

class qubefit.qube.Qube

Initate the Qube class.

Class for handeling of a variety of data cubes. This includes reading in the data cube, and preforming basic operations on the data cube. It also allows for a couple of higher level operations, such as moment generation and simple position-velocity diagrams. The cube has been primarily tested and written for (sub)-mm data cubes, but can work with a range of optical / NIR data cubes as well.

calculate_moment(moment=0, channels=None, restfreq=None, use_model=False, **kwargs)

Calculate the moments for the data.

This method will determine the moments for the data. These are the moments w.r.t. the spectral axis, as is typical. A detailed description of the different moments can be found elsewhere. In summary, the moment 0 will yield an integrated measurement of the flux density, the moment 1 will give an estimate of the velocity field and the moment 2 is an estimate of the velocity dispersion.

Parameters:
  • moment (INT, optional) – Determines which moment to return. Currently only moments up to and including 2 are defined. The default is 0.

  • channels (TUPLE or NUMPY.ARRAY, optional) – The channels to use in creating the moment If not set, the full data cube is used. The channel ranges can be given as a numpy array or 2-element tuple with the first and last channel, where the last channel is NOT included. The default is None.

  • restfreq (FLOAT, optional) – The rest frequency to use (in Hz). It is recommended that this is defined directly in the Qube instance, i.e., when it it read in or right after. However, for convenience this can be (re-)defined here. The default is None.

  • use_model (BOOLEAN, optional) – If set to true, the moment will be calculated from the model data attribute instead of the data attribute. The default is False.

  • **kwargs (VARIED, optional) – This method will take in the keywords defined in the method get_velocity. In particular the convention keyword which can change the output unit of the moments (see below)

Raises:

NotImplementedError – Will raise a NotImplementedError if the moment keyword is set to something else than 0, 1, 2.

Returns:

To output is a full instance of Qube where the data attribute containts the moment calculation. The units of the output varies and depend on the input and the convention used to get the velocities. In general the first moment will be an integrated flux density. For example if the data of the cube is in Jy/beam and the velocities are in km/s (the default), then the moment-0 will have the units of Jy * km/s / beam. The moment-1 and moment-2 will have the same unit as the velocity convention which for the default is km/s.

Return type:

Qube

calculate_sigma(ignorezero=True, fullarray=False, channels=None, use_residual=False, plot=False, plotfile='./sigma_estimate.pdf', **kwargs)

Calculate the Gaussian noise of the data.

This method will take the measurements of each channel, create a histogram and fits a Gaussian function to the histogram. It will return the Gaussian sigma of this fit as a noise estimate. This assumes that 1) the nosie is Gaussian and 2) the signal does not signficantly affect the noise property. If it does, then the data should be masked first.

Parameters:
  • ignorezero (BOOLEAN, optional) – If set, this will ignore any zero values in the array for the gaussian fitting. This is useful if masked values are set to zero. The default is True.

  • fullarray (BOOLEAN, optional) – If set, the returned array will be a numpy array with the same shape as the input data, where each point corresponds to the sigma value of that channel. The default is False.

  • channels (NUMPY.ARRAY, optional) – The list of channels to calculate the sigma value for. This can be either a list or a numpy array. When not specfied or set to None, the full range is taken. The default is None.

  • plot (BOOLEAN, optional) – If set a QA file is made of the Gaussian fits, to determine the fit of the Gaussian for each channel. The default is False.

  • plotfile (TYPE, optional) – The file in which to save the sigma QA plot. The default is ‘./sigma_estimate.pdf’.

  • doguess (BOOLEAN, optional) – If set, the code will calculate the bin size and the intial estimates for the gaussian fit. The default is True.

  • gausspar (LIST, optional) – If a list is given it is taken as the initial guesses for the Gaussian fit. This has to be set with the number of bins for the histogram, otherwise doguess has to be set to True, which will overwrite the values given here. The default is None.

  • bins (INT, optional) – The number of bins to use in the histogram. If not set, doguess needs to be set, which calculates this value.

Returns:

The Gaussian noise (sigma) for each requested channel. If full array is set, it will generate a new instance of Qube and will return the array in the ‘data’ attribute for this Qube.

Return type:

FLOAT, NUMPY.ARRAY or Qube instance

classmethod from_fits(fitsfile, extension=0, **kwargs)

Instantiate a Qube class from a fits file.

Read in a fits file. This is the primary way to load in a data cube. The fits file will be parsed for a header and data. These will be stored in the Qube class. Several other attributes will be defined (most notably the beam, shape and instrument).

Parameters:
  • cls (QUBE) – Qube Class instance.

  • fitsfile (STRING) – The name of the fits file to be read in.

  • extension (INT, optional) – The extension of the fots file to read in. The default is 0.

  • **kwargs (DICT, optional) – keyword arguments that can be directly passed into the Qube class.

Returns:

Will return an instance of the qubefit.qube.Qube class. The instance should have the data and header attribute populated, as well as several other minor data attributes.

Return type:

Qube

gaussian_moment(mom1=None, mom2=None, channels=None, use_model=False, return_amp=False, **kwargs)

Calculate the Gaussian ‘moments’ of the cube.

This method will fit a Gaussian to the spectrum of each spatial pixel. The velocity field and velocity dispersion field are then estimated from the velocity shift of the Gaussian and the width of the Gaussian, respectively. This method is often more robust in the case of low S/N and lower resolution (see the reference paper). To help provide better convergence of the fitting routine, it is recommended to supply initial guesses to the velocity field and dispersion field, which probably come from the method ‘calculate_moment’.

This method can be very slow for a large data cube in the spatial direction.

Parameters:
  • mom1 (Qube, optional) – This is a qube instance that contains the initial guess for the velocity field. If not given this estimate will be calculated from the data cube (not recommended). The default is None.

  • mom2 (Qube, optional) – This is a qube instance that contains the initial guess for the velocity dispersion field. If not given this estimate will be calculated from the data cube (not recommended). The default is None.

  • channels (NUMPY.ARRAY or TUPLE, optional) – The channels to use in creating the gaussian ‘moments’ If not set, the full data cube is used. The channel ranges can be given as a numpy array or 2-element tuple with the first and last channel, where the last channel is NOT included. It is recommended to use a large range in channels, so the zero level can be accurately estimated, which is one of the strengths of this appraoch. The default is None.

  • use_model (BOOLEAN, optional) – If set to true, the moment will be calculated from the model data attribute instead of the data attribute. The default is False.

  • **kwargs (VARIED , optional) – This method will take in the keywords defined in the method get_velocity. In particular the convention keyword which can change the output unit of the moments.

Returns:

  • mom1 (Qube) – A qube instance where the data attribute contains the velocity field estimate. Typically in km/s.

  • mom2 (Qube) – A qube instance where the data attribute contains the velocity dispersion field estimate. Typically in km/s.

get_slice(xindex=None, yindex=None, zindex=None)

Slice the data cube.

This method will slice the data cube to extract smaller cubes or individual channels. For generality, the two spatial dimensions are the x and y indices. The frequency/velocity/wavelength dimension is denoted by the zindex. This method also updates the header files with the new information.

Parameters:
  • xindex (TUPLE or NUMPY.ARRAY, optional) – The indices to use to slice in the x-direction. This can either be a tuple, which will be intepreted as a range starting at the first up to but NOT including the second value, or a numpy array of indices. The default is None.

  • yindex (TUPLE or NUMPY.ARRAY, optional) – The indices to use to slice in the y-direction. This can either be a tuple, which will be intepreted as a range starting at the first up to but NOT including the second value, or a numpy array of indices. The default is None.

  • zindex (TUPLE or NUMPY.ARRAY, optional) – The indices to use to slice in the z-direction. This can either be a tuple, which will be intepreted as a range starting at the first up to but NOT including the second value, or a numpy array of indices. The default is None.

Raises:

ValueError – Can only crop two and three dimensional data. If another dimension is selected, this will result in a ValueError.

Returns:

Return a Qube instance with the updated data and header of the sliced data set (if present, the model, sig and variance are will also be siced).

Return type:

Qube

get_spec1d(continuum_correct=False, limits=None, beam_correct=True, use_model=False, **kwargs)

Generate a 1D spectrum from the data cube.

This method will extract a spectrum from the data cube integrated over an area This is a simple sum over all of the pixels in the data cube and therefore it probably is most useful after some masking has been applied using the method ‘mask_region’. Standard is to correct for the beam to get a flux density from the region, i.e., if the data cube is in Jy/beam then the result will be a flux density in Jy. It is also possible to correct for any residual continuum.

Parameters:
  • continuum_correct (BOOLEAN, optional) – If set, this will correct the spectrum for any potential continuum flux by fitting a second order polynomial to the spectrum outside the channels given by limits. The default is False.

  • limits (LIST, optional) – Two-element list which contains the limits outside which the continuum will be fit, if set. The default is None.

  • beam_correct (BOOLEAN, optional) – Correct the flux by dividing by the area of the beam. This should be the default when dealing with interferometric data in Jy/beam. The default is True.

  • use_model (BOOLEAN, optional) – If set to true, the moment will be calculated from the model data attribute instead of the data attribute. The default is False.

  • **kwargs (VARIED , optional) – This method will take in the keywords defined in the method get_velocity. In particular the convention keyword which can change the output unit of the moments.

Returns:

Summed flux and ‘Velocity’ are returned. The first is in the same units as the data cube, whereas the latter quantity has the units as defined by the convention. This is either a velocity, frequency or wavelength array. The default is to return a velocity array in km/s.

Return type:

NUMPY.ARRAY, NUMPY.ARRAY

get_velocity(convention='radio', channels=None, as_quantity=False)

Get ‘velocities’ of the channels in a data cube.

This function will take the header information of the third dimension and the rest frequency defined also in the header, and convert the frequency values (using astropy’s units and equivalencies package) to velocities.

Parameters:
  • convention (STRING, optional) – The convention used for converting the frequencies to velocities. For possible choices see the astropy.equivalencies documentation. They are ‘radio’, ‘optical’, ‘relativistic’. In addition, the ‘frequency’ and ‘wavelength’ can be give, which will return the array in frequency (Hz) or wavelength (m). The default is ‘radio’.

  • channels (NUMPY.ARRAY, optional) – If None then get velocities for all of the channels in the data array. otherwise only the velocities of those channels that are specfied are returned. The default is None.

  • as_quantity (BOOLEAN, optional) – If set, it will return the array as a astropy.quantity instead of a unitless numpy.array. The default is False.

Raises:

ValueError – Method will raise a ValueError if the ‘CTYPE3’ in the header has a value that is unknown.

Returns:

‘Velocity’ array (in km/s or Hz, m)

Return type:

NUMPY.ARRAY or ASTROPY.QUANTITY

get_velocitywidth(**kwargs)

Calculate the channel width.

This is a small function that will get the velocities and calculate the median distance between them (i.e., width of the velocity channel). It inherits the same keywords as the get_velocity method.

mask_region(ellipse=None, rectangle=None, value=None, moment=None, channels=None, mask=None, applymask=True)

Mask a region of the data cube.

This method will mask the data cube. Several options for masking are available. Either a rectangle or ellipse for regions can be chosen or a specfiic value in the data cube. Finally, an option is to mask the region from a summed cube (moment-zero). A mask can also be read in and finally either the mask is returned as a numpy array, or the mask is applied and another Qube instance is returned. It should be noted that the masks are multiplicative.

Parameters:
  • ellipse (TUPLE, optional) – The 5-element tuple describing an ellipse (xc, yc, maj, min, ang). Each of the values in the tuple correspond to: the center of the ellipse in the x- and y-direction, (xc and yc), the length of the major and minor axis (maj and min) and the angle of the major axis (ang). The default is None.

  • rectangle (TUPLE, optional) – The 4-element tuple describing the rectangles bottom left (xb, yb) and uppper right (xt, yt) corners (xb, yb, xt, yt). The default is None.

  • value (FLOAT or LIST, optional) – Value used to mask everything below this value. This value can also be a list with the same size as the number of channels, in which case the comparison per channel is done piece-wise. The default is None.

  • moment (FLOAT, optional) – First a temporary moment image is created, using the channels keyword. Then the cube will be masked for all pixels below the value set by moment, which is given in terms of the Gaussian noise sigma. The default is None.

  • channels (TUPLE or NUMPY.ARRAY, optional) – Values used to slice the data in the spectral direction. This only needs to be set if the moment mask is requested, and if not set, it will use the full spetral range. The default is None.

  • mask (NUMPY.ARRAY, optional) – Predefined numpy array to use as a mask. The default is None.

  • applymask (BOOLEAN, optional) – This will apply the mask and return a Qube instance. If this is not set the mask will be returned as a numpy.array. The default is True.

Returns:

If applymask is true the mask data set will be returned as a Qubefit instance, if false, a numpy array of the mask will be returned.

Return type:

NUMPY.ARRAY or Qube instance

pvdiagram(PA, center, width=3.0, vshift=0.0, scale=1.0, use_model=False, **kwargs)

Create the PV diagram along the given line.

This method will take a data or model attribute from the qube instance, and create a 2D image with position along the axis and velocity on the second axis. It uses scipy’s ‘rotate’ to rotate the datacube along the axis defined by the PA and center and then extracts the total flux along the ‘slit’ with a width given by the width keyword. The final 2D output has been flipped (if needed) to have the velocity increasing upward and the extent of the pvdata is returned.

NOTE: Rotate introduces a problem in that it does not conserve the total flux of a cube when rotated. This effect is often small (<10%), but should be kept in mind when looking at these pv diagrams.

Parameters:
  • PA (FLOAT) – The position angle of the axis used to generate the PV diagram in degrees east of north.

  • center (Tuple) – Two-element tuple describing the center of the PV line, where the first value is the x position and the second value is the y position.

  • width (INT, optional) – The ‘width’ of the ‘slit’. The slit width is actually defined as two times this width plus 1 (2 x width + 1). Because of the issues with rotation, it is not advisable to set this value to 0 (i.e., a width of exactly 1 pixel). The default is 3.

  • vshift (FLOAT, optional) – This keyword allows the velocity to be shifted w.r.t. to the zero velocity as defined by the rest frequency of the cube. This is useful to make small velocity corrections.

  • scale (FLOAT, optional) – The scale to use for the pixels (x and y). This can provide a convenient way to convert the distance along the PV line from the unit pixels to more useful quantities such as kpc or arcsec. The default is 1.

  • **kwargs (VARIED , optional) – This method will take in the keywords defined in the method get_velocity. In particular the convention keyword which can change the output unit of the moments.

Returns:

The method returns a dictionary with four items defined. ‘pvdata’, which is a numpy.array that contains the PV array of values’. ‘extent’, which is a four-element tuple with the extrema of the extent of the pv array, i.e., (xmin, xmax, vmin, vmax). This can be used with the extent keyword of matplotlib.pyplot.imshow. ‘position’, which is a numpy.array of positions defining each pixel along the PV line. Finally, ‘velocity’ is a numpy.array of velocities defining each channel of the spectrum.

Return type:

DICT

to_fits(fitsfile='./cube.fits', overwrite_beam=True)

Write the qube instance to fits file.

Parameters:

fitsfile (STRING, optional) – The name of the fits file to save the qube instance to. The default is ‘./cube.fits’.

Return type:

None.