From 41f2b38e829c40f7c0f72bcb54f01decb6f7d74c Mon Sep 17 00:00:00 2001 From: e-koch Date: Wed, 17 Jan 2018 15:56:55 -0700 Subject: [PATCH 1/2] Test constant velocity offsets in stacking --- .../tests/test_analysis_functions.py | 51 ++++++++++++++++++- spectral_cube/tests/utilities.py | 8 ++- 2 files changed, 55 insertions(+), 4 deletions(-) diff --git a/spectral_cube/tests/test_analysis_functions.py b/spectral_cube/tests/test_analysis_functions.py index 5bfca9a7a..8b65e2826 100644 --- a/spectral_cube/tests/test_analysis_functions.py +++ b/spectral_cube/tests/test_analysis_functions.py @@ -4,8 +4,8 @@ import astropy.units as u # from astropy.modeling import models, fitting -from ..analysis_utilities import stack_spectra, fourier_shift -from .utilities import generate_gaussian_cube, gaussian +from spectral_cube.analysis_utilities import stack_spectra, fourier_shift +from spectral_cube.tests.utilities import generate_gaussian_cube, gaussian def test_shift(): @@ -128,6 +128,53 @@ def test_stacking_wpadding(): or (stacked.size == stack_shape + 1) +def test_stacking_woffset(): + ''' + Use a set of identical Gaussian profiles randomly offset to ensure the + shifted spectrum has the correct properties. + + Make sure the operations aren't affected by absolute velocity offsets + + ''' + + amp = 1. + sigma = 8. + v0 = 100. * u.km / u.s + noise = None + shape = (100, 25, 25) + + test_cube, test_vels = \ + generate_gaussian_cube(shape=shape, amp=amp, sigma=sigma, noise=noise, + v0=v0.value) + + # Stack the spectra in the cube + stacked = \ + stack_spectra(test_cube, test_vels, v0=v0, + stack_function=np.nanmean, + xy_posns=None, num_cores=1, + chunk_size=-1, + progressbar=False, pad_edges=True) + + true_spectrum = gaussian(stacked.spectral_axis.value, + amp, v0.value, sigma) + + # Calculate residuals + resid = np.abs(stacked.value - true_spectrum) + assert np.std(resid) <= 1e-3 + + # Now fit a Gaussian to the mean stacked profile. + # fit_vals = fit_gaussian(stacked.spectral_axis.value, stacked.value)[0] + + # np.testing.assert_allclose(fit_vals, np.array([amp, 0.0, sigma]), + # atol=1e-3) + + # The spectral axis should be padded by ~25% on each side + stack_shape = int(test_cube.shape[0] * 1.5) + # This is rounded, so the shape could be +/- 1 + assert (stacked.size == stack_shape) or (stacked.size == stack_shape - 1) \ + or (stacked.size == stack_shape + 1) + + def test_stacking_noisy(): # Test stack w/ S/N of 0.2 diff --git a/spectral_cube/tests/utilities.py b/spectral_cube/tests/utilities.py index 4f332b180..370bebbf8 100644 --- a/spectral_cube/tests/utilities.py +++ b/spectral_cube/tests/utilities.py @@ -55,6 +55,7 @@ def generate_gaussian_cube(shape=(100, 25, 25), sigma=8., amp=1., noise=None, spec_scale=1 * u.km / u.s, pixel_scale=1 * u.arcsec, beamfwhm=3 * u.arcsec, + v0=None, vel_surface=None, seed=247825498): ''' @@ -79,6 +80,9 @@ def generate_gaussian_cube(shape=(100, 25, 25), sigma=8., amp=1., spec_middle = int(shape[0] / 2) spec_quarter = int(shape[0] / 4) + if v0 is None: + v0 = 0 + with NumpyRNGContext(seed): spec_inds = np.mgrid[-spec_middle:spec_middle] * spec_scale.value @@ -88,13 +92,13 @@ def generate_gaussian_cube(shape=(100, 25, 25), sigma=8., amp=1., mean_pos = \ np.random.uniform(low=spec_inds[spec_quarter], high=spec_inds[spec_quarter + spec_middle]) - mean_positions[y, x] = mean_pos test_cube[:, y, x] = gaussian(spec_inds, amp, mean_pos, sigma) + mean_positions[y, x] = mean_pos + v0 if noise is not None: test_cube[:, y, x] += np.random.normal(0, noise, shape[0]) test_hdu = generate_hdu(test_cube, pixel_scale, spec_scale, beamfwhm, - spec_inds[0]) + spec_inds[0] + v0) spec_cube = SpectralCube.read(test_hdu) From 4ef4b5dffbee080b493dd8afdce17c2286dee0bd Mon Sep 17 00:00:00 2001 From: e-koch Date: Wed, 17 Jan 2018 16:05:03 -0700 Subject: [PATCH 2/2] Fix stacking with a decreasing spectral axis; add test for this case --- spectral_cube/analysis_utilities.py | 10 ++-- .../tests/test_analysis_functions.py | 46 +++++++++++++++---- 2 files changed, 45 insertions(+), 11 deletions(-) diff --git a/spectral_cube/analysis_utilities.py b/spectral_cube/analysis_utilities.py index 154719026..44a4b51fd 100644 --- a/spectral_cube/analysis_utilities.py +++ b/spectral_cube/analysis_utilities.py @@ -205,7 +205,11 @@ def stack_spectra(cube, velocity_surface, v0=None, "axis of the cube.") # Calculate the pixel shifts that will be applied. - vdiff = np.abs(np.diff(cube.spectral_axis[:2])[0]) + spec_size = np.diff(cube.spectral_axis[:2])[0] + # Assign the correct +/- for pixel shifts based on whether the spectral + # axis is increasing (-1) or decreasing (+1) + vdiff_sign = -1. if spec_size.value > 0. else 1. + vdiff = np.abs(spec_size) vel_unit = vdiff.unit # Check to make sure vdiff doesn't change more than the allowed tolerance @@ -214,8 +218,8 @@ def stack_spectra(cube, velocity_surface, v0=None, if not np.isclose(vdiff2.value, vdiff.value, rtol=vdiff_tol): raise ValueError("Cannot shift spectra on a non-linear axes") - pix_shifts = -((velocity_surface.to(vel_unit) - - v0.to(vel_unit)) / vdiff).value[xy_posns] + pix_shifts = vdiff_sign * ((velocity_surface.to(vel_unit) - + v0.to(vel_unit)) / vdiff).value[xy_posns] # May a header copy so we can start altering new_header = cube[:, 0, 0].header.copy() diff --git a/spectral_cube/tests/test_analysis_functions.py b/spectral_cube/tests/test_analysis_functions.py index 8b65e2826..4d46b2e77 100644 --- a/spectral_cube/tests/test_analysis_functions.py +++ b/spectral_cube/tests/test_analysis_functions.py @@ -4,8 +4,8 @@ import astropy.units as u # from astropy.modeling import models, fitting -from spectral_cube.analysis_utilities import stack_spectra, fourier_shift -from spectral_cube.tests.utilities import generate_gaussian_cube, gaussian +from ..analysis_utilities import stack_spectra, fourier_shift +from .utilities import generate_gaussian_cube, gaussian def test_shift(): @@ -85,6 +85,42 @@ def test_stacking(): test_cube.spectral_axis.value) +def test_stacking_reversed_specaxis(): + ''' + Use a set of identical Gaussian profiles randomly offset to ensure the + shifted spectrum has the correct properties. + ''' + + amp = 1. + v0 = 0. * u.km / u.s + sigma = 8. + noise = None + shape = (100, 25, 25) + + test_cube, test_vels = \ + generate_gaussian_cube(amp=amp, sigma=sigma, noise=noise, + shape=shape, spec_scale=-1. * u.km / u.s) + + true_spectrum = gaussian(test_cube.spectral_axis.value, + amp, v0.value, sigma) + + # Stack the spectra in the cube + stacked = \ + stack_spectra(test_cube, test_vels, v0=v0, + stack_function=np.nanmean, + xy_posns=None, num_cores=1, + chunk_size=-1, + progressbar=False, pad_edges=False) + + # Calculate residuals + resid = np.abs(stacked.value - true_spectrum) + assert np.std(resid) <= 1e-3 + + # The stacked spectrum should have the same spectral axis + np.testing.assert_allclose(stacked.spectral_axis.value, + test_cube.spectral_axis.value) + + def test_stacking_wpadding(): ''' Use a set of identical Gaussian profiles randomly offset to ensure the @@ -162,12 +198,6 @@ def test_stacking_woffset(): resid = np.abs(stacked.value - true_spectrum) assert np.std(resid) <= 1e-3 - # Now fit a Gaussian to the mean stacked profile. - # fit_vals = fit_gaussian(stacked.spectral_axis.value, stacked.value)[0] - - # np.testing.assert_allclose(fit_vals, np.array([amp, 0.0, sigma]), - # atol=1e-3) - # The spectral axis should be padded by ~25% on each side stack_shape = int(test_cube.shape[0] * 1.5) # This is rounded, so the shape could be +/- 1