diff --git a/.travis.yml b/.travis.yml index 0052f762..fa4e96ff 100644 --- a/.travis.yml +++ b/.travis.yml @@ -59,7 +59,7 @@ install: # Use pip to install development versions - pip install git+https://github.com/dendrograms/astrodendro.git - pip install git+https://github.com/radio-astro-tools/spectral-cube.git - - pip install git+https://github.com/radio-astro-tools/signal-id.git + # - pip install git+https://github.com/radio-astro-tools/signal-id.git - pip install git+https://github.com/radio-astro-tools/radio_beam.git - pip install emcee diff --git a/turbustat/cube_tools/obs_cubes.py b/turbustat/cube_tools/obs_cubes.py index dbb3b741..7b4cac58 100644 --- a/turbustat/cube_tools/obs_cubes.py +++ b/turbustat/cube_tools/obs_cubes.py @@ -9,146 +9,147 @@ from spectral_cube import SpectralCube, CompositeMask +from cube_utils import _check_mask, _check_beam, _get_int_intensity +# from cleaning_algs import * + try: from signal_id import Noise -except: - print("No signal_id package!") - pass - prefix = "/srv/astro/erickoch/" # Adjust if you're not me! - execfile(prefix + "Dropbox/code_development/signal-id/signal_id/noise.py") + class ObsCube(object): + """ + A wrapping class of SpectralCube which prepares observational data + cubes to be compared to any other data cube. -from cube_utils import _check_mask, _check_beam, _get_int_intensity -# from cleaning_algs import * + Parameters + ---------- + cube : str + Path to file. + mask : numpy.ndarray, any mask class from spectral_cube, optional + Mask to be applied to the cube. + algorithm : {NAME HERE}, optional + Name of the cleaning algorithm to use. -class ObsCube(object): - """ - A wrapping class of SpectralCube which prepares observational data cubes - to be compared to any other data cube. + Example + ------- + ``` + from turbustat.cube_tools import ObsCube - Parameters - ---------- + cube = ObsCube("data.fits") + cube.apply_cleaning(algorithm="SUPERAWESOMECLEANING") - cube : str - Path to file. - mask : numpy.ndarray, any mask class from spectral_cube, optional - Mask to be applied to the cube. - algorithm : {NAME HERE}, optional - Name of the cleaning algorithm to use. + ``` + """ + def __init__(self, cube, mask=None, algorithm=None, beam=None): - Example - ------- - ``` - from turbustat.cube_tools import ObsCube + raise NotImplementedError("ObsCube is not yet implemented for " + "general use.") - cube = ObsCube("data.fits") - cube.apply_cleaning(algorithm="SUPERAWESOMECLEANING") + super(ObsCube, self).__init__() + self.cube = SpectralCube.read(cube) - ``` - """ - def __init__(self, cube, mask=None, algorithm=None, beam=None): - super(ObsCube, self).__init__() - self.cube = sc.SpectralCube.read(cube) + self.algorithm = algorithm - self.algorithm = algorithm + # Make sure mask is an accepted type + if mask is not None: + _check_mask(mask) + self.mask = mask - # Make sure mask is an accepted type - if mask is not None: - _check_mask(mask) - self.mask = mask + if beam is not None: + _check_beam(beam) + self.noise = Noise(self.cube, beam=beam) - if beam is not None: - _check_beam(beam) - self.noise = Noise(self.cube, beam=beam) + def clean_cube(self, algorithm=None): + raise NotImplementedError("") - def clean_cube(self, algorithm=None): - raise NotImplementedError("") + def apply_mask(self, mask=None): + ''' + Check if the given mask is acceptable abd apply to + SpectralCube. + ''' - def apply_mask(self, mask=None): - ''' - Check if the given mask is acceptable abd apply to - SpectralCube. - ''' + # Update mask + if mask is not None: + _check_mask(mask) + self.mask = mask - # Update mask - if mask is not None: - _check_mask(mask) - self.mask = mask + # Create the mask, auto masking nan values + default_mask = np.isfinite(self.cube.filled_data[:]) + if self.mask is not None: + self.mask = CompositeMask(default_mask, self.mask) + else: + self.mask = default_mask - # Create the mask, auto masking nan values - default_mask = np.isfinite(self.cube.filled_data[:]) - if self.mask is not None: - self.mask = CompositeMask(default_mask, self.mask) - else: - self.mask = default_mask + # Apply mask to spectral cube object + self.cube = self.cube.with_mask(mask) - # Apply mask to spectral cube object - self.cube = self.cube.with_mask(mask) + return self - return self + def _update(self, data=None, wcs=None, beam=None, method="MAD"): + ''' + Helper function to update classes. + ''' - def _update(self, data=None, wcs=None, beam=None, method="MAD"): - ''' - Helper function to update classes. - ''' + # Check if we need a new SpectralCube + if data is None and wcs is None: + pass + else: + if data is None: + data = self.cube.unmasked_data[:] + if wcs is None: + wcs = self.cube.wcs + # Make new SpectralCube object + self.cube = SpectralCube(data=data, wcs=wcs) - # Check if we need a new SpectralCube - if data is None and wcs is None: - pass - else: - if data is None: - data = self.cube.unmasked_data[:] - if wcs is None: - wcs = self.cube.wcs - # Make new SpectralCube object - self.cube = SpectralCube(data=data, wcs=wcs) + if beam is not None: + _check_beam(beam) + self.noise = Noise(self.cube, beam=beam, method=method) - if beam is not None: - _check_beam(beam) - self.noise = Noise(self.cube, beam=beam, method=method) + def compute_properties(self): + ''' + Use SpectralCube to compute the moments. Also compute the integrated + intensity based on the noise properties from Noise. + ''' - def compute_properties(self): - ''' - Use SpectralCube to compute the moments. Also compute the integrated - intensity based on the noise properties from Noise. - ''' + self._moment0 = self.cube.moment0().value - self._moment0 = self.cube.moment0().value + self._moment1 = self.cube.moment1().value - self._moment1 = self.cube.moment1().value + self._moment2 = self.cube.moment2().value - self._moment2 = self.cube.moment2().value + _get_int_intensity(self) - _get_int_intensity(self) + return self - return self + @property + def moment0(self): + return self._moment0 - @property - def moment0(self): - return self._moment0 + @property + def moment1(self): + return self._moment1 - @property - def moment1(self): - return self._moment1 + @property + def moment2(self): + return self._moment2 - @property - def moment2(self): - return self._moment2 + @property + def intint(self): + return self._intint - @property - def intint(self): - return self._intint + def prep(self, mask=None, algorithm=None): + ''' + Prepares the cube to be compared to another cube. + ''' - def prep(self, mask=None, algorithm=None): - ''' - Prepares the cube to be compared to another cube. - ''' + if not mask is None: + self.apply_mask() - if not mask is None: - self.apply_mask() + self.clean_cube() + self.compute_properties() - self.clean_cube() - self.compute_properties() + return self - return self +except ImportError: + # raise ImportError("The signal_id package must be installed") + pass diff --git a/turbustat/cube_tools/sim_cubes.py b/turbustat/cube_tools/sim_cubes.py index 54eb6349..bf1deb15 100644 --- a/turbustat/cube_tools/sim_cubes.py +++ b/turbustat/cube_tools/sim_cubes.py @@ -9,159 +9,162 @@ from spectral_cube import SpectralCube, CompositeMask +from cube_utils import _check_mask, _check_beam, _get_int_intensity + try: from signal_id import Noise -except ImportError: - pass - prefix = "/srv/astro/erickoch/" # Adjust if you're not me! - execfile(prefix + "Dropbox/code_development/signal-id/signal_id/noise.py") - -from cube_utils import _check_mask, _check_beam, _get_int_intensity + class SimCube(object): + ''' + A wrapping class to prepare a simulated spectral data cube for + comparison with another cube. + ''' -class SimCube(object): - ''' - A wrapping class to prepare a simulated spectral data cube for - comparison with another cube. - ''' + def __init__(self, cube, beam=None, mask=None, method="MAD", + compute=True): - def __init__(self, cube, beam=None, mask=None, method="MAD", compute=True): + raise NotImplementedError("SimCube is not yet implemented for " + "general use.") - # Initialize cube object - self.cube = SpectralCube.read(cube) + # Initialize cube object + self.cube = SpectralCube.read(cube) - if mask is not None: - _check_mask(mask) - self.mask = mask + if mask is not None: + _check_mask(mask) + self.mask = mask - if beam is not None: - _check_beam(mask) + if beam is not None: + _check_beam(mask) - # Initialize noise object - self.noise = Noise(self.cube, beam=beam, method=method) + # Initialize noise object + self.noise = Noise(self.cube, beam=beam, method=method) - def add_noise(self): - ''' - Use Noise to add synthetic noise to the data. Then update - SpectralCube. - ''' + def add_noise(self): + ''' + Use Noise to add synthetic noise to the data. Then update + SpectralCube. + ''' - # Create the noisy cube - self.noise.get_noise_cube() - noise_data = self.noise.noise_cube +\ - self.cube.filled_data[:] + # Create the noisy cube + self.noise.get_noise_cube() + noise_data = self.noise.noise_cube +\ + self.cube.filled_data[:] - # Update SpectralCube object - self._update(data=noise_data) + # Update SpectralCube object + self._update(data=noise_data) - return self + return self + + def clean_cube(self, algorithm=None): + raise NotImplementedError("") + + def apply_mask(self, mask=None): + ''' + Check if the given mask is acceptable abd apply to + SpectralCube. + ''' - def clean_cube(self, algorithm=None): - raise NotImplementedError("") + # Update mask + if mask is not None: + _check_mask(mask) + self.mask = mask - def apply_mask(self, mask=None): - ''' - Check if the given mask is acceptable abd apply to - SpectralCube. - ''' + # Create the mask, auto masking nan values + default_mask = np.isfinite(self.cube.filled_data[:]) + if self.mask is not None: + self.mask = CompositeMask(default_mask, self.mask) + else: + self.mask = default_mask - # Update mask - if mask is not None: - _check_mask(mask) - self.mask = mask + # Apply mask to spectral cube object + self.cube = self.cube.with_mask(mask) - # Create the mask, auto masking nan values - default_mask = np.isfinite(self.cube.filled_data[:]) - if self.mask is not None: - self.mask = CompositeMask(default_mask, self.mask) - else: - self.mask = default_mask + return self - # Apply mask to spectral cube object - self.cube = self.cube.with_mask(mask) + def _update(self, data=None, wcs=None, beam=None, method="MAD"): + ''' + Helper function to update classes. + ''' - return self + # Check if we need a new SpectralCube + if data is None and wcs is None: + pass + else: + if data is None: + data = self.cube.unmasked_data[:] + if wcs is None: + wcs = self.cube.wcs + # Make new SpectralCube object + self.cube = SpectralCube(data=data, wcs=wcs) - def _update(self, data=None, wcs=None, beam=None, method="MAD"): - ''' - Helper function to update classes. - ''' + if beam is not None: + _check_beam(beam) + self.noise = Noise(self.cube, beam=beam, method=method) - # Check if we need a new SpectralCube - if data is None and wcs is None: - pass - else: - if data is None: - data = self.cube.unmasked_data[:] - if wcs is None: - wcs = self.cube.wcs - # Make new SpectralCube object - self.cube = SpectralCube(data=data, wcs=wcs) - - if beam is not None: - _check_beam(beam) - self.noise = Noise(self.cube, beam=beam, method=method) + def compute_properties(self): + ''' + Use SpectralCube to compute the moments. Also compute the + integrated intensity based on the noise properties from Noise. + ''' - def compute_properties(self): - ''' - Use SpectralCube to compute the moments. Also compute the integrated - intensity based on the noise properties from Noise. - ''' + self._moment0 = self.cube.moment0().value - self._moment0 = self.cube.moment0().value + self._moment1 = self.cube.moment1().value - self._moment1 = self.cube.moment1().value + self._moment2 = self.cube.moment2().value - self._moment2 = self.cube.moment2().value + _get_int_intensity(self) - _get_int_intensity(self) + return self - return self + @property + def moment0(self): + return self._moment0 - @property - def moment0(self): - return self._moment0 + @property + def moment1(self): + return self._moment1 - @property - def moment1(self): - return self._moment1 + @property + def moment2(self): + return self._moment2 - @property - def moment2(self): - return self._moment2 + @property + def intint(self): + return self._intint - @property - def intint(self): - return self._intint + def sim_prep(self, mask=None): + ''' + Prepares the cube when being compared to another simulation. + This entails: + * Optionally applying a mask to the data. + * Computing the cube's property arrays + ''' - def sim_prep(self, mask=None): - ''' - Prepares the cube when being compared to another simulation. - This entails: - * Optionally applying a mask to the data. - * Computing the cube's property arrays - ''' + if not mask is None: + self.apply_mask() - if not mask is None: - self.apply_mask() + self.compute_properties() - self.compute_properties() + return self - return self + def obs_prep(self, mask=None): + ''' + Prepares the cube when being compared to observational data. + This entails: + * Optionally applying a mask to the data. + * Add synthetic noise based on the cube's properties. + * Computing the cube's property arrays + ''' - def obs_prep(self, mask=None): - ''' - Prepares the cube when being compared to observational data. - This entails: - * Optionally applying a mask to the data. - * Add synthetic noise based on the cube's properties. - * Computing the cube's property arrays - ''' + if not mask is None: + self.apply_mask() - if not mask is None: - self.apply_mask() + self.add_noise() + self.compute_properties() - self.add_noise() - self.compute_properties() + return self - return self +except ImportError: + # raise ImportError("The signal_id package must be installed") + pass diff --git a/turbustat/data_reduction/make_moments.py b/turbustat/data_reduction/make_moments.py index 62da3ba0..48e14608 100644 --- a/turbustat/data_reduction/make_moments.py +++ b/turbustat/data_reduction/make_moments.py @@ -44,7 +44,7 @@ class Mask_and_Moments(object): for deriving the noise level. clip : float, optional Sigma level to clip data at. - scale : float, optional + scale : `~astropy.units.Quantity`, optional The noise level in the cube. Overrides estimation using `signal_id `_ moment_method : {'slice', 'cube', 'ray'}, optional @@ -369,7 +369,7 @@ def to_fits(self, save_name=None): @staticmethod def from_fits(fits_name, moments_prefix=None, moments_path=None, mask_name=None, moment0=None, centroid=None, linewidth=None, - intint=None): + intint=None, scale=None): ''' Load pre-made moment arrays given a cube name. Saved moments must match the naming of the cube for the automatic loading to work @@ -401,6 +401,9 @@ def from_fits(fits_name, moments_prefix=None, moments_path=None, intint : str, optional Filename of the integrated intensity array. Use if naming scheme is not valid for automatic loading. + scale : `~astropy.units.Quantity`, optional + The noise level in the cube. Overrides estimation using + `signal_id `_ ''' if not spectral_cube_flag: @@ -416,7 +419,7 @@ def from_fits(fits_name, moments_prefix=None, moments_path=None, else: root_name = moments_prefix - self = Mask_and_Moments(fits_name) + self = Mask_and_Moments(fits_name, scale=scale) if mask_name is not None: mask = fits.getdata(mask_name) diff --git a/turbustat/tests/_testing_data.py b/turbustat/tests/_testing_data.py index da56e654..2d089468 100644 --- a/turbustat/tests/_testing_data.py +++ b/turbustat/tests/_testing_data.py @@ -21,6 +21,7 @@ from astropy.io.fits import getheader from astropy.wcs import WCS from astropy.io import fits +import astropy.units as u # Set seed for adding noise. ra.seed(121212) @@ -54,7 +55,8 @@ sc1 = SpectralCube(data=cube1, wcs=WCS(header)) mask = LazyMask(np.isfinite, sc1) sc1 = sc1.with_mask(mask) -props1 = Mask_and_Moments(sc1) +# Set the scale for the purposes of the tests +props1 = Mask_and_Moments(sc1, scale=0.003031065017916262 * u.Unit("")) props1.make_mask(mask=mask) props1.make_moments() props1.make_moment_errors() @@ -90,7 +92,8 @@ sc2 = SpectralCube(data=cube2, wcs=WCS(header)) mask = LazyMask(np.isfinite, sc2) sc2 = sc2.with_mask(mask) -props2 = Mask_and_Moments(sc2) +# Set the scale for the purposes of the tests +props2 = Mask_and_Moments(sc2, scale=0.003029658694658428 * u.Unit("")) props2.make_moments() props2.make_moment_errors() diff --git a/turbustat/tests/test_make_moments.py b/turbustat/tests/test_make_moments.py index 94c4ced7..4c5bada5 100644 --- a/turbustat/tests/test_make_moments.py +++ b/turbustat/tests/test_make_moments.py @@ -1,32 +1,28 @@ # Licensed under an MIT open source license - see LICENSE -from unittest import TestCase - import numpy.testing as npt - +import astropy.units as u from ..data_reduction import Mask_and_Moments from ._testing_data import dataset1, sc1 -class test_Mask_and_Moments(TestCase): - """docstring for test_Mask_and_Moments""" - - def test_loading(self): - - # Try loading the files. +def test_loading(): - test = Mask_and_Moments.from_fits(sc1, moments_prefix="dataset1", - moments_path=".") + # Try loading the files. + # Set the scale to the assumed value. + test = Mask_and_Moments.from_fits(sc1, moments_prefix="dataset1", + moments_path=".", + scale=0.003031065017916262 * u.Unit("")) - npt.assert_allclose(test.moment0, dataset1["moment0"][0]) - npt.assert_allclose(test.moment1, dataset1["centroid"][0]) - npt.assert_allclose(test.linewidth, dataset1["linewidth"][0]) - npt.assert_allclose(test.intint, - dataset1["integrated_intensity"][0]) + npt.assert_allclose(test.moment0, dataset1["moment0"][0]) + npt.assert_allclose(test.moment1, dataset1["centroid"][0]) + npt.assert_allclose(test.linewidth, dataset1["linewidth"][0]) + npt.assert_allclose(test.intint, + dataset1["integrated_intensity"][0]) - npt.assert_allclose(test.moment0_err, dataset1["moment0_error"][0]) - npt.assert_allclose(test.moment1_err, dataset1["centroid_error"][0]) - npt.assert_allclose(test.linewidth_err, dataset1["linewidth_error"][0]) - npt.assert_allclose(test.intint_err, - dataset1["integrated_intensity_error"][0]) + npt.assert_allclose(test.moment0_err, dataset1["moment0_error"][0]) + npt.assert_allclose(test.moment1_err, dataset1["centroid_error"][0]) + npt.assert_allclose(test.linewidth_err, dataset1["linewidth_error"][0]) + npt.assert_allclose(test.intint_err, + dataset1["integrated_intensity_error"][0])