Source code for pyspeckit.cubes.SpectralCube

"""
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 <git://github.com/astropy/pyregion.git>`_


The 'grunt work' is performed by the :py:mod:`cubes` module


"""
from __future__ import print_function

import time
import sys
import traceback
import numpy as np
import types
import copy
import itertools
from ..specwarnings import warn,PyspeckitWarning

import astropy
from astropy.io import fits
from astropy import log
from astropy import wcs
from astropy import units
from astropy.utils.console import ProgressBar
from six import iteritems, string_types
from functools import wraps

# import parent package
from .. import spectrum
from ..spectrum import smooth
from ..spectrum.units import (generate_xarr, SpectroscopicAxis,
                              SpectroscopicAxes)
from ..parallel_map import parallel_map
from ..spectrum import history

# import local things
from . import mapplot
from . import cubes

def not_for_cubes(func):

    @wraps(func)
    def wrapper(*args):
        warn("This operation ({0}) operates on the spectrum selected "
             "from the cube, e.g. with `set_spectrum` or `set_apspec`"
             ", it does not operate on the whole cube.", PyspeckitWarning)
        return func(*args)
    return wrapper

[docs]class Cube(spectrum.Spectrum): def __init__(self, filename=None, cube=None, xarr=None, xunit=None, errorcube=None, header=None, x0=0, y0=0, maskmap=None, **kwargs): """ 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 `~astropy.units.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 `~astropy.units.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`. """ if filename is not None: self.load_fits(filename) return else: if hasattr(cube, 'spectral_axis'): # Load from a SpectralCube instance self.cube = cube.hdu.data if (cube.unit in ('undefined', units.dimensionless_unscaled) and 'BUNIT' in cube._meta): self.unit = cube._meta['BUNIT'] else: self.unit = cube.unit log.debug("Self.unit: {0}".format(self.unit)) if xarr is None: if cube.spectral_axis.flags['OWNDATA']: xarr = SpectroscopicAxis(cube.spectral_axis, unit=cube.spectral_axis.unit, refX=cube.wcs.wcs.restfrq, refX_unit='Hz') else: xarr = SpectroscopicAxis(cube.spectral_axis.copy(), unit=cube.spectral_axis.unit, refX=cube.wcs.wcs.restfrq, refX_unit='Hz') if header is None: header = cube.header elif hasattr(cube, 'unit'): self.cube = cube.value self.unit = cube.unit else: self.cube = cube if hasattr(errorcube, 'spectral_axis'): # Load from a SpectralCube instance self.errorcube = errorcube.hdu.data elif hasattr(errorcube, 'unit'): self.errorcube = errorcube.value else: self.errorcube = errorcube if hasattr(xarr, 'flags'): log.debug("XARR flags: {0}".format(xarr.flags)) self.xarr = generate_xarr(xarr, unit=xunit) if hasattr(xarr, 'flags'): log.debug("self.xarr flags: {0}".format(xarr.flags)) self.header = header self.error = None if self.cube is not None: self.data = self.cube[:,int(y0),int(x0)] if not hasattr(self, '_unit'): self.unit = units.dimensionless_unscaled log.debug("Self.unit before header: {0}".format(self.unit)) if self.header is not None: self.parse_header(self.header) else: log.debug("self.header is None: {0}".format(self.header)) self.unit = 'undefined' self.header = fits.Header() log.debug("Self.unit after header: {0}".format(self.unit)) if maskmap is not None: if maskmap.ndim != 2: raise ValueError("Mask map must be two-dimensional.") self.maskmap = maskmap else: self.maskmap = np.ones(self.cube.shape[1:],dtype='bool') if isinstance(filename,str): self.fileprefix = filename.rsplit('.', 1)[0] # Everything prior to .fits or .txt else: self.fileprefix = "pyfitsHDU" self.plotter = spectrum.plotters.Plotter(self) self._register_fitters() self.specfit = spectrum.fitters.Specfit(self,Registry=self.Registry) self.baseline = spectrum.baseline.Baseline(self) self.speclines = spectrum.speclines # Initialize writers self.writer = {} for writer in spectrum.writers.writers: self.writer[writer] = spectrum.writers.writers[writer](self) # Special. This needs to be modified to be more flexible; for now I need it to work for nh3 self.plot_special = None self.plot_special_kwargs = {} self._modelcube = None if self.header: self.wcs = wcs.WCS(self.header) self.wcs.wcs.fix() self._spectral_axis_number = self.wcs.wcs.spec+1 self._first_cel_axis_num = np.where(self.wcs.wcs.axis_types // 1000 == 2)[0][0]+1 # TODO: Improve this!!! self.system = ('galactic' if ('CTYPE{0}'.format(self._first_cel_axis_num) in self.header and 'GLON' in self.header['CTYPE{0}'.format(self._first_cel_axis_num)]) else 'celestial') else: self._spectral_axis_number = 2 self._first_cel_axis_num = 0 self.system = 'PIXEL' self.mapplot = mapplot.MapPlotter(self)
[docs] def load_fits(self, fitsfile): try: from spectral_cube import SpectralCube except ImportError: raise ImportError("Could not import spectral_cube. As of pyspeckit" " 0.17, spectral_cube is required for cube reading. " "It can be pip installed or acquired from " "spectral-cube.rtfd.org.") mycube = SpectralCube.read(fitsfile) return self.load_spectral_cube(mycube)
[docs] def load_spectral_cube(self, cube): """ Load the cube from a spectral_cube.SpectralCube object """ self.__init__(cube=cube)
def __repr__(self): return (r'<Cube object over spectral range %6.5g :' ' %6.5g %s and flux range = [%2.1f, %2.1f]' ' %s with shape %r at %s>' % (self.xarr.min().value, self.xarr.max().value, self.xarr.unit, self.data.min(), self.data.max(), self.unit, self.cube.shape, str(hex(self.__hash__()))))
[docs] def copy(self,deep=True): """ Create a copy of the spectral cube with its own plotter, fitter, etc. Useful for, e.g., comparing smoothed to unsmoothed data """ newcube = copy.copy(self) newcube.header = copy.copy(self.header) deep_attr_lst = ['xarr', 'data', 'cube', 'maskmap', 'error', 'errorcube'] if deep: for attr in deep_attr_lst: setattr(newcube, attr, copy.copy(getattr(self, attr))) if hasattr(self, 'wcs'): newcube.wcs = self.wcs.deepcopy() newcube.header = self.header.copy() newcube.plotter = self.plotter.copy(parent=newcube) newcube._register_fitters() newcube.specfit = self.specfit.copy(parent=newcube) newcube.specfit.Spectrum.plotter = newcube.plotter newcube.baseline = self.baseline.copy(parent=newcube) newcube.baseline.Spectrum.plotter = newcube.plotter newcube.mapplot = self.mapplot.copy(parent=newcube) newcube.mapplot.Cube = newcube return newcube
def _update_header_from_xarr(self): """Uses SpectroscopiAxis' _make_header method to update Cube header""" self.header['NAXIS3'] = self.xarr.size self.xarr._make_header() sp_naxis = self._spectral_axis_number # change keywords in xarr._make_header from, e.g., CRPIX1 to CRPIX3 newhead = {(key.replace('1', str(sp_naxis)) if key.endswith('1') else key): val for key, val in iteritems(self.xarr.wcshead)} for key, val in iteritems(newhead): if isinstance(val, units.Quantity): newhead[key] = val.value elif (isinstance(val, units.CompositeUnit) or isinstance(val, units.Unit)): newhead[key] = val.to_string() log.debug("Updating header: {}: {}".format(key, val)) self.header.update(newhead)
[docs] def slice(self, start=None, stop=None, unit='pixel', preserve_fits=False, copy=True, update_header=False): """ 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' update_header : bool modifies the header of the spectral cube according to the slice """ x_in_units = self.xarr.as_unit(unit) start_ind = x_in_units.x_to_pix(start) stop_ind = x_in_units.x_to_pix(stop) if start_ind > stop_ind: start_ind, stop_ind = stop_ind, start_ind spectrum_slice = slice(start_ind,stop_ind) if not copy: raise NotImplementedError("Must copy when slicing a cube.") newcube = self.copy() newcube.cube = newcube.cube[spectrum_slice,:,:] if hasattr(newcube,'errcube'): newcube.errcube = newcube.errcube[spectrum_slice,:,:] newcube.data = newcube.data[spectrum_slice] if newcube.error is not None: newcube.error = newcube.error[spectrum_slice] newcube.xarr = newcube.xarr[spectrum_slice] # create new specfit / baseline instances (otherwise they'll be the wrong length) newcube._register_fitters() newcube.baseline = spectrum.baseline.Baseline(newcube) newcube.specfit = spectrum.fitters.Specfit(newcube,Registry=newcube.Registry) if preserve_fits: newcube.specfit.modelpars = self.specfit.modelpars newcube.specfit.parinfo = self.specfit.parinfo newcube.baseline.baselinepars = self.baseline.baselinepars newcube.baseline.order = self.baseline.order # modify the header in the new cube if update_header: newcube._update_header_from_xarr() # create a new wcs instance from the updated header newcube.wcs = wcs.WCS(newcube.header) newcube.wcs.wcs.fix() newcube._spectral_axis_number = newcube.wcs.wcs.spec + 1 newcube._first_cel_axis_num = np.where(newcube.wcs.wcs.axis_types // 1000 == 2)[0][0] + 1 return newcube
def __getitem__(self, indx): """ If [] is used on a cube, slice on the cube and use the first dimension to slice on the xarr and the data """ return Cube(xarr=self.xarr.__getitem__(indx[0]), cube=self.cube[indx], errorcube=self.errorcube[indx] if self.errorcube else None, maskmap=self.maskmap[indx[1:]] if self.maskmap is not None else None, header=self.header )
[docs] def set_spectrum(self, x, y): self.data = self.cube[:,int(y),int(x)] if self.errorcube is not None: self.error = self.errorcube[:,int(y),int(x)]
[docs] def plot_spectrum(self, x, y, plot_fit=False, **kwargs): """ Fill the .data array with a real spectrum and plot it """ self.set_spectrum(x,y) if self.plot_special is None: self.plotter(**kwargs) if plot_fit: self.plot_fit(x,y) self.plotted_spectrum = self else: sp = self.get_spectrum(x,y) sp.plot_special = types.MethodType(self.plot_special, sp) combined_kwargs = dict(kwargs.items()) combined_kwargs.update(self.plot_special_kwargs) self._spdict = sp.plot_special(**combined_kwargs) self.plotted_spectrum = sp self.plotter = sp.plotter self.plotter.refresh = lambda: [spi.plotter.refresh() for spi in self._spdict.values()] self.specfit.modelplot = [comp for spi in self._spdict.values() for comp in spi.specfit.modelplot] self.specfit._plotted_components = [comp for spi in self._spdict.values() for comp in spi.specfit._plotted_components]
[docs] def plot_fit(self, 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) """ if not hasattr(self,'parcube'): if not silent: log.info("Must run fiteach before plotting a fit. " "If you want to fit a single spectrum, " "use plot_spectrum() and specfit() directly.") return if self.plot_special is not None: # don't try to overplot a fit on a "special" plot # this is already handled in plot_spectrum return if not self.has_fit[int(y), int(x)]: # no fit to plot return self.specfit.modelpars = self.parcube[:,int(y),int(x)] if np.any(np.isnan(self.specfit.modelpars)): log.exception("Attempted to plot a model with NaN parameters.") return self.specfit.npeaks = self.specfit.fitter.npeaks self.specfit.model = self.specfit.fitter.n_modelfunc(self.specfit.modelpars, **self.specfit.fitter.modelfunc_kwargs)(self.xarr) # set the parinfo values correctly for annotations self.specfit.parinfo.values = self.parcube[:,int(y),int(x)] self.specfit.parinfo.errors = self.errcube[:,int(y),int(x)] self.specfit.fitter.parinfo.values = self.parcube[:,int(y),int(x)] self.specfit.fitter.parinfo.errors = self.errcube[:,int(y),int(x)] #for pi,p,e in zip(self.specfit.parinfo, # self.specfit.modelpars, # self.errcube[:,int(y),int(x)]): # try: # pi['value'] = p # pi['error'] = e # except ValueError: # # likely to happen for failed fits # pass self.specfit.plot_fit(**kwargs)
[docs] def plot_apspec(self, 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. """ if self.plot_special is None: self.set_apspec(aperture, coordsys=coordsys, method=method) self.plotter(reset_ylimits=reset_ylimits, **kwargs) else: #self.plot_special(reset_ylimits=reset_ylimits, **dict(kwargs.items()+self.plot_special_kwargs.items())) sp = self.get_apspec(aperture, coordsys=coordsys, wunit=wunit, method=method) sp.plot_special = types.MethodType(self.plot_special, sp) combined_kwargs = dict(kwargs.items()) combined_kwargs.update(self.plot_special_kwargs) sp.plot_special(reset_ylimits=reset_ylimits, **combined_kwargs)
[docs] def get_spectrum(self, x, y): """ Very simple: get the spectrum at coordinates x,y (inherits fitter from self) Returns a SpectroscopicAxis instance """ ct = 'CTYPE{0}'.format(self._first_cel_axis_num) header = cubes.speccen_header(fits.Header(cards=[(k,v) for k,v in iteritems(self.header) if k != 'HISTORY']), lon=x, lat=y, system=self.system, proj=(self.header[ct][-3:] if ct in self.header else 'CAR')) sp = spectrum.Spectrum(xarr=self.xarr.copy(), data=self.cube[:,int(y),int(x)], header=header, error=(self.errorcube[:,int(y),int(x)] if self.errorcube is not None else None), unit=self.unit, model_registry=self.Registry, ) sp.specfit = self.specfit.copy(parent=sp, registry=sp.Registry) # explicitly re-do this (test) sp.specfit.includemask = self.specfit.includemask.copy() sp.specfit.Spectrum = sp if hasattr(self, 'parcube'): if self.has_fit[int(y),int(x)]: # only set parameters if they're valid sp.specfit.modelpars = self.parcube[:,int(y),int(x)] if hasattr(self.specfit,'parinfo') and self.specfit.parinfo is not None: # set the parinfo values correctly for annotations for pi,p,e in zip(sp.specfit.parinfo, sp.specfit.modelpars, self.errcube[:,int(y),int(x)]): try: pi['value'] = p pi['error'] = e except ValueError: pass if hasattr(self.specfit,'fitter') and self.specfit.fitter is not None: sp.specfit.fitter.mpp = sp.specfit.modelpars # also for annotations (differs depending on which function... sigh... need to unify) sp.specfit.npeaks = self.specfit.fitter.npeaks sp.specfit.fitter.npeaks = len(sp.specfit.modelpars) / sp.specfit.fitter.npars sp.specfit.fitter.parinfo = sp.specfit.parinfo try: sp.specfit.model = sp.specfit.fitter.n_modelfunc(sp.specfit.modelpars, **sp.specfit.fitter.modelfunc_kwargs)(sp.xarr) except ValueError: # possibly invalid model parameters, just skip sp.specfit.model = np.zeros_like(sp.data) return sp
[docs] def get_apspec(self, 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 """ if coordsys is not None: wcs = self.mapplot.wcs else: wcs = None data = cubes.extract_aperture(self.cube, aperture, coordsys=coordsys, wcs=wcs, method=method, **kwargs) if self.errorcube is not None: error = cubes.extract_aperture(self.errorcube, aperture, coordsys=coordsys, wcs=self.mapplot.wcs, method='error', **kwargs) else: error = None ct = 'CTYPE{0}'.format(self._first_cel_axis_num) header = cubes.speccen_header(fits.Header(cards=[(k,v) for k,v in iteritems(self.header) if k != 'HISTORY']), lon=aperture[0], lat=aperture[1], system=self.system, proj=self.header[ct][-3:]) if len(aperture) == 3: header['APRADIUS'] = aperture[2] if len(aperture) == 5: header['APMAJ'] = aperture[2] header['APMIN'] = aperture[3] header['APREFF'] = (aperture[2]*aperture[3])**0.5 header['APPA'] = aperture[4] sp = spectrum.Spectrum(xarr=self.xarr.copy(), data=data, error=error, header=header, model_registry=self.Registry, ) sp.specfit = self.specfit.copy(parent=sp, registry=sp.Registry) return sp
[docs] def set_apspec(self, aperture, coordsys=None, method='mean'): """ Extract an aperture using cubes.extract_aperture (defaults to Cube coordinates) """ if coordsys is not None: self.data = cubes.extract_aperture(self.cube, aperture, coordsys=coordsys, wcs=self.mapplot.wcs, method=method) else: self.data = cubes.extract_aperture(self.cube, aperture, coordsys=None, method=method)
[docs] def get_modelcube(self, update=False, multicore=1): """ Return or generate a "model cube", which will have the same shape as the ``.cube`` but will have spectra generated from the fitted model. If the model cube does not yet exist, one will be generated Parameters ---------- update : bool If the cube has already been computed, set this to ``True`` to recompute the model. multicore: int if >1, try to use multiprocessing via parallel_map to run on multiple cores """ if self._modelcube is None or update: yy,xx = np.indices(self.parcube.shape[1:]) nanvals = np.any(~np.isfinite(self.parcube),axis=0) isvalid = np.any(self.parcube, axis=0) & ~nanvals valid_pixels = zip(xx[isvalid], yy[isvalid]) self._modelcube = np.full_like(self.cube, np.nan) def model_a_pixel(xy): x,y = int(xy[0]), int(xy[1]) self._modelcube[:,y,x] = self.specfit.get_full_model(pars=self.parcube[:,y,x]) return ((x,y), self._modelcube[:,y,x]) if multicore > 1: sequence = [(x,y) for x,y in valid_pixels] result = parallel_map(model_a_pixel, sequence, numcores=multicore) merged_result = [core_result for core_result in result if core_result is not None] for mr in merged_result: ((x,y), model) = mr x = int(x) y = int(y) self._modelcube[:,y,x] = model else: # progressbar doesn't work with zip; I'm therefore giving up on # "efficiency" in memory by making a list here. for xy in ProgressBar(list(valid_pixels)): model_a_pixel(xy) return self._modelcube
[docs] def fiteach(self, errspec=None, errmap=None, guesses=(), verbose=True, verbose_level=1, quiet=True, signal_cut=3, usemomentcube=None, blank_value=0, integral=False, direct_integral=False, absorption=False, use_nearest_as_guess=False, use_neighbor_as_guess=False, start_from_point=(0,0), multicore=1, position_order=None, continuum_map=None, prevalidate_guesses=False, maskmap=None, skip_failed_fits=False, **fitkwargs): """ Fit a spectrum to each valid pixel in the cube For guesses, priority is *use_nearest_as_guess*, *usemomentcube*, *guesses*, None Once you have successfully run this function, the results will be stored in the ``.parcube`` and ``.errcube`` attributes, which are each cubes of shape ``[npars, ny, nx]``, where npars is the number of fitted parameters and ``nx``, ``ny`` are the shape of the map. ``errcube`` contains the errors on the fitted parameters (1-sigma, as returned from the Levenberg-Marquardt fit's covariance matrix). You can use the attribute ``has_fit``, which is a map of shape ``[ny,nx]`` to find which pixels have been successfully fit. 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 use_neighbor_as_guess: bool Set this keyword to use the average best-fit parameters from neighboring positions with successful fits 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. position_order: ndarray[naxis=2] 2D map of region with pixel values indicating the order in which to carry out the fitting. Any type with increasing pixel values. 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. errmap: ndarray[naxis=2] or ndarray[naxis=3] A map of errors used for the individual pixels of the spectral cube. 2D errmap results in an equal weighting of each given spectrum, while a 3D array sets individual weights of each channel 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 >1, 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. prevalidate_guesses: bool An extra check before fitting is run to make sure the guesses are all within the specified limits. May be slow, so it is off by default. It also should not be necessary, since careful checking is performed before each fit. 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`. If ``None``, will use ``self.maskmap``. integral : bool If set, the integral of each spectral fit will be computed and stored in the attribute ``.integralmap`` direct_integral : bool Return the integral of the *spectrum* (as opposed to the fitted model) over a range defined by the `integration_limits` if specified or `threshold` otherwise skip_failed_fits : bool Flag to forcibly skip failed fits that fail with "unknown error". Generally, you do not want this on, but this is the 'finger-over-the-engine-light' approach that will allow these incomprehensible failures to go by and just ignore them. Keep an eye on how many of these you get: if it's just one or two out of hundreds, then maybe those are just pathological cases that can be ignored. If it's a significant fraction, you probably want to take a different approach. """ if 'multifit' in fitkwargs: warn("The multifit keyword is no longer required. All fits " "allow for multiple components.", DeprecationWarning) if not hasattr(self.mapplot,'plane'): self.mapplot.makeplane() if maskmap is None: maskmap = self.maskmap yy,xx = np.indices(self.mapplot.plane.shape) if isinstance(self.mapplot.plane, np.ma.core.MaskedArray): OK = ((~self.mapplot.plane.mask) & maskmap.astype('bool')).astype('bool') else: OK = (np.isfinite(self.mapplot.plane) & maskmap.astype('bool')).astype('bool') # NAN guesses rule out the model too if hasattr(guesses,'shape') and guesses.shape[1:] == self.cube.shape[1:]: bad = np.isnan(guesses).sum(axis=0).astype('bool') OK &= (~bad) log.info("Fitting up to {0} spectra".format(OK.sum())) if start_from_point == 'center': start_from_point = (xx.max()/2., yy.max()/2.) if hasattr(position_order,'shape') and position_order.shape == self.cube.shape[1:]: sort_distance = np.argsort(position_order.flat) else: d_from_start = ((xx-start_from_point[1])**2 + (yy-start_from_point[0])**2)**0.5 sort_distance = np.argsort(d_from_start.flat) if use_neighbor_as_guess or use_nearest_as_guess: distance = ((xx)**2 + (yy)**2)**0.5 valid_pixels = list(zip(xx.flat[sort_distance][OK.flat[sort_distance]], yy.flat[sort_distance][OK.flat[sort_distance]])) if len(valid_pixels) != len(set(valid_pixels)): raise ValueError("There are non-unique pixels in the 'valid pixel' list. " "This should not be possible and indicates a major error.") elif len(valid_pixels) == 0: raise ValueError("No valid pixels selected.") if start_from_point not in valid_pixels: raise ValueError("The starting fit position is not among the valid" " pixels. Check your selection criteria to make " "sure you have not unintentionally excluded " "this first fit pixel.") if verbose_level > 0: log.debug("Number of valid pixels: %i" % len(valid_pixels)) guesses_are_moments = (isinstance(guesses, string_types) and guesses in ('moment','moments')) if guesses_are_moments or (usemomentcube and len(guesses)): if not hasattr(self, 'momentcube') and guesses_are_moments: self.momenteach() npars = self.momentcube.shape[0] else: npars = len(guesses) if npars == 0: raise ValueError("Parameter guesses are required.") self.parcube = np.zeros((npars,)+self.mapplot.plane.shape) self.errcube = np.zeros((npars,)+self.mapplot.plane.shape) if integral: self.integralmap = np.zeros((2,)+self.mapplot.plane.shape) # newly needed as of March 27, 2012. Don't know why. if 'fittype' in fitkwargs: self.specfit.fittype = fitkwargs['fittype'] self.specfit.fitter = self.specfit.Registry.multifitters[self.specfit.fittype] # TODO: VALIDATE THAT ALL GUESSES ARE WITHIN RANGE GIVEN THE # FITKWARG LIMITS # array to store whether pixels have fits self.has_fit = np.zeros(self.mapplot.plane.shape, dtype='bool') self._counter = 0 self._tracebacks = {} t0 = time.time() def fit_a_pixel(iixy): ii,x,y = iixy sp = self.get_spectrum(x,y) # very annoying - cannot use min/max without checking type # maybe can use np.asarray here? # cannot use sp.data.mask because it can be a scalar boolean, # which does unpredictable things. if hasattr(sp.data, 'mask') and not isinstance(sp.data.mask, (bool, np.bool_)): sp.data[sp.data.mask] = np.nan sp.error[sp.data.mask] = np.nan sp.data = np.array(sp.data) sp.error = np.array(sp.error) if errspec is not None: sp.error = errspec elif errmap is not None: if self.errorcube is not None: raise ValueError("Either the 'errmap' argument or" " self.errorcube attribute should be" " specified, but not both.") if errmap.shape == self.cube.shape[1:]: sp.error = np.ones(sp.data.shape) * errmap[int(y),int(x)] elif errmap.shape == self.cube.shape: sp.error = errmap[:, int(y), int(x)] elif self.errorcube is not None: sp.error = self.errorcube[:, int(y), int(x)] else: if ii==0: # issue the warning only once (ii==0), but always issue warn("Using data std() as error. " "If signal_cut is set, this can result in " "some pixels not being fit.", PyspeckitWarning) sp.error[:] = sp.data[sp.data==sp.data].std() if sp.error is None: raise TypeError("The Spectrum's error is unset. This should " "not be possible. Please raise an Issue.") if signal_cut > 0 and not all(sp.error == 0): if continuum_map is not None: with np.errstate(divide='raise'): snr = (sp.data-continuum_map[int(y),int(x)]) / sp.error else: with np.errstate(divide='raise'): snr = sp.data / sp.error if absorption: max_sn = np.nanmax(-1*snr) else: max_sn = np.nanmax(snr) if max_sn < signal_cut: if verbose_level > 1: log.info("Skipped %4i,%4i (s/n=%0.2g)" % (x,y,max_sn)) return elif np.isnan(max_sn): if verbose_level > 1: log.info("Skipped %4i,%4i (s/n is nan; max(data)=%0.2g, min(error)=%0.2g)" % (x,y,np.nanmax(sp.data),np.nanmin(sp.error))) return if verbose_level > 2: log.info("Fitting %4i,%4i (s/n=%0.2g)" % (x,y,max_sn)) else: max_sn = None sp.specfit.Registry = self.Registry # copy over fitter registry # Do some homework for local fits # Exclude out of bounds points xpatch, ypatch = get_neighbors(x,y,self.has_fit.shape) local_fits = self.has_fit[ypatch+y,xpatch+x] if use_nearest_as_guess and self.has_fit.sum() > 0: if verbose_level > 1 and ii == 0 or verbose_level > 4: log.info("Using nearest fit as guess") rolled_distance = np.roll(np.roll(distance, x, 0), y, 1) # If there's no fit, set its distance to be unreasonably large # so it will be ignored by argmin nearest_ind = np.argmin(rolled_distance+1e10*(~self.has_fit)) nearest_x, nearest_y = xx.flat[nearest_ind],yy.flat[nearest_ind] if np.all(np.isfinite(self.parcube[:,nearest_y,nearest_x])): gg = self.parcube[:,nearest_y,nearest_x] else: log.exception("Pixel {0},{1} had a fit including a NaN: {2}" " so it will not be used as a guess for {3},{4}" .format(nearest_x, nearest_y, self.parcube[:, nearest_y, nearest_x], x, y)) gg = guesses elif use_neighbor_as_guess and np.any(local_fits): # Array is N_guess X Nvalid_nbrs so averaging over # Axis=1 is the axis of all valid neighbors gg = np.mean(self.parcube[:, (ypatch+y)[local_fits], (xpatch+x)[local_fits]], axis=1) if np.any(~np.isfinite(gg)): log.exception("Pixel {0},{1} neighbors had non-finite guess: {2}" .format(x, y, gg)) gg = guesses elif guesses_are_moments and usemomentcube is False: raise ValueError("usemomentcube must be set to True") elif guesses_are_moments or (usemomentcube and len(guesses)): if not guesses_are_moments and ii == 0: log.warn("guesses will be ignored because usemomentcube " "was set to True.", PyspeckitWarning) if verbose_level > 1 and ii == 0: log.info("Using moment cube") gg = self.momentcube[:,int(y),int(x)] elif hasattr(guesses,'shape') and guesses.shape[1:] == self.cube.shape[1:]: if verbose_level > 1 and ii == 0: log.info("Using input guess cube") gg = guesses[:,int(y),int(x)] elif isinstance(guesses, dict): if verbose_level > 1 and ii == 0: log.info("Using input guess dict") gg = guesses[(int(y),int(x))] else: if verbose_level > 1 and ii == 0: log.info("Using input guess") gg = guesses if np.all(np.isfinite(gg)): try: with np.errstate(divide='raise'): sp.specfit(guesses=gg, quiet=verbose_level<=3, verbose=verbose_level>3, **fitkwargs) self.parcube[:,int(y),int(x)] = sp.specfit.modelpars self.errcube[:,int(y),int(x)] = sp.specfit.modelerrs if np.any(~np.isfinite(sp.specfit.modelpars)): log.exception("Fit result included nan for pixel {0},{1}: " "{2}".format(x, y, sp.specfit.modelpars)) success = False # this is basically a debug statement to try to get the # code to crash here raise KeyboardInterrupt else: success = True except Exception as ex: exc_traceback = sys.exc_info()[2] self._tracebacks[(ii,x,y)] = exc_traceback log.exception("Fit number %i at %i,%i failed on error %s" % (ii,x,y, str(ex))) log.exception("Failure was in file {0} at line {1}".format( exc_traceback.tb_frame.f_code.co_filename, exc_traceback.tb_lineno,)) traceback.print_tb(exc_traceback) log.exception("Guesses were: {0}".format(str(gg))) log.exception("Fitkwargs were: {0}".format(str(fitkwargs))) success = False if isinstance(ex, KeyboardInterrupt): raise ex # keep this out of the 'try' statement if integral and success: self.integralmap[:,int(y),int(x)] = sp.specfit.integral(direct=direct_integral, return_error=True) self.has_fit[int(y),int(x)] = success else: log.exception("Fit number {0} at {1},{2} had non-finite guesses {3}" .format(ii, x, y, guesses)) self.has_fit[int(y),int(x)] = False self.parcube[:,int(y),int(x)] = blank_value self.errcube[:,int(y),int(x)] = blank_value if integral: self.integralmap[:,int(y),int(x)] = blank_value self._counter += 1 if verbose: if ii % (min(10**(3-verbose_level),1)) == 0: snmsg = " s/n=%5.1f" % (max_sn) if max_sn is not None else "" npix = len(valid_pixels) pct = 100 * (ii+1.0)/float(npix) log.info("Finished fit %6i of %6i at (%4i,%4i)%s. Elapsed time is %0.1f seconds. %%%01.f" % (ii+1, npix, x, y, snmsg, time.time()-t0, pct)) if sp.specfit.modelerrs is None: log.exception("Fit number %i at %i,%i failed with no specific error." % (ii,x,y)) if hasattr(sp.specfit, 'mpfit_status'): log.exception("mpfit status is {0}".format(sp.specfit.mpfit_status)) log.exception("The problem is that the model errors were never set, " "which implies that the fit simply failed to finish.") log.exception("The string representation of `sp.specfit.parinfo` is: {0}" .format(sp.specfit.parinfo)) log.exception("The string representation of `sp.specfit.fitter.parinfo` is: {0}" .format(sp.specfit.fitter.parinfo)) log.exception("modelpars is: {0}".format(str(sp.specfit.modelpars))) log.exception("cube modelpars are: {0}".format(str(self.parcube[:,int(y),int(x)]))) log.exception("cube modelerrs are: {0}".format(str(self.errcube[:,int(y),int(x)]))) log.exception("Guesses were: {0}".format(str(gg))) log.exception("Fitkwargs were: {0}".format(str(fitkwargs))) if skip_failed_fits: # turn the flag into a count log.exception("The fit never completed; something has gone wrong. Failed fits = {0}".format(int(skip_failed_fits))) else: raise TypeError("The fit never completed; something has gone wrong.") # blank out the errors (and possibly the values) wherever they are zero = assumed bad # this is done after the above exception to make sure we can inspect these values if blank_value != 0: self.parcube[self.parcube == 0] = blank_value self.errcube[self.parcube == 0] = blank_value if integral: return ((x,y), sp.specfit.modelpars, sp.specfit.modelerrs, self.integralmap[:,int(y),int(x)]) else: return ((x,y), sp.specfit.modelpars, sp.specfit.modelerrs) #### BEGIN TEST BLOCK #### # This test block is to make sure you don't run a 30 hour fitting # session that's just going to crash at the end. # try a first fit for exception-catching if len(start_from_point) == 2: try0 = fit_a_pixel((0,start_from_point[0], start_from_point[1])) else: try0 = fit_a_pixel((0,valid_pixels[0][0],valid_pixels[0][1])) try: len_guesses = len(self.momentcube) if (usemomentcube or guesses_are_moments) else len(guesses) assert len(try0[1]) == len_guesses == len(self.parcube) == len(self.errcube) assert len(try0[2]) == len_guesses == len(self.parcube) == len(self.errcube) except TypeError as ex: if try0 is None: raise AssertionError("The first fitted pixel did not yield a " "fit. Please try starting from a " "different pixel.") else: raise ex except AssertionError: raise AssertionError("The first pixel had the wrong fit " "parameter shape. This is probably " "a bug; please report it.") # This is a secondary test... I'm not sure it's necessary, but it # replicates what's inside the fit_a_pixel code and so should be a # useful sanity check x,y = valid_pixels[0] sp = self.get_spectrum(x,y) sp.specfit.Registry = self.Registry # copy over fitter registry # this reproduced code is needed because the functional wrapping # required for the multicore case prevents gg from being set earlier if usemomentcube or guesses_are_moments: gg = self.momentcube[:,int(y),int(x)] elif hasattr(guesses,'shape') and guesses.shape[1:] == self.cube.shape[1:]: gg = guesses[:,int(y),int(x)] else: gg = guesses # This is NOT in a try/except block because we want to raise the # exception here if an exception is going to happen sp.specfit(guesses=gg, **fitkwargs) if prevalidate_guesses: if guesses.ndim == 3: for ii,(x,y) in ProgressBar(tuple(enumerate(valid_pixels))): pinf, _ = sp.specfit.fitter._make_parinfo(parvalues=guesses[:,int(y),int(x)], **fitkwargs) sp.specfit._validate_parinfo(pinf, 'raise') else: pinf, _ = sp.specfit.fitter._make_parinfo(parvalues=guesses, **fitkwargs) sp.specfit._validate_parinfo(pinf, 'raise') #### END TEST BLOCK #### if multicore > 1: sequence = [(ii,x,y) for ii,(x,y) in tuple(enumerate(valid_pixels))] with np.errstate(divide='raise'): result = parallel_map(fit_a_pixel, sequence, numcores=multicore) self._result = result # backup - don't want to lose data in the case of a failure # a lot of ugly hacking to deal with the way parallel_map returns # its results needs TWO levels of None-filtering, because any # individual result can be None (I guess?) but apparently (and this # part I don't believe) any individual *fit* result can be None as # well (apparently the x,y pairs can also be None?) merged_result = [core_result for core_result in result if core_result is not None] # for some reason, every other time I run this code, merged_result # ends up with a different intrinsic shape. This is an attempt to # force it to maintain a sensible shape. try: if integral: ((x,y), m1, m2, intgl) = merged_result[0] else: ((x,y), m1, m2) = merged_result[0] except ValueError: if verbose > 1: log.exception("ERROR: merged_result[0] is {0} which has the" " wrong shape".format(merged_result[0])) merged_result = itertools.chain.from_iterable(merged_result) for TEMP in merged_result: if TEMP is None: # this shouldn't be possible, but it appears to happen # anyway. parallel_map is great, up to a limit that was # reached long before this level of complexity log.debug("Skipped a None entry: {0}".format(str(TEMP))) continue try: if integral: ((x,y), modelpars, modelerrs, intgl) = TEMP else: ((x,y), modelpars, modelerrs) = TEMP except TypeError: # implies that TEMP does not have the shape ((a,b),c,d) # as above, shouldn't be possible, but it happens... log.debug("Skipped a misshapen entry: {0}".format(str(TEMP))) continue if ((len(modelpars) != len(modelerrs)) or (len(modelpars) != len(self.parcube))): raise ValueError("There was a serious problem; modelpar and" " error shape don't match that of the " "parameter cubes") if ((any([x is None for x in modelpars]) or np.any(np.isnan(modelpars)) or any([x is None for x in modelerrs]) or np.any(np.isnan(modelerrs)))): self.parcube[:,int(y),int(x)] = np.nan self.errcube[:,int(y),int(x)] = np.nan self.has_fit[int(y),int(x)] = False else: self.parcube[:,int(y),int(x)] = modelpars self.errcube[:,int(y),int(x)] = modelerrs self.has_fit[int(y),int(x)] = max(modelpars) > 0 if integral: self.integralmap[:,int(y),int(x)] = intgl else: for ii,(x,y) in enumerate(valid_pixels): fit_a_pixel((ii,x,y)) # March 27, 2014: This is EXTREMELY confusing. This isn't in a loop... # make sure the fitter / fittype are set for the cube # this has to be done within the loop because skipped-over spectra # don't ever get their fittypes set self.specfit.fitter = sp.specfit.fitter self.specfit.fittype = sp.specfit.fittype self.specfit.parinfo = sp.specfit.parinfo if verbose: log.info("Finished final fit %i. " "Elapsed time was %0.1f seconds" % (len(valid_pixels), time.time()-t0)) pars_are_finite = np.all(np.isfinite(self.parcube), axis=0) # if you see one of these exceptions, please try to produce a minimum # working example and report it as a bug. # all non-finite fit parameters should be has_fit=False assert np.all(~self.has_fit[~pars_are_finite]), "Non-finite parameters found in fits"
[docs] def momenteach(self, verbose=True, verbose_level=1, multicore=1, **kwargs): """ Return a cube of the moments of each pixel Parameters ---------- multicore: int if >1, try to use multiprocessing via parallel_map to run on multiple cores """ if not hasattr(self.mapplot,'plane'): self.mapplot.makeplane() if 'vheight' not in kwargs: kwargs['vheight'] = False yy,xx = np.indices(self.mapplot.plane.shape) if isinstance(self.mapplot.plane, np.ma.core.MaskedArray): OK = (~self.mapplot.plane.mask) * self.maskmap else: OK = np.isfinite(self.mapplot.plane) * self.maskmap valid_pixels = zip(xx[OK],yy[OK]) # run the moment process to find out how many elements are in a moment _temp_moment = self.get_spectrum(yy[OK][0],xx[OK][0]).moments(**kwargs) self.momentcube = np.zeros((len(_temp_moment),)+self.mapplot.plane.shape) t0 = time.time() def moment_a_pixel(iixy): ii,x,y = iixy sp = self.get_spectrum(x,y) self.momentcube[:,int(y),int(x)] = sp.moments(**kwargs) if verbose: if ii % 10**(3-verbose_level) == 0: log.info("Finished moment %i. " "Elapsed time is %0.1f seconds" % (ii, time.time()-t0)) return ((x,y), self.momentcube[:,int(y),int(x)]) if multicore > 1: sequence = [(ii,x,y) for ii,(x,y) in tuple(enumerate(valid_pixels))] result = parallel_map(moment_a_pixel, sequence, numcores=multicore) merged_result = [core_result.tolist() for core_result in result if core_result is not None] for TEMP in merged_result: ((x,y), moments) = TEMP self.momentcube[:,int(y),int(x)] = moments else: for ii,(x,y) in enumerate(valid_pixels): moment_a_pixel((ii,x,y)) if verbose: log.info("Finished final moment %i. " "Elapsed time was %0.1f seconds" % (OK.sum(), time.time()-t0))
[docs] def show_moment(self, momentnumber, **kwargs): """ If moments have been computed, display them in the mapplot window """ if not hasattr(self,'momentcube'): raise ValueError("Compute moments first") self.mapplot.plane = self.momentcube[momentnumber,:,:].squeeze() self.mapplot(estimator=None, **kwargs)
[docs] def show_fit_param(self, 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 """ if not hasattr(self,'parcube'): raise ValueError("Compute fit parameters first") self.mapplot.plane = self.parcube[parnumber,:,:].squeeze() self.mapplot(estimator=None, **kwargs)
[docs] def load_model_fit(self, fitsfilename, npars, npeaks=1, fittype=None, _temp_fit_loc=(0,0)): """ Load a parameter + error cube into the ``.parcube`` and ``.errcube`` attributes. The models can then be examined and plotted using ``.mapplot`` as if you had run ``.fiteach``. 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. """ try: import astropy.io.fits as pyfits except ImportError: import pyfits cubefile = pyfits.open(fitsfilename,ignore_missing_end=True) cube = cubefile[0].data if cube.shape[0] != npars * npeaks * 2: raise ValueError("The cube shape is not correct. The cube has " "first dimension = {0}, but it should be {1}. " "The keyword npars = number of parameters per " "model component, and npeaks = number of " "independent peaks. You gave npars={2} and " "npeaks={3}".format(cube.shape[0], npars*npeaks*2, npars, npeaks)) # grab a spectrum and fit it however badly you want # this is just to __init__ the relevant data structures if fittype is None: if cubefile[0].header.get('FITTYPE'): fittype = cubefile[0].header.get('FITTYPE') else: raise KeyError("Must specify FITTYPE or include it in cube header.") self.parcube = cube[:npars*npeaks,:,:] self.errcube = cube[npars*npeaks:npars*npeaks*2,:,:] if np.any(np.all(self.parcube == 0, axis=(1,2))): # there are some slices where all parameters are zero, we should # ignore this when establishing whether there's a fit (some # parameters, like fortho, can be locked to zero) self.has_fit = np.all((np.isfinite(self.parcube)), axis=0) else: self.has_fit = np.all((self.parcube != 0) & (np.isfinite(self.parcube)), axis=0) nanvals = ~np.isfinite(self.parcube) nanvals_flat = np.any(nanvals, axis=0) if np.any(nanvals): warn("NaN or infinite values encountered in parameter cube.", PyspeckitWarning) # make sure params are within limits fitter = self.specfit.Registry.multifitters[fittype] guesses,throwaway = fitter._make_parinfo(npeaks=npeaks) try: x,y = _temp_fit_loc sp = self.get_spectrum(x,y) guesses.values = self.parcube[:,int(y),int(x)] sp.specfit(fittype=fittype, guesses=guesses.values) except Exception as ex1: try: OKmask = np.any(self.parcube, axis=0) & ~nanvals_flat whereOK = np.where(OKmask) x,y = whereOK[1][0],whereOK[0][0] sp = self.get_spectrum(x,y) guesses.values = self.parcube[:,int(y),int(x)] sp.specfit(fittype=fittype, guesses=guesses.values) except Exception as ex2: log.error("Fitting the pixel at location {0} failed with error: {1}. " "Re-trying at location {2} failed with error {3}. " "Try setting _temp_fit_loc to a valid pixel".format(_temp_fit_loc, ex1, (x,y), ex2)) self.specfit.fitter = sp.specfit.fitter self.specfit.fittype = sp.specfit.fittype self.specfit.parinfo = sp.specfit.parinfo
[docs] def smooth(self,factor,**kwargs): """ Smooth the spectrum by factor `factor`. Documentation from the :mod:`cubes.spectral_smooth` module: """ factor = round(factor) self.cube = cubes.spectral_smooth(self.cube,factor,**kwargs) self.xarr = self.xarr[::factor] if hasattr(self,'data'): self.data = smooth.smooth(self.data,factor,**kwargs) if len(self.xarr) != self.cube.shape[0]: raise ValueError("Convolution resulted in different X and Y array lengths. Convmode should be 'same'.") if self.errorcube is not None: self.errorcube = cubes.spectral_smooth(self.errorcube,factor,**kwargs) self._smooth_header(factor)
__doc__ += "cubes.spectral_smooth doc: \n" + cubes.spectral_smooth.__doc__ def _smooth_header(self,factor): """ Internal - correct the FITS header parameters when smoothing """ if self.header.get('CDELT3') is not None and self.header.get('CRPIX3') is not None: self.header['CDELT3'] = self.header.get('CDELT3') * float(factor) self.header['CRPIX3'] = self.header.get('CRPIX3') / float(factor) history.write_history(self.header,"SMOOTH: Smoothed and downsampled spectrum by factor %i" % (factor)) history.write_history(self.header,"SMOOTH: Changed CRPIX3 from %f to %f" % (self.header.get('CRPIX3')*float(factor),self.header.get('CRPIX3'))) history.write_history(self.header,"SMOOTH: Changed CDELT3 from %f to %f" % (self.header.get('CRPIX3')/float(factor),self.header.get('CRPIX3')))
[docs] def write_fit(self, fitcubefilename, overwrite=False): """ Write out a fit cube containing the ``.parcube`` and ``.errcube`` using the information in the fit's parinfo to set the header keywords. The ``PLANE#`` keywords will be used to indicate the content of each plane in the data cube written to the FITS file. All of the fitted parameters will be written first, followed by all of the errors on those parameters. So, for example, if you have fitted a single gaussian to each pixel, the dimensions of the saved cube will be ``[6, ny, nx]``, and they will be the amplitude, centroid, width, error on amplitude, error on centroid, and error on width, respectively. To load such a file back in for plotting purposes, see `SpectralCube.load_model_fit`. Parameters ---------- fitcubefilename: string Filename to write to overwrite: bool Overwrite file if it exists? """ try: import astropy.io.fits as pyfits except ImportError: import pyfits try: fitcubefile = pyfits.PrimaryHDU(data=np.concatenate([self.parcube,self.errcube]), header=self.header) fitcubefile.header['FITTYPE'] = self.specfit.fittype for ii,par in enumerate(self.specfit.parinfo): kw = "PLANE%i" % ii parname = par['parname'].strip('0123456789') fitcubefile.header[kw] = parname # set error parameters for jj,par in enumerate(self.specfit.parinfo): kw = "PLANE%i" % (ii+jj+1) parname = "e"+par['parname'].strip('0123456789') fitcubefile.header[kw] = parname # overwrite the WCS fitcubefile.header['CDELT3'] = 1 fitcubefile.header['CTYPE3'] = 'FITPAR' fitcubefile.header['CRVAL3'] = 0 fitcubefile.header['CRPIX3'] = 1 except AttributeError: log.exception("Make sure you run the cube fitter first.") return if astropy.version.major >= 2 or (astropy.version.major==1 and astropy.version.minor>=3): fitcubefile.writeto(fitcubefilename, overwrite=overwrite) else: fitcubefile.writeto(fitcubefilename, clobber=overwrite)
[docs] def write_cube(self): raise NotImplementedError
[docs]class CubeStack(Cube): """ The Cube equivalent of Spectra: for stitching multiple cubes with the same spatial grid but different frequencies together """ def __init__(self, cubelist, xunit='GHz', x0=0, y0=0, maskmap=None, **kwargs): """ Initialize the Cube. Accepts FITS files. x0,y0 - initial spectrum to use (defaults to lower-left corner) """ log.info("Creating Cube Stack") cubelist = list(cubelist) for ii,cube in enumerate(cubelist): if type(cube) is str: cube = Cube(cube) cubelist[ii] = cube if cube.xarr.unit != xunit: # convert all inputs to same (non-velocity) unit cube.xarr.convert_to_unit(xunit, **kwargs) self.cubelist = cubelist log.info("Concatenating data") self.xarr = SpectroscopicAxes([sp.xarr for sp in cubelist]) self.cube = np.ma.concatenate([icube.cube for icube in cubelist]) if np.any([icube.errorcube is not None for icube in cubelist]): if all([icube.errorcube is not None for icube in cubelist]): self.errorcube = np.ma.concatenate([icube.errorcube for icube in cubelist]) else: raise ValueError("Mismatched error cubes.") else: self.errorcube = None if hasattr(self.cube,'mask'): try: if self.cube.mask in (False,np.bool_(False)): # mask causes major problems internally for numpy... self.cube = np.array(self.cube) except ValueError: # this means that self.cube.mask is an array; # techically that's alright pass self._sort() self.data = self.cube[:,int(y0),int(x0)] self.error = self.errorcube[:,int(y0),int(x0)] if self.errorcube is not None else None self.header = cubelist[0].header.copy() for cube in cubelist: for key,value in cube.header.items(): if key in ['HISTORY', 'COMMENT']: continue self.header[key] = value if self.header: self.wcs = wcs.WCS(self.header) self.wcs.wcs.fix() self._spectral_axis_number = self.wcs.wcs.spec+1 self._first_cel_axis_num = np.where(self.wcs.wcs.axis_types // 1000 == 2)[0][0]+1 # TODO: Improve this!!! self.system = ('galactic' if ('CTYPE{0}'.format(self._first_cel_axis_num) in self.header and 'GLON' in self.header['CTYPE{0}'.format(self._first_cel_axis_num)]) else 'celestial') else: self._spectral_axis_number = 3 self._first_cel_axis_num = 1 self.system = 'PIXEL' self.unit = cubelist[0].unit for cube in cubelist: if cube.unit != self.unit: raise ValueError("Mismatched units " "{0} and {1}".format(cube.unit, self.unit)) self.fileprefix = cubelist[0].fileprefix # first is the best? if maskmap is not None: self.maskmap = maskmap else: self.maskmap = np.ones(self.cube.shape[1:],dtype='bool') self._register_fitters() self.plotter = spectrum.plotters.Plotter(self) self.specfit = spectrum.fitters.Specfit(self,Registry=self.Registry) self.baseline = spectrum.baseline.Baseline(self) self.speclines = spectrum.speclines # Initialize writers TO DO: DO WRITERS WORK FOR CUBES? self.writer = {} for writer in spectrum.writers.writers: self.writer[writer] = spectrum.writers.writers[writer](self) # Special. This needs to be modified to be more flexible; for now I need it to work for nh3 self.plot_special = None self.plot_special_kwargs = {} self._modelcube = None self.mapplot = mapplot.MapPlotter(self) def _sort(self): """ Sort the data in order of increasing X (could be decreasing, but must be monotonic for plotting reasons) """ indices = self.xarr.argsort() self.xarr = self.xarr[indices] self.cube = self.cube[indices,:,:] if self.errorcube is not None: self.errorcube = self.errorcube[indices,:,:]
def get_neighbors(x, y, shape): """ Find the 9 nearest neighbors, excluding self and any out of bounds points """ ysh, xsh = shape xpyp = [(ii,jj) for ii,jj in itertools.product((-1,0,1), (-1,0,1)) if (ii+x < xsh) and (ii+x >= 0) and (jj+y < ysh) and (jj+y >= 0) and not (ii==0 and jj==0)] xpatch, ypatch = zip(*xpyp) return np.array(xpatch, dtype='int'), np.array(ypatch, dtype='int') def test_get_neighbors(): xp,yp = get_neighbors(0,0,[10,10]) assert set(xp) == {0,1} assert set(yp) == {0,1} xp,yp = get_neighbors(0,1,[10,10]) assert set(xp) == {0,1} assert set(yp) == {-1,0,1} xp,yp = get_neighbors(5,6,[10,10]) assert set(xp) == {-1,0,1} assert set(yp) == {-1,0,1} xp,yp = get_neighbors(9,9,[10,10]) assert set(xp) == {0,-1} assert set(yp) == {0,-1} xp,yp = get_neighbors(9,8,[10,10]) assert set(xp) == {-1,0} assert set(yp) == {-1,0,1}