import warnings import numpy as np from numpy.ma.core import nomask import dask.array as da from astropy import convolution from astropy import units as u from astropy import wcs #from astropy import log from astropy.io.fits import Header, HDUList, PrimaryHDU, BinTableHDU, FITS_rec from radio_beam import Beam, Beams from astropy.io.registry import UnifiedReadWriteMethod from . import spectral_axis from .io.core import LowerDimensionalObjectWrite from .utils import SliceWarning, BeamWarning, SmoothingWarning, FITSWarning, BeamUnitsError from . import cube_utils from . import wcs_utils from .masks import BooleanArrayMask, MaskBase from .base_class import (BaseNDClass, SpectralAxisMixinClass, SpatialCoordMixinClass, MaskableArrayMixinClass, MultiBeamMixinClass, BeamMixinClass, HeaderMixinClass ) __all__ = ['LowerDimensionalObject', 'Projection', 'Slice', 'OneDSpectrum'] class LowerDimensionalObject(u.Quantity, BaseNDClass, HeaderMixinClass): """ Generic class for 1D and 2D objects. """ def _new_view(self, obj=None, unit=None, finalize=True, **kwargs): """ kwargs are passed to _new_view of other object; only known keyword as of June 2023 is ``propagate_info`` """ # FORCE finalization to hack around https://github.com/astropy/astropy/issues/14514#issuecomment-1463935711 try: return super(LowerDimensionalObject, self)._new_view(obj=obj, unit=unit, finalize=True, **kwargs) except TypeError: return super(LowerDimensionalObject, self)._new_view(obj=obj, unit=unit, **kwargs) @property def hdu(self): if self.wcs is None: hdu = PrimaryHDU(self.value) else: hdu = PrimaryHDU(self.value, header=self.header) hdu.header['BUNIT'] = self.unit.to_string(format='fits') if self.meta is not None and 'beam' in self.meta: hdu.header.update(self.meta['beam'].to_header_keywords()) return hdu def read(self, *args, **kwargs): raise NotImplementedError() write = UnifiedReadWriteMethod(LowerDimensionalObjectWrite) def __getslice__(self, start, end, increment=None): # I don't know why this is needed, but apparently one of the inherited # classes implements getslice, which forces us to overwrite it # I can't find any examples where __getslice__ is actually implemented, # though, so this seems like a deep and frightening bug. #log.debug("Getting a slice from {0} to {1}".format(start,end)) return self.__getitem__(slice(start, end, increment)) def __getitem__(self, key, **kwargs): """ Return a new `~spectral_cube.lower_dimensional_structures.LowerDimensionalObject` of the same class while keeping other properties fixed. """ new_qty = super(LowerDimensionalObject, self).__getitem__(key) if new_qty.ndim < 2: # do not return a projection return u.Quantity(new_qty) if self._wcs is not None: if ((isinstance(key, tuple) and any(isinstance(k, slice) for k in key) and len(key) > self.ndim)): # Example cases include: indexing tricks like [:,:,None] warnings.warn("Slice {0} cannot be used on this {1}-dimensional" " array's WCS. If this is intentional, you " " should use this {2}'s ``array`` or ``quantity``" " attribute." .format(key, self.ndim, type(self)), SliceWarning ) return self.quantity[key] else: newwcs = self._wcs[key] else: newwcs = None new = self.__class__(value=new_qty.value, unit=new_qty.unit, copy=False, wcs=newwcs, meta=self._meta, mask=(self._mask[key] if self._mask is not nomask else None), header=self._header, **kwargs) new._wcs = newwcs new._meta = self._meta new._mask=(self._mask[key] if self._mask is not nomask else nomask) new._header = self._header return new def __array_finalize__(self, obj): self._wcs = getattr(obj, '_wcs', None) self._meta = getattr(obj, '_meta', None) self._mask = getattr(obj, '_mask', None) self._header = getattr(obj, '_header', None) self._spectral_unit = getattr(obj, '_spectral_unit', None) self._fill_value = getattr(obj, '_fill_value', np.nan) self._wcs_tolerance = getattr(obj, '_wcs_tolerance', 0.0) if isinstance(obj, VaryingResolutionOneDSpectrum): self._beams = getattr(obj, '_beams', None) else: self._beam = getattr(obj, '_beam', None) super(LowerDimensionalObject, self).__array_finalize__(obj) @property def __array_priority__(self): return super(LowerDimensionalObject, self).__array_priority__*2 @property def array(self): """ Get a pure array representation of the LDO. Useful when multiplying and using numpy indexing tricks. """ return np.asarray(self) @property def _data(self): # the _data property is required by several other mixins # (which probably means defining it here is a bad design) return self.array @property def quantity(self): """ Get a pure `~astropy.units.Quantity` representation of the LDO. """ return u.Quantity(self) def to(self, unit, equivalencies=[], freq=None): """ Return a new `~spectral_cube.lower_dimensional_structures.Projection` of the same class with the specified unit. See `astropy.units.Quantity.to` for further details. """ if not isinstance(unit, u.Unit): unit = u.Unit(unit) if unit == self.unit: # No copying return self if hasattr(self, 'with_spectral_unit'): freq = self.with_spectral_unit(u.Hz).spectral_axis if freq is None and 'RESTFRQ' in self.header: freq = self.header['RESTFRQ'] * u.Hz # Create the tuple of unit conversions needed. factor = cube_utils.bunit_converters(self, unit, equivalencies=equivalencies, freq=freq) converted_array = (self.quantity * factor).value # use private versions of variables, not the generated property # versions # Not entirely sure the use of __class__ here is kosher, but we do want # self.__class__, not super() new = self.__class__(value=converted_array, unit=unit, copy=True, wcs=self._wcs, meta=self._meta, mask=self._mask, header=self._header) return new @property def _mask(self): """ Annoying hack to deal with np.ma.core.is_mask failures (I don't like using __ but I think it's necessary here)""" if self.__mask is None: # need this to be *exactly* the numpy boolean False return nomask return self.__mask @_mask.setter def _mask(self, value): self.__mask = value def shrink_mask(self): """ Copy of the numpy masked_array shrink_mask method. This is essentially a hack needed for matplotlib to show images. """ m = self._mask if m.ndim and not m.any(): self._mask = nomask return self def _initial_set_mask(self, mask): """ Helper tool to validate mask when originally setting it in __new__ Note that because this is intended to be used in __new__, order matters: ``self`` must have ``_wcs``, for example. """ if mask is None: mask = BooleanArrayMask(np.ones_like(self.value, dtype=bool), self._wcs, shape=self.value.shape) elif isinstance(mask, np.ndarray): if mask.shape != self.value.shape: raise ValueError("Mask shape must match the {0} shape." .format(self.__class__.__name__) ) mask = BooleanArrayMask(mask, self._wcs, shape=self.value.shape) elif isinstance(mask, MaskBase): pass else: raise TypeError("mask of type {} is not a supported mask " "type.".format(type(mask))) # Validate the mask before setting mask._validate_wcs(new_data=self.value, new_wcs=self._wcs, wcs_tolerance=self._wcs_tolerance) self._mask = mask class Projection(LowerDimensionalObject, SpatialCoordMixinClass, MaskableArrayMixinClass, BeamMixinClass): def __new__(cls, value, unit=None, dtype=None, copy=True, wcs=None, meta=None, mask=None, header=None, beam=None, fill_value=np.nan, read_beam=False, wcs_tolerance=0.0): if np.asarray(value).ndim != 2: raise ValueError("value should be a 2-d array") if wcs is not None and wcs.wcs.naxis != 2: raise ValueError("wcs should have two dimensions") self = u.Quantity.__new__(cls, value, unit=unit, dtype=dtype, copy=copy).view(cls) self._wcs = wcs self._meta = {} if meta is None else meta self._wcs_tolerance = wcs_tolerance self._initial_set_mask(mask) self._fill_value = fill_value if header is not None: self._header = header else: self._header = Header() if beam is None: if "beam" in self.meta: beam = self.meta['beam'] elif read_beam: beam = cube_utils.try_load_beam(header) if beam is None: warnings.warn("Cannot load beam from header.", BeamWarning ) if beam is not None: self.beam = beam self.meta['beam'] = beam # TODO: Enable header updating when non-celestial slices are # properly handled in the WCS object. # self._header.update(beam.to_header_keywords()) self._cache = {} return self def with_beam(self, beam, raise_error_jybm=True): ''' Attach a new beam object to the Projection. Parameters ---------- beam : `~radio_beam.Beam` A new beam object. ''' if not isinstance(beam, Beam): raise TypeError("beam must be a radio_beam.Beam object.") self.check_jybeam_smoothing(raise_error_jybm=raise_error_jybm) meta = self.meta.copy() meta['beam'] = beam return self._new_projection_with(beam=beam, meta=meta) def with_fill_value(self, fill_value): """ Create a new :class:`Projection` or :class:`Slice` with a different ``fill_value``. """ return self._new_projection_with(fill_value=fill_value) @property def _new_thing_with(self): return self._new_projection_with def _new_projection_with(self, data=None, wcs=None, mask=None, meta=None, fill_value=None, spectral_unit=None, unit=None, header=None, wcs_tolerance=None, beam=None, **kwargs): data = self._data if data is None else data if unit is None and hasattr(data, 'unit'): if data.unit != self.unit: raise u.UnitsError("New data unit '{0}' does not" " match unit '{1}'. You can" " override this by specifying the" " `unit` keyword." .format(data.unit, self.unit)) unit = data.unit elif unit is None: unit = self.unit elif unit is not None: # convert string units to Units if not isinstance(unit, u.Unit): unit = u.Unit(unit) if hasattr(data, 'unit'): if u.Unit(unit) != data.unit: raise u.UnitsError("The specified new cube unit '{0}' " "does not match the input unit '{1}'." .format(unit, data.unit)) else: data = u.Quantity(data, unit=unit, copy=False) wcs = self._wcs if wcs is None else wcs mask = self._mask if mask is None else mask if meta is None: meta = {} meta.update(self._meta) if unit is not None: meta['BUNIT'] = unit.to_string(format='FITS') fill_value = self._fill_value if fill_value is None else fill_value if beam is None: if hasattr(self, 'beam'): beam = self.beam newproj = self.__class__(value=data, wcs=wcs, mask=mask, meta=meta, unit=unit, fill_value=fill_value, header=header or self._header, wcs_tolerance=wcs_tolerance or self._wcs_tolerance, beam=beam, **kwargs) return newproj @staticmethod def from_hdu(hdu, ext=0): ''' Return a projection from a FITS HDU. Parameters ----------- ext : int The integer index to load when given an :class:`astropy.io.fits.HDUList`. Default is 0 (the first HDU in the list. ''' if isinstance(hdu, HDUList): hdul = hdu hdu = hdul[ext] if not len(hdu.data.shape) == 2: raise ValueError("HDU must contain two-dimensional data.") meta = {} mywcs = wcs.WCS(hdu.header) if "BUNIT" in hdu.header: unit = cube_utils.convert_bunit(hdu.header["BUNIT"]) meta["BUNIT"] = hdu.header["BUNIT"] else: unit = None beam = cube_utils.try_load_beam(hdu.header) self = Projection(hdu.data, unit=unit, wcs=mywcs, meta=meta, header=hdu.header, beam=beam) return self def quicklook(self, filename=None, use_aplpy=True, aplpy_kwargs={}): """ Use `APLpy `_ to make a quick-look image of the projection. This will make the ``FITSFigure`` attribute available. If there are unmatched celestial axes, this will instead show an image without axis labels. Parameters ---------- filename : str or Non Optional - the filename to save the quicklook to. """ if use_aplpy: try: if not hasattr(self, 'FITSFigure'): import aplpy self.FITSFigure = aplpy.FITSFigure(self.hdu, **aplpy_kwargs) self.FITSFigure.show_grayscale() self.FITSFigure.add_colorbar() if filename is not None: self.FITSFigure.save(filename) except (wcs.InconsistentAxisTypesError, ImportError): self._quicklook_mpl(filename=filename) else: self._quicklook_mpl(filename=filename) def _quicklook_mpl(self, filename=None): from matplotlib import pyplot self.figure = pyplot.gcf() self.image = pyplot.imshow(self.value) if filename is not None: self.figure.savefig(filename) def convolve_to(self, beam, convolve=convolution.convolve_fft, **kwargs): """ Convolve the image to a specified beam. Parameters ---------- beam : `radio_beam.Beam` The beam to convolve to convolve : function The astropy convolution function to use, either `astropy.convolution.convolve` or `astropy.convolution.convolve_fft` Returns ------- proj : `Projection` A Projection convolved to the given ``beam`` object. """ self._raise_wcs_no_celestial() if not hasattr(self, 'beam'): raise ValueError("No beam is contained in Projection.meta.") # Check if the beams are the same. if beam == self.beam: warnings.warn("The given beam is identical to the current beam. " "Skipping convolution.") return self pixscale = wcs.utils.proj_plane_pixel_area(self.wcs.celestial)**0.5 * u.deg convolution_kernel = \ beam.deconvolve(self.beam).as_kernel(pixscale) newdata = convolve(self.value, convolution_kernel, normalize_kernel=True, **kwargs) self = Projection(newdata, unit=self.unit, wcs=self.wcs, meta=self.meta, header=self.header, beam=beam) return self def reproject(self, header, order='bilinear'): """ Reproject the image into a new header. Parameters ---------- header : `astropy.io.fits.Header` A header specifying a cube in valid WCS order : int or str, optional The order of the interpolation (if ``mode`` is set to ``'interpolation'``). This can be either one of the following strings: * 'nearest-neighbor' * 'bilinear' * 'biquadratic' * 'bicubic' or an integer. A value of ``0`` indicates nearest neighbor interpolation. """ self._raise_wcs_no_celestial() try: from reproject.version import version except ImportError: raise ImportError("Requires the reproject package to be" " installed.") # Need version > 0.2 to work with cubes from packaging.version import Version, parse if parse(version) < Version("0.3"): raise Warning("Requires version >=0.3 of reproject. The current " "version is: {}".format(version)) from reproject import reproject_interp # TODO: Find the minimal footprint that contains the header and only reproject that # (see FITS_tools.regrid_cube for a guide on how to do this) newwcs = wcs.WCS(header) shape_out = [header['NAXIS{0}'.format(i + 1)] for i in range(header['NAXIS'])][::-1] newproj, newproj_valid = reproject_interp((self.value, self.header), newwcs, shape_out=shape_out, order=order) self = Projection(newproj, unit=self.unit, wcs=newwcs, meta=self.meta, header=header, read_beam=True) return self def subimage(self, xlo='min', xhi='max', ylo='min', yhi='max'): """ Extract a region spatially. When spatial WCS dimensions are given as an `~astropy.units.Quantity`, the spatial coordinates of the 'lo' and 'hi' corners are solved together. This minimizes WCS variations due to the sky curvature when slicing from a large (>1 deg) image. Parameters ---------- [xy]lo/[xy]hi : int or `astropy.units.Quantity` or `min`/`max` The endpoints to extract. If given as a quantity, will be interpreted as World coordinates. If given as a string or int, will be interpreted as pixel coordinates. """ self._raise_wcs_no_celestial() # Solve for the spatial pixel indices together limit_dict = wcs_utils.find_spatial_pixel_index(self, xlo, xhi, ylo, yhi) slices = [slice(limit_dict[xx + 'lo'], limit_dict[xx + 'hi']) for xx in 'yx'] return self[tuple(slices)] def to(self, unit, equivalencies=[], freq=None): """ Return a new `~spectral_cube.lower_dimensional_structures.Projection` of the same class with the specified unit. See `astropy.units.Quantity.to` for further details. """ return super(Projection, self).to(unit, equivalencies, freq) # A slice is just like a projection in every way class Slice(Projection): pass class BaseOneDSpectrum(LowerDimensionalObject, MaskableArrayMixinClass, SpectralAxisMixinClass): """ Properties shared between OneDSpectrum and VaryingResolutionOneDSpectrum. """ def __new__(cls, value, unit=None, dtype=None, copy=True, wcs=None, meta=None, mask=None, header=None, spectral_unit=None, fill_value=np.nan, wcs_tolerance=0.0): #log.debug("Creating a OneDSpectrum with __new__") if np.asarray(value).ndim != 1: raise ValueError("value should be a 1-d array") if wcs is not None and wcs.wcs.naxis != 1: raise ValueError("wcs should have one dimension") self = u.Quantity.__new__(cls, value, unit=unit, dtype=dtype, copy=copy).view(cls) self._wcs = wcs self._meta = {} if meta is None else meta self._wcs_tolerance = wcs_tolerance self._initial_set_mask(mask) self._fill_value = fill_value if header is not None: self._header = header else: self._header = Header() self._spectral_unit = spectral_unit if spectral_unit is None: if 'CUNIT1' in self._header: self._spectral_unit = u.Unit(self._header['CUNIT1']) elif self._wcs is not None: self._spectral_unit = u.Unit(self._wcs.wcs.cunit[0]) return self def __repr__(self): prefixstr = '<' + self.__class__.__name__ + ' ' arrstr = np.array2string(self.filled_data[:].value, separator=',', prefix=prefixstr) return '{0}{1}{2:s}>'.format(prefixstr, arrstr, self._unitstr) @staticmethod def from_hdu(hdu, ext=0): ''' Return a OneDSpectrum from a FITS HDU or HDU list. Parameters ----------- ext : int The integer index to load when given an :class:`astropy.io.fits.HDUList`. Default is 0 (the first HDU in the list. ''' if isinstance(hdu, HDUList): hdul = hdu hdu = hdul[ext] else: hdul = HDUList([hdu]) if not len(hdu.data.shape) == 1: raise ValueError("HDU must contain one-dimensional data.") meta = {} mywcs = wcs.WCS(hdu.header) if "BUNIT" in hdu.header: unit = cube_utils.convert_bunit(hdu.header["BUNIT"]) meta["BUNIT"] = hdu.header["BUNIT"] else: unit = None with warnings.catch_warnings(): warnings.filterwarnings('ignore', category=FITSWarning) beam = cube_utils.try_load_beams(hdul) try: beams = beam _ = len(beams) except TypeError: # beam is scalar and has no len() beams = None if beams is not None: self = VaryingResolutionOneDSpectrum(hdu.data, unit=unit, wcs=mywcs, meta=meta, header=hdu.header, beams=beams) else: beam = cube_utils.try_load_beam(hdu.header) self = OneDSpectrum(hdu.data, unit=unit, wcs=mywcs, meta=meta, header=hdu.header, beam=beam) return self @property def header(self): header = super(BaseOneDSpectrum, self).header # Preserve the spectrum's spectral units if 'CUNIT1' in header and self._spectral_unit != u.Unit(header['CUNIT1']): spectral_scale = spectral_axis.wcs_unit_scale(self._spectral_unit) header['CDELT1'] *= spectral_scale header['CRVAL1'] *= spectral_scale header['CUNIT1'] = self.spectral_axis.unit.to_string(format='FITS') return header @property def spectral_axis(self): """ A `~astropy.units.Quantity` array containing the central values of each channel along the spectral axis. """ if self._wcs is None: spec_axis = np.arange(self.size) * u.one else: spec_axis = self.wcs.wcs_pix2world(np.arange(self.size), 0)[0] * \ u.Unit(self.wcs.wcs.cunit[0]) if self._spectral_unit is not None: spec_axis = spec_axis.to(self._spectral_unit) return spec_axis def quicklook(self, filename=None, drawstyle='steps-mid', **kwargs): """ Plot the spectrum with current spectral units in the currently open figure kwargs are passed to `matplotlib.pyplot.plot` Parameters ---------- filename : str or Non Optional - the filename to save the quicklook to. """ from matplotlib import pyplot ax = pyplot.gca() ax.plot(self.spectral_axis, self.filled_data[:].value, drawstyle=drawstyle, **kwargs) ax.set_xlabel(self.spectral_axis.unit.to_string(format='latex')) ax.set_ylabel(self.unit) if filename is not None: pyplot.gcf().savefig(filename) def with_spectral_unit(self, unit, velocity_convention=None, rest_value=None): newwcs, newmeta = self._new_spectral_wcs(unit, velocity_convention=velocity_convention, rest_value=rest_value) newheader = self._nowcs_header.copy() newheader.update(newwcs.to_header()) wcs_cunit = u.Unit(newheader['CUNIT1']) newheader['CUNIT1'] = unit.to_string(format='FITS') newheader['CDELT1'] *= wcs_cunit.to(unit) if self._mask is not None: newmask = self._mask.with_spectral_unit(unit, velocity_convention=velocity_convention, rest_value=rest_value) newmask._wcs = newwcs else: newmask = None return self._new_spectrum_with(wcs=newwcs, spectral_unit=unit, mask=newmask, meta=newmeta, header=newheader) def __getitem__(self, key, **kwargs): # Ideally, this could just be in VaryingResolutionOneDSpectrum, # but it's about the code is about the same length by just # keeping it here. try: kwargs['beams'] = self.beams[key] except (AttributeError, TypeError): pass new_qty = super(BaseOneDSpectrum, self).__getitem__(key) if isinstance(key, slice): new = self.__class__(value=new_qty.value, unit=new_qty.unit, copy=False, wcs=wcs_utils.slice_wcs(self._wcs, key, shape=self.shape), meta=self._meta, mask=(self._mask[key] if self._mask is not nomask else nomask), header=self._header, wcs_tolerance=self._wcs_tolerance, fill_value=self.fill_value, **kwargs) return new else: if self._mask is not nomask: # Kind of a hack; this is probably inefficient bad = self._mask.exclude()[key] if isinstance(bad, da.Array): bad = bad.compute() new_qty[bad] = np.nan return new_qty def __getattribute__(self, attrname): # This is a hack to handle dimensionality-reducing functions # We want spectrum.max() to return a Quantity, not a spectrum # Long-term, we really want `OneDSpectrum` to not inherit from # `Quantity`, but for now this approach works.... we just have # to add more functions to this list. if attrname in ('min', 'max', 'std', 'mean', 'sum', 'cumsum', 'var'): return getattr(self.quantity, attrname) else: return super(BaseOneDSpectrum, self).__getattribute__(attrname) def spectral_interpolate(self, spectral_grid, suppress_smooth_warning=False, fill_value=None): """ Resample the spectrum onto a specific grid Parameters ---------- spectral_grid : array An array of the spectral positions to regrid onto suppress_smooth_warning : bool If disabled, a warning will be raised when interpolating onto a grid that does not nyquist sample the existing grid. Disable this if you have already appropriately smoothed the data. fill_value : float Value for extrapolated spectral values that lie outside of the spectral range defined in the original data. The default is to use the nearest spectral channel in the cube. Returns ------- spectrum : OneDSpectrum """ assert spectral_grid.ndim == 1 inaxis = self.spectral_axis.to(spectral_grid.unit) indiff = np.mean(np.diff(inaxis)) outdiff = np.mean(np.diff(spectral_grid)) # account for reversed axes if outdiff < 0: spectral_grid = spectral_grid[::-1] outdiff = np.mean(np.diff(spectral_grid)) outslice = slice(None, None, -1) else: outslice = slice(None, None, 1) specslice = slice(None) if indiff >= 0 else slice(None, None, -1) inaxis = inaxis[specslice] indiff = np.mean(np.diff(inaxis)) # insanity checks if indiff < 0 or outdiff < 0: raise ValueError("impossible.") assert np.all(np.diff(spectral_grid) > 0) assert np.all(np.diff(inaxis) > 0) np.testing.assert_allclose(np.diff(spectral_grid), outdiff, err_msg="Output grid must be linear") if outdiff > 2 * indiff and not suppress_smooth_warning: warnings.warn("Input grid has too small a spacing. The data should " "be smoothed prior to resampling.", SmoothingWarning ) newspec = np.empty([spectral_grid.size], dtype=self.dtype) newmask = np.empty([spectral_grid.size], dtype='bool') newspec[outslice] = np.interp(spectral_grid.value, inaxis.value, self.filled_data[specslice].value, left=fill_value, right=fill_value) mask = self.mask.include() if all(mask): newmask = np.ones([spectral_grid.size], dtype='bool') else: interped = np.interp(spectral_grid.value, inaxis.value, mask[specslice]) > 0 newmask[outslice] = interped newwcs = self.wcs.deepcopy() newwcs.wcs.crpix[0] = 1 newwcs.wcs.crval[0] = spectral_grid[0].value if outslice.step > 0 \ else spectral_grid[-1].value newwcs.wcs.cunit[0] = spectral_grid.unit.to_string(format='FITS') newwcs.wcs.cdelt[0] = outdiff.value if outslice.step > 0 \ else -outdiff.value newwcs.wcs.set() newheader = self._nowcs_header.copy() newheader.update(newwcs.to_header()) wcs_cunit = u.Unit(newheader['CUNIT1']) newheader['CUNIT1'] = spectral_grid.unit.to_string(format='FITS') newheader['CDELT1'] *= wcs_cunit.to(spectral_grid.unit) newbmask = BooleanArrayMask(newmask, wcs=newwcs) return self._new_spectrum_with(data=newspec, wcs=newwcs, mask=newbmask, header=newheader, spectral_unit=spectral_grid.unit) def spectral_smooth(self, kernel, convolve=convolution.convolve, **kwargs): """ Smooth the spectrum Parameters ---------- kernel : `~astropy.convolution.Kernel1D` A 1D kernel from astropy convolve : function The astropy convolution function to use, either `astropy.convolution.convolve` or `astropy.convolution.convolve_fft` kwargs : dict Passed to the convolve function """ newspec = convolve(self.value, kernel, normalize_kernel=True, **kwargs) return self._new_spectrum_with(data=newspec) def to(self, unit, equivalencies=[]): """ Return a new `~spectral_cube.lower_dimensional_structures.OneDSpectrum` of the same class with the specified unit. See `astropy.units.Quantity.to` for further details. """ return super(BaseOneDSpectrum, self).to(unit, equivalencies, freq=None) def with_fill_value(self, fill_value): """ Create a new :class:`OneDSpectrum` with a different ``fill_value``. """ return self._new_spectrum_with(fill_value=fill_value) @property def _new_thing_with(self): return self._new_spectrum_with def _new_spectrum_with(self, data=None, wcs=None, mask=None, meta=None, fill_value=None, spectral_unit=None, unit=None, header=None, wcs_tolerance=None, **kwargs): data = self._data if data is None else data if unit is None and hasattr(data, 'unit'): if data.unit != self.unit: raise u.UnitsError("New data unit '{0}' does not" " match unit '{1}'. You can" " override this by specifying the" " `unit` keyword." .format(data.unit, self.unit)) unit = data.unit elif unit is None: unit = self.unit elif unit is not None: # convert string units to Units if not isinstance(unit, u.Unit): unit = u.Unit(unit) if hasattr(data, 'unit'): if u.Unit(unit) != data.unit: raise u.UnitsError("The specified new cube unit '{0}' " "does not match the input unit '{1}'." .format(unit, data.unit)) else: data = u.Quantity(data, unit=unit, copy=False) wcs = self._wcs if wcs is None else wcs mask = self._mask if mask is None else mask if meta is None: meta = {} meta.update(self._meta) if unit is not None: meta['BUNIT'] = unit.to_string(format='FITS') fill_value = self._fill_value if fill_value is None else fill_value spectral_unit = self._spectral_unit if spectral_unit is None else u.Unit(spectral_unit) spectrum = self.__class__(value=data, wcs=wcs, mask=mask, meta=meta, unit=unit, fill_value=fill_value, header=header or self._header, wcs_tolerance=wcs_tolerance or self._wcs_tolerance, **kwargs) spectrum._spectral_unit = spectral_unit return spectrum class OneDSpectrum(BaseOneDSpectrum, BeamMixinClass): def __new__(cls, value, beam=None, read_beam=False, **kwargs): self = super(OneDSpectrum, cls).__new__(cls, value, **kwargs) if beam is None: if "beam" in self.meta: beam = self.meta['beam'] elif read_beam: beam = cube_utils.try_load_beam(self.header) if beam is None: warnings.warn("Cannot load beam from header.", BeamWarning ) if beam is not None: self.beam = beam self.meta['beam'] = beam self._cache = {} return self def _new_spectrum_with(self, **kwargs): beam = kwargs.pop('beam', None) if 'beam' in self._meta and beam is None: beam = self.beam out = super(OneDSpectrum, self)._new_spectrum_with(beam=beam, **kwargs) return out def with_beam(self, beam, raise_error_jybm=True): ''' Attach a new beam object to the OneDSpectrum. Parameters ---------- beam : `~radio_beam.Beam` A new beam object. ''' if not isinstance(beam, Beam): raise TypeError("beam must be a radio_beam.Beam object.") self.check_jybeam_smoothing(raise_error_jybm=raise_error_jybm) meta = self.meta.copy() meta['beam'] = beam return self._new_spectrum_with(beam=beam, meta=meta) class VaryingResolutionOneDSpectrum(BaseOneDSpectrum, MultiBeamMixinClass): def __new__(cls, value, beams=None, read_beam=False, goodbeams_mask=None, **kwargs): self = super(VaryingResolutionOneDSpectrum, cls).__new__(cls, value, **kwargs) assert hasattr(self, '_fill_value') if beams is None: if "beams" in self.meta: beams = self.meta['beams'] elif read_beam: beams = cube_utils.try_load_beams(self.header) if beams is None: warnings.warn("Cannot load beams table from header.", BeamWarning ) if beams is not None: if isinstance(beams, BinTableHDU): beam_data_table = beams.data elif isinstance(beams, FITS_rec): beam_data_table = beams else: beam_data_table = None if beam_data_table is not None: beams = Beams(major=u.Quantity(beam_data_table['BMAJ'], u.arcsec), minor=u.Quantity(beam_data_table['BMIN'], u.arcsec), pa=u.Quantity(beam_data_table['BPA'], u.deg), meta=[{key: row[key] for key in beam_data_table.names if key not in ('BMAJ','BPA', 'BMIN')} for row in beam_data_table],) self.beams = beams self.meta['beams'] = beams if goodbeams_mask is not None: self.goodbeams_mask = goodbeams_mask self._cache = {} return self @property def hdu(self): warnings.warn("There are multiple beams for this spectrum that " "are being ignored when creating the HDU.", BeamWarning ) return super(VaryingResolutionOneDSpectrum, self).hdu @property def hdulist(self): with warnings.catch_warnings(): warnings.simplefilter("ignore") hdu = self.hdu beamhdu = cube_utils.beams_to_bintable(self.beams) return HDUList([hdu, beamhdu]) def _new_spectrum_with(self, **kwargs): beams = kwargs.pop('beams', self.beams) if beams is None: beams = self.beams VRODS = VaryingResolutionOneDSpectrum out = super(VRODS, self)._new_spectrum_with(beams=beams, **kwargs) return out def __array_finalize__(self, obj): super(VaryingResolutionOneDSpectrum, self).__array_finalize__(obj) self._beams = getattr(obj, '_beams', None) if getattr(obj, 'goodbeams_mask', None) is not None: # do NOT use the setter here, because we sometimes need to write # intermediate size-mismatch things that later get fixed, e.g., in # __getitem__ below self._goodbeams_mask = getattr(obj, 'goodbeams_mask', None) def __getitem__(self, key): new_qty = super(VaryingResolutionOneDSpectrum, self).__getitem__(key) # use the goodbeams_mask setter here because it checks size new_qty.goodbeams_mask = self.goodbeams_mask[key] new_qty.beams = self.unmasked_beams[key] return new_qty