From 494f6188bc745367474b8d529e22973d69382ca0 Mon Sep 17 00:00:00 2001 From: Clare Shanahan Date: Mon, 22 Apr 2024 09:16:51 -0400 Subject: [PATCH] . --- specreduce/background.py | 25 ++- specreduce/core.py | 1 - specreduce/tests/test_background.py | 83 +++---- specreduce/tests/test_extract.py | 60 +++++ specreduce/tests/test_tracing.py | 327 ++++++++++++++++++---------- specreduce/tracing.py | 76 +++++-- 6 files changed, 384 insertions(+), 188 deletions(-) diff --git a/specreduce/background.py b/specreduce/background.py index db8ced1..db272a9 100644 --- a/specreduce/background.py +++ b/specreduce/background.py @@ -98,16 +98,27 @@ def __post_init__(self): dispersion axis crossdisp_axis : int cross-dispersion axis + mask_treatment : string + The method for handling masked or non-finite data. Choice of `filter`, + `omit`, or `zero-fill`. If `filter` is chosen, masked/non-finite data + will be filtered during the fit to each bin/column (along disp. axis) to + find the peak. If `omit` is chosen, columns along disp_axis with any + masked/non-finite data values will be fully masked (i.e, 2D mask is + collapsed to 1D and applied). If `zero-fill` is chosen, masked/non-finite + data will be replaced with 0.0 in the input image, and the mask will then + be dropped. For all three options, the input mask (optional on input + NDData object) will be combined with a mask generated from any non-finite + values in the image data. """ # Parse image, including masked/nonfinite data handling based on - # choice of `mask_treatment`. returns a Spectrum1D + # choice of `mask_treatment`. Any uncaught nonfinte data values will be + # masked as well. Returns a Spectrum1D. self.image = self._parse_image(self.image) - # _parse_image returns a Spectrum1D. convert this to a masked array - # for ease of calculations here (even if there is no masked data). - # Note: uncertainties are dropped, this should also be addressed at - # some point probably across the package. + # always work with masked array, even if there is no masked + # or nonfinite data, in case padding is needed. if not, mask will be + # dropped at the end and a regular array will be returned. img = np.ma.masked_array(self.image.data, self.image.mask) if self.width < 0: @@ -120,6 +131,9 @@ def __post_init__(self): bkg_wimage = np.zeros_like(self.image.data, dtype=np.float64) for trace in self.traces: + # note: ArrayTrace can have masked values, but if it does a MaskedArray + # will be returned so this should be reflected in the window size here + # (i.e, np.nanmax is not required.) windows_max = trace.trace.data.max() + self.width/2 windows_min = trace.trace.data.min() - self.width/2 if windows_max >= self.image.shape[self.crossdisp_axis]: @@ -136,6 +150,7 @@ def __post_init__(self): self.crossdisp_axis, self.image.shape) + if np.any(bkg_wimage > 1): raise ValueError("background regions overlapped") if np.any(np.sum(bkg_wimage, axis=self.crossdisp_axis) == 0): diff --git a/specreduce/core.py b/specreduce/core.py index 87b801d..9b6eae5 100644 --- a/specreduce/core.py +++ b/specreduce/core.py @@ -139,7 +139,6 @@ def _mask_and_nonfinite_data_handling(self, image, mask): # nothing needs to be done. input mask (combined with nonfinite data) # remains with data as-is. - if mask_treatment == 'zero-fill': # make a copy of the input image since we will be modifying it image = deepcopy(image) diff --git a/specreduce/tests/test_background.py b/specreduce/tests/test_background.py index 753f881..99a53d2 100644 --- a/specreduce/tests/test_background.py +++ b/specreduce/tests/test_background.py @@ -24,13 +24,13 @@ def test_background(mk_test_img_raw, mk_test_spec_no_spectral_axis, # all the following should be equivalent, whether image's spectral axis # is in pixels or physical units: - bg1 = Background(image, [trace-bkg_sep, trace+bkg_sep], width=bkg_width) + bg1 = Background(image, [trace - bkg_sep, trace + bkg_sep], width=bkg_width) bg2 = Background.two_sided(image, trace, bkg_sep, width=bkg_width) bg3 = Background.two_sided(image, trace_pos, bkg_sep, width=bkg_width) assert np.allclose(bg1.bkg_image().flux, bg2.bkg_image().flux) assert np.allclose(bg1.bkg_image().flux, bg3.bkg_image().flux) - bg4 = Background(image_um, [trace-bkg_sep, trace+bkg_sep], width=bkg_width) + bg4 = Background(image_um, [trace - bkg_sep, trace + bkg_sep], width=bkg_width) bg5 = Background.two_sided(image_um, trace, bkg_sep, width=bkg_width) bg6 = Background.two_sided(image_um, trace_pos, bkg_sep, width=bkg_width) assert np.allclose(bg1.bkg_image().flux, bg4.bkg_image().flux) @@ -74,7 +74,7 @@ def test_background(mk_test_img_raw, mk_test_spec_no_spectral_axis, stats = ['average', 'median'] for st in stats: - bg = Background(img, trace-bkg_sep, width=bkg_width, statistic=st) + bg = Background(img, trace - bkg_sep, width=bkg_width, statistic=st) assert np.isnan(bg.image.flux).sum() == 2 assert np.isnan(bg._bkg_array).sum() == 0 assert np.isnan(bg.bkg_spectrum().flux).sum() == 0 @@ -101,7 +101,7 @@ def test_warnings_errors(mk_test_spec_no_spectral_axis): with pytest.warns(match="background window extends beyond image boundaries"): Background.two_sided(image, 7, 5, width=6) - trace = ArrayTrace(image, trace=np.arange(10)+20) # from 20 to 29 + trace = ArrayTrace(image, trace=np.arange(10) + 20) # from 20 to 29 with pytest.warns(match="background window extends beyond image boundaries"): with pytest.raises(ValueError, match="background window does not remain in bounds across entire dispersion axis"): # noqa @@ -112,8 +112,8 @@ def test_warnings_errors(mk_test_spec_no_spectral_axis): with pytest.raises(ValueError, match="width must be positive"): Background.two_sided(image, 25, 2, width=-1) -def test_trace_inputs(mk_test_img_raw): +def test_trace_inputs(mk_test_img_raw): """ Tests for the input argument 'traces' to `Background`. This should accept a list of or a single Trace object, or a list of or a single (positive) @@ -143,20 +143,23 @@ def test_trace_inputs(mk_test_img_raw): with pytest.raises(ValueError, match=match_str): Background(image, 'non_valid_trace_pos') + class TestMasksBackground(): """ Various test functions to test how masked and non-finite data is handled - in `Background. + in `Background. There are three currently implemented options for masking + in Background: filter, omit, and zero-fill. """ def mk_img(self, nrows=4, ncols=5, nan_slices=None): """ - Make a simpleimage to test masking in Background. - Optionally add NaNs to data. Returned array is in u.DN. + Make a simple gradient image to test masking in Background. + Optionally add NaNs to data with `nan_slices`. Returned array is in + u.DN. """ - img = np.tile((np.arange(1., ncols+1)), (nrows, 1)) + img = np.tile((np.arange(1., ncols + 1)), (nrows, 1)) if nan_slices: # add nans in data for s in nan_slices: @@ -164,22 +167,23 @@ def mk_img(self, nrows=4, ncols=5, nan_slices=None): return img * u.DN - def test_fully_masked_column(self): + @pytest.mark.parametrize("mask", ["filter", "omit", "zero-fill"]) + def test_fully_masked_column(self, mask): """ - Test what happens when a full column is masked, not the entire - image. In this case, the background value for that fully-masked - column should be 0.0, with no error or warning raised. + Test background with some fully-masked columns (not fully masked image). + In this case, the background value for that fully-masked column should + be 0.0, with no error or warning raised. This is the case for + mask_treatment=filter, omit, or zero-fill. """ - img = np.ones((12, 12)) + img = self.mk_img(nrows=10, ncols=10) img[:, 0:1] = np.nan - bkg = Background(img, traces=FlatTrace(img, 6)) - + bkg = Background(img, traces=FlatTrace(img, 6), mask_treatment=mask) assert np.all(bkg.bkg_image().data[:, 0:1] == 0.0) - - def test_fully_masked(self): + @pytest.mark.parametrize("mask", ["filter", "omit", "zero-fill"]) + def test_fully_masked_image(self, mask): """ Test that the appropriate error is raised by `Background` when image is fully masked/NaN. @@ -187,16 +191,18 @@ def test_fully_masked(self): with pytest.raises(ValueError, match='Image is fully masked.'): # fully NaN image - img = np.zeros((4, 5)) * np.nan - Background(img, traces=FlatTrace(self.mk_img(), 2)) + img = self.mk_img() * np.nan + Background(img, traces=FlatTrace(self.mk_img(), 2), mask_treatment=mask) with pytest.raises(ValueError, match='Image is fully masked.'): - # fully masked image (should be equivilant) + # fully masked image (should be equivalent) img = NDData(np.ones((4, 5)), mask=np.ones((4, 5))) - Background(img, traces=FlatTrace(self.mk_img(), 2)) + Background(img, traces=FlatTrace(self.mk_img(), 2), mask_treatment=mask) # Now test that an image that isn't fully masked, but is fully masked - # within the window determined by `width`, produces the correct result + # within the window determined by `width`, produces the correct result. + # only applicable for mask_treatment=filter, because this is the only + # option that allows a slice of masked values that don't span all rows. msg = 'Image is fully masked within background window determined by `width`.' with pytest.raises(ValueError, match=msg): img = self.mk_img(nrows=12, ncols=12, nan_slices=[np.s_[3:10, :]]) @@ -204,17 +210,17 @@ def test_fully_masked(self): @pytest.mark.filterwarnings("ignore:background window extends beyond image boundaries") @pytest.mark.parametrize("method,expected", - [("filter", np.array([1., 2., 3., 4., 5., 6., 7., + [("filter", np.array([1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12.])), - ("omit", np.array([0., 2., 3., 0., 5., 6., - 7., 0., 9., 10., 11., 12.])), - ("zero-fill", np.array([ 0.58333333, 2., 3., - 2.33333333, 5., 6., 7., - 7.33333333, 9., 10., 11., - 12.]))]) + ("omit", np.array([0., 2., 3., 0., 5., 6., + 7., 0., 9., 10., 11., 12.])), + ("zero-fill", np.array([0.58333333, 2., 3., + 2.33333333, 5., 6., 7., + 7.33333333, 9., 10., 11., + 12.]))]) def test_mask_treatment_bkg_img_spectrum(self, method, expected): - """ - This test function tests `Backgroud.bkg_image` and + """ + This test function tests `Backgroud.bkg_image` and `Background.bkg_spectrum` when there is masked data. It also tests background subtracting the image, and returning the spectrum of the background subtracted image. This test is parameterized over all @@ -227,8 +233,8 @@ def test_mask_treatment_bkg_img_spectrum(self, method, expected): # make image, set some value to nan, which will be masked in the function image1 = self.mk_img(nrows=img_size, ncols=img_size, - nan_slices=[np.s_[5:10, 0], np.s_[7:12, 3], - np.s_[2, 7]]) + nan_slices=[np.s_[5:10, 0], np.s_[7:12, 3], + np.s_[2, 7]]) # also make an image that doesn't have nonf data values, but has # masked values at the same locations, to make sure they give the same @@ -240,13 +246,13 @@ def test_mask_treatment_bkg_img_spectrum(self, method, expected): for image in [image1, image2]: # construct a flat trace in center of image - trace = FlatTrace(image, img_size/2) + trace = FlatTrace(image, img_size / 2) # create 'Background' object with `mask_treatment` set # 'width' should be > size of image to use all pix (but warning will # be raised, which we ignore.) background = Background(image, mask_treatment=method, - traces=trace, width=img_size+1) + traces=trace, width=img_size + 1) # test background image matches 'expected' bk_img = background.bkg_image() @@ -265,7 +271,7 @@ def test_sub_bkg_image(self): """ Test that masked and nonfinite data is handled correctly when subtracting background from image, for all currently implemented masking - options ('filter', 'omit', and 'zero-fill'). + options ('filter', 'omit', and 'zero-fill'). """ # make image, set some value to nan, which will be masked in the function @@ -288,7 +294,7 @@ def test_sub_bkg_image(self): # 2d mask is reduced to a 1d mask to mask out full columns in the # presence of any nans - this means that (as tested above in # `test_mask_treatment_bkg_img_spectrum`) those columns will have 0.0 - # background. In this case, image.mask is expanded to mask full + # background. In this case, image.mask is expanded to mask full # columns - the image itself will not have full columns set to np.nan, # so there are still valid background subtracted data values in this # case, but the corresponding mask for that entire column will be masked. @@ -314,4 +320,3 @@ def test_sub_bkg_image(self): assert np.all(np.isfinite(subtracted_img_zero_fill.data)) assert np.all(subtracted_img_zero_fill.mask == 0) - diff --git a/specreduce/tests/test_extract.py b/specreduce/tests/test_extract.py index 47c5716..ece5a4d 100644 --- a/specreduce/tests/test_extract.py +++ b/specreduce/tests/test_extract.py @@ -81,6 +81,7 @@ def test_boxcar_extraction(mk_test_img): assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 67.15)) + def test_boxcar_outside_image_condition(mk_test_img): # # Trace is such that extraction aperture lays partially outside the image @@ -326,3 +327,62 @@ def test_horne_interpolated_nbins_fails(mk_test_img): spatial_profile={'name': 'interpolated_profile', 'n_bins_interpolated_profile': 100}) ex.spectrum + +class TestMasksExtract(): + + def mk_flat_gauss_img(nrows=200, ncols=160, nan_slices=None, add_noise=True): + + """ + Makes a flat gaussian image for testing, with optional added gaussian + nosie and optional data values set to NaN. Variance is included, which + is required by HorneExtract. Returns a Spectrum1D with flux, spectral + axis, and uncertainty. + """ + + sigma_pix = 4 + col_model = models.Gaussian1D(amplitude=1, mean=nrows/2, + stddev=sigma_pix) + spec2dvar = np.ones((nrows, ncols)) + noise = 0 + if add_noise: + np.random.seed(7) + sigma_noise = 1 + noise = np.random.normal(scale=sigma_noise, size=(nrows, ncols)) + + index_arr = np.tile(np.arange(nrows), (ncols, 1)) + img = col_model(index_arr.T) + noise + + if nan_slices: # add nans in data + for s in nan_slices: + img[s] = np.nan + + wave = np.arange(0, img.shape[1], 1) + objectspec = Spectrum1D(spectral_axis=wave*u.m, flux=img*u.Jy, + uncertainty=VarianceUncertainty(spec2dvar*u.Jy*u.Jy)) + + return objectspec + + def test_boxcar_fully_masked(): + """ + Test that the appropriate error is raised by `BoxcarExtract` when image + is fully masked/NaN. + """ + + img = mk_img() + + with pytest.raises(ValueError, match='Image is fully masked.'): + # fully NaN image + img = np.zeros((4, 5)) * np.nan + Background(img, traces=FlatTrace(self.mk_img(), 2)) + + with pytest.raises(ValueError, match='Image is fully masked.'): + # fully masked image (should be equivalent) + img = NDData(np.ones((4, 5)), mask=np.ones((4, 5))) + Background(img, traces=FlatTrace(self.mk_img(), 2)) + + # Now test that an image that isn't fully masked, but is fully masked + # within the window determined by `width`, produces the correct result + msg = 'Image is fully masked within background window determined by `width`.' + with pytest.raises(ValueError, match=msg): + img = self.mk_img(nrows=12, ncols=12, nan_slices=[np.s_[3:10, :]]) + Background(img, traces=FlatTrace(img, 6), width=7) \ No newline at end of file diff --git a/specreduce/tests/test_tracing.py b/specreduce/tests/test_tracing.py index 2e1af8b..03cc294 100644 --- a/specreduce/tests/test_tracing.py +++ b/specreduce/tests/test_tracing.py @@ -9,6 +9,34 @@ IM = make_2d_trace_image() +def mk_img(nrows=200, ncols=160, nan_slices=None, add_noise=True): + + """ + Makes a gaussian image for testing, with optional added gaussian + nosie and optional data values set to NaN. + """ + + # NOTE: Will move this to a fixture at some point. + + sigma_pix = 4 + col_model = models.Gaussian1D(amplitude=1, mean=nrows/2, + stddev=sigma_pix) + noise = 0 + if add_noise: + np.random.seed(7) + sigma_noise = 1 + noise = np.random.normal(scale=sigma_noise, size=(nrows, ncols)) + + index_arr = np.tile(np.arange(nrows), (ncols, 1)) + img = col_model(index_arr.T) + noise + + if nan_slices: # add nans in data + for s in nan_slices: + img[s] = np.nan + + return img * u.DN + + # test basic trace class def test_basic_trace(): t_pos = IM.shape[0] / 2 @@ -42,6 +70,8 @@ def test_negative_flat_trace_err(): # negative trace_pos with pytest.raises(ValueError, match='must be positive.'): FlatTrace(IM, trace_pos=-1) + with pytest.raises(ValueError, match='must be positive.'): + FlatTrace(IM, trace_pos=0) # test array traces @@ -73,8 +103,8 @@ def test_array_trace(): # make sure nonfinite data in input `trace` is masked arr[0:5] = np.nan t = ArrayTrace(IM, arr) - assert np.all(t.trace.mask[0:5] == True) - assert np.all(t.trace.mask[5:] == False) + assert np.all(t.trace.mask[0:5]) + assert np.all(t.trace.mask[5:] == 0) # test fitted traces @@ -152,46 +182,100 @@ def test_fit_trace(): @pytest.mark.filterwarnings("ignore:The fit may be unsuccessful") @pytest.mark.filterwarnings("ignore:Model is linear in parameters") class TestMasksTracing(): - - def mk_img(self, nrows=200, ncols=160, nan_slices=None, add_noise=True): - + """ + There are three currently implemented options for masking in FitTrace: filter, + omit, and zero-fill. Trace, FlatTrace, and ArrayTrace do not have + `mask_treatment` options as input because masked/nonfinite values in the data + are not relevant for those trace types as they are not affected by masked + input data. The tests in this class test masking options for FitTrace, as + well as some basic tests (errors, etc) for the other trace types. + """ + + def test_flat_and_basic_trace_mask(self): """ - Makes a gaussian image for testing, with optional added gaussian - nosie and optional data values set to NaN. + Mask handling is not relevant for basic and flat trace - nans or masked + values in the input image will not impact the trace value. The attribute + should be initialized though, and be one of the valid options ([None] + in this case), for consistancy with all other Specreduce operations. + Note that unlike FitTrace, a fully-masked image should NOT result in an + error raised because the trace does not depend on the data. """ - sigma_pix = 4 - col_model = models.Gaussian1D(amplitude=1, mean=nrows/2, - stddev=sigma_pix) - noise = 0 - if add_noise: - np.random.seed(7) - sigma_noise = 1 - noise = np.random.normal(scale=sigma_noise, size=(nrows, ncols)) + img = mk_img(nrows=5, ncols=5) - index_arr = np.tile(np.arange(nrows), (ncols, 1)) - img = col_model(index_arr.T) + noise + basic_trace = Trace(img) + assert basic_trace.mask_treatment is None - if nan_slices: # add nans in data - for s in nan_slices: - img[s] = np.nan + flat_trace = FlatTrace(img, trace_pos=2) + assert flat_trace.mask_treatment is None - return img * u.DN + arr = [1, 2, np.nan, 3, 4] + array_trace = ArrayTrace(img, arr) + assert array_trace.mask_treatment is None - def test_window_fit_trace(self): + def test_array_trace_masking(self): + """ + The `trace` input to ArrayTrace can be a masked array, or an array + containing nonfinite data which will be converted to a masked array. + Additionally, if any padding needs to be added, the returned trace will + be a masked array. Otherwise, it should be a regular array. + + Even though an ArrayTrace may have nans or masked values + in the input 1D array for the trace, `mask_treatment_method` refers + to how masked values in the input data should be treated. Nans / masked + values passed in the array trace should be considered intentional, so + also test that `mask_treatment` is initialized to None. """ - This test function will test that masked values are treated correctly in - FitTrace, and produce the correct results and warning messages based on - `peak_method`. + img = mk_img(nrows=10, ncols=10) + + # non-finite data in input trace should be masked out + trace_arr = np.array((1, 2, np.nan, 4, 5)) + array_trace = ArrayTrace(img, trace_arr) + assert array_trace.trace.mask[2] + assert isinstance(array_trace.trace, np.ma.MaskedArray) + + # and combined with any input masked data + trace_arr = np.ma.MaskedArray([1, 2, np.nan, 4, 5], mask=[1, 0, 0, 0, 0]) + array_trace = ArrayTrace(img, trace_arr) + assert array_trace.trace.mask[0] + assert array_trace.trace.mask[2] + assert isinstance(array_trace.trace, np.ma.MaskedArray) + + # check that mask_treatment is None as there are no valid choices + assert array_trace.mask_treatment is None + + # check that if array is fully finite and not masked, that the returned + # trace is a notrmal array, not a masked array + trace = ArrayTrace(img, np.ones(100)) + assert isinstance(trace.trace, np.ndarray) + assert not isinstance(trace.trace, np.ma.MaskedArray) + + # ensure correct warning is raised when entire trace is masked. + trace_arr = np.ma.MaskedArray([1, 2, np.nan, 4, 5], mask=[1, 1, 0, 1, 1]) + with pytest.raises(UserWarning, match=r'Entire trace array is masked.'): + array_trace = ArrayTrace(img, trace_arr) + + def test_fit_trace_fully_masked_image(self): + """ + Test that the correct warning is raised when a fully maksed image is + encountered. Also test that when a non-fully masked image is provided, + but `window` is set and the image is fully masked within that window, + that the correct error is raised. """ - img = self.mk_img() + # make simple gaussian image. + img = mk_img() - # create same-shaped variations of image with invalid values + # create same-shaped variations of image with nans in data array + # which will be masked within FitTrace. nrows = 200 ncols = 160 img_all_nans = np.tile(np.nan, (nrows, ncols)) + # error on trace of all-nan image + with pytest.raises(ValueError, match=r'Image is fully masked. Check for invalid values.'): + FitTrace(img_all_nans) + window = 10 guess = int(nrows/2) img_win_nans = img.copy() @@ -201,27 +285,68 @@ def test_window_fit_trace(self): with pytest.raises(ValueError, match='pixels in window region are masked'): FitTrace(img_win_nans, guess=guess, window=window) - # error on trace of all-nan image - with pytest.raises(ValueError, match=r'Image is fully masked. Check for invalid values.'): - FitTrace(img_all_nans) + def test_fit_trace_fully_masked_columns_warn_msg(self): + """ + Test that the correct warning message is raised when fully masked columns + (in a not-fully-masked image) are encountered in FitTrace. These columns + will be set to NaN and filtered from the final all-bin fit (as tested in + test_fit_trace_fully_masked_cols), but a warning message is raised. This + should happen for mask_treatment='filter' and 'omit' (for 'zero-fill', + all NaN columns become all-zero columns). + """ + img = mk_img() + + # test that warning (dependent on choice of `peak_method`) is raised when a + # few bins are masked, and that theyre listed individually + mask = np.zeros(img.shape) + mask[:, 100] = 1 + mask[:, 20] = 1 + mask[:, 30] = 1 + nddat = NDData(data=img, mask=mask, unit=u.DN) + + match_str = 'All pixels in bins 20, 30, 100 are fully masked. Setting bin peaks to NaN.' + + with pytest.warns(UserWarning, match=match_str): + FitTrace(nddat, peak_method='max') + + with pytest.warns(UserWarning, match=match_str): + FitTrace(nddat, peak_method='centroid') + + with pytest.warns(UserWarning, match=match_str): + FitTrace(nddat, peak_method='gaussian') + + # and when many bins are masked, that the message is consolidated + mask = np.zeros(img.shape) + mask[:, 0:21] = 1 + nddat = NDData(data=img, mask=mask, unit=u.DN) + with pytest.warns(UserWarning, match='All pixels in bins ' + '0, 1, 2, 3, 4, 5, 6, 7, 8, 9, ..., 20 are ' + 'fully masked. Setting bin peaks to NaN.'): + FitTrace(nddat) @pytest.mark.filterwarnings("ignore:All pixels in bins") @pytest.mark.parametrize("mask_treatment", ['filter', 'omit']) - def test_fit_trace_all_nan_cols(self, mask_treatment): - """ - Create a test image that has some fully-masked columns, and test that + def test_fit_trace_fully_masked_cols(self, mask_treatment): + """ + Create a test image with some fully-nan/masked columns, and test that when the final fit to all bin peaks is done for the trace, that these - fully-masked columns are set to NaN and filtered during the fit. This - should happen for mask_treatment = 'filter' and 'omit' (for 'zero-fill', - all NaN columns become all-zero columns). + fully-masked columns are set to NaN and filtered during the final all-bin + fit. This should happen for mask_treatment = 'filter' and 'omit' + (for 'zero-fill', all NaN columns become all-zero columns). Ignore the + warning that is produced when this case is encountered (that is tested + in `test_fit_trace_fully_masked_cols_warn_msg`.) """ - img = self.mk_img(nrows=10, ncols=11) + img = mk_img(nrows=10, ncols=11) + # set some columns fully to nan, which will be masked out img[:, 7] = np.nan img[:, 4] = np.nan img[:, 0] = np.nan - # peak_method = 'max' + # also create an image that doesn't have nans in the data, but + # is masked in the same locations, to make sure that is equivilant. + + # test peak_method = 'max' truth = [1.6346154, 2.2371795, 2.8397436, 3.4423077, 4.0448718, 4.6474359, 5.25, 5.8525641, 6.4551282, 7.0576923, 7.6602564] @@ -244,55 +369,16 @@ def test_fit_trace_all_nan_cols(self, mask_treatment): max_trace = FitTrace(img, peak_method='centroid', mask_treatment=mask_treatment) np.testing.assert_allclose(truth, max_trace.trace) - - def test_warn_msg_fit_trace_all_nan_cols(self): - """ - Test that the correct warning is raised when fully masked columns - are encountered in FitTrace. These columns will be set to NaN and filtered - from the final all-bin fit. This should happen for - mask_treatment='filter' and 'omit' (for 'zero-fill', all NaN columns - become all-zero columns). - """ - img = self.mk_img() - - # test that warning (dependent on choice of `peak_method`) is raised when a - # few bins are masked, and that theyre listed individually - mask = np.zeros(img.shape) - mask[:, 100] = 1 - mask[:, 20] = 1 - mask[:, 30] = 1 - nddat = NDData(data=img, mask=mask, unit=u.DN) - - match_str = 'All pixels in bins 20, 30, 100 are fully masked. Setting bin peaks to NaN.' - - with pytest.warns(UserWarning, match=match_str): - FitTrace(nddat, peak_method='max') - - with pytest.warns(UserWarning, match=match_str): - FitTrace(nddat, peak_method='centroid') - - with pytest.warns(UserWarning, match=match_str): - FitTrace(nddat, peak_method='gaussian') - - # and when many bins are masked, that the message is consolidated - mask = np.zeros(img.shape) - mask[:, 0:21] = 1 - nddat = NDData(data=img, mask=mask, unit=u.DN) - with pytest.warns(UserWarning, match='All pixels in bins ' - '0, 1, 2, 3, 4, 5, 6, 7, 8, 9, ..., 20 are ' - 'fully masked. Setting bin peaks to NaN.'): - FitTrace(nddat) - @pytest.mark.filterwarnings("ignore:All pixels in bins") @pytest.mark.parametrize("peak_method,expected", - [("max", [ 5., 3., 5., 5., 7., 5., - np.nan, 5., 5., 5., 2., 5.]), - ("gaussian", [5., 2.10936004, 5., 5., 7.80744334, - 5., np.nan, 5., 5., 5., 1.28216332, 5.]), - ("centroid", [4.27108332, 2.24060342, 4.27108332, - 4.27108332, 6.66827608, 4.27108332, - np.nan, 4.27108332, 4.27108332, - 4.27108332, 1.19673467, 4.27108332])]) + [("max", [5., 3., 5., 5., 7., 5., + np.nan, 5., 5., 5., 2., 5.]), + ("gaussian", [5., 2.10936004, 5., 5., 7.80744334, + 5., np.nan, 5., 5., 5., 1.28216332, 5.]), + ("centroid", [4.27108332, 2.24060342, 4.27108332, + 4.27108332, 6.66827608, 4.27108332, + np.nan, 4.27108332, 4.27108332, + 4.27108332, 1.19673467, 4.27108332])]) def test_mask_treatment_filter(self, peak_method, expected): """ Test for mask_treatment=filter for FitTrace. @@ -303,14 +389,14 @@ def test_mask_treatment_filter(self, peak_method, expected): """ # Make an image with some nonfinite values. - image1 = self.mk_img(nan_slices=[np.s_[4:8, 1:2], np.s_[2:7, 4:5], - np.s_[:, 6:7], np.s_[3:9, 10:11]], - nrows=10, ncols=12, add_noise = False) + image1 = mk_img(nan_slices=[np.s_[4:8, 1:2], np.s_[2:7, 4:5], + np.s_[:, 6:7], np.s_[3:9, 10:11]], + nrows=10, ncols=12, add_noise=False) # Also make an image that doesn't have nonf data values, but has masked # values at the same locations, to make sure they give the same results. mask = ~np.isfinite(image1) - dat = self.mk_img(nrows=10, ncols=12, add_noise=False) + dat = mk_img(nrows=10, ncols=12, add_noise=False) image2 = NDData(dat, mask=mask) for imgg in [image1, image2]: @@ -327,19 +413,19 @@ def test_mask_treatment_filter(self, peak_method, expected): mask = np.isfinite(y_bins) all_bin_fit = fitter(trace.trace_model, x_bins[mask], y_bins[mask]) all_bin_fit = all_bin_fit((np.arange(12))) - + np.testing.assert_allclose(trace.trace, all_bin_fit) @pytest.mark.filterwarnings("ignore:All pixels in bins") @pytest.mark.parametrize("peak_method,expected", - [("max", [ 5., 3., 5., 5., 7., 5., - 0., 5., 5., 5., 2., 5.]), - ("gaussian", [5., 1.09382384, 5., 5., 7.81282206, - 5., 0., 5., 5., 5., 1.28216332, 5.]), - ("centroid", [4.27108332, 2.24060342, 4.27108332, - 4.27108332, 6.66827608, 4.27108332, - 9., 4.27108332, 4.27108332, - 4.27108332, 1.19673467, 4.27108332])]) + [("max", [5., 3., 5., 5., 7., 5., + 0., 5., 5., 5., 2., 5.]), + ("gaussian", [5., 1.09382384, 5., 5., 7.81282206, + 5., 0., 5., 5., 5., 1.28216332, 5.]), + ("centroid", [4.27108332, 2.24060342, 4.27108332, + 4.27108332, 6.66827608, 4.27108332, + 9., 4.27108332, 4.27108332, + 4.27108332, 1.19673467, 4.27108332])]) def test_mask_treatment_zero_fill(self, peak_method, expected): """ Test for mask_treatment=`zero_fill` for FitTrace. @@ -349,14 +435,14 @@ def test_mask_treatment_zero_fill(self, peak_method, expected): """ # Make an image with some nonfinite values. - image1 = self.mk_img(nan_slices=[np.s_[4:8, 1:2], np.s_[2:7, 4:5], - np.s_[:, 6:7], np.s_[3:9, 10:11]], - nrows=10, ncols=12, add_noise = False) + image1 = mk_img(nan_slices=[np.s_[4:8, 1:2], np.s_[2:7, 4:5], + np.s_[:, 6:7], np.s_[3:9, 10:11]], + nrows=10, ncols=12, add_noise=False) # Also make an image that doesn't have nonf data values, but has masked # values at the same locations, to make sure they give the same results. mask = ~np.isfinite(image1) - dat = self.mk_img(nrows=10, ncols=12, add_noise=False) + dat = mk_img(nrows=10, ncols=12, add_noise=False) image2 = NDData(dat, mask=mask) for imgg in [image1, image2]: @@ -374,35 +460,36 @@ def test_mask_treatment_zero_fill(self, peak_method, expected): mask = np.isfinite(y_bins) all_bin_fit = fitter(trace.trace_model, x_bins[mask], y_bins[mask]) all_bin_fit = all_bin_fit((np.arange(12))) - + np.testing.assert_allclose(trace.trace, all_bin_fit) @pytest.mark.filterwarnings("ignore:All pixels in bins") @pytest.mark.parametrize("peak_method,expected", - [("max", [ 5., np.nan, 5., 5., np.nan, 5., - np.nan, 5., 5., 5., np.nan, 5.]), - ("gaussian", [ 5., np.nan, 5., 5., np.nan, 5., - np.nan, 5., 5., 5., np.nan, 5.]), - ("centroid", [4.27108332, np.nan, 4.27108332, - 4.27108332, np.nan, 4.27108332, - np.nan, 4.27108332, 4.27108332, - 4.27108332, np.nan, 4.27108332])]) + [("max", [5., np.nan, 5., 5., np.nan, 5., + np.nan, 5., 5., 5., np.nan, 5.]), + ("gaussian", [5., np.nan, 5., 5., np.nan, 5., + np.nan, 5., 5., 5., np.nan, 5.]), + ("centroid", [4.27108332, np.nan, 4.27108332, + 4.27108332, np.nan, 4.27108332, + np.nan, 4.27108332, 4.27108332, + 4.27108332, np.nan, 4.27108332])]) def test_mask_treatment_omit(self, peak_method, expected): """ Test for mask_treatment=`omit` for FitTrace. Columns (assuming - disp_axis==1) with any masked data values will be fully masked and + disp_axis==1) with any masked data values will be fully masked and therefore not contribute to the bin peaks. Parametrized over different - `peak_method` options.""" + `peak_method` options. + """ # Make an image with some nonfinite values. - image1 = self.mk_img(nan_slices=[np.s_[4:8, 1:2], np.s_[2:7, 4:5], - np.s_[:, 6:7], np.s_[3:9, 10:11]], - nrows=10, ncols=12, add_noise = False) + image1 = mk_img(nan_slices=[np.s_[4:8, 1:2], np.s_[2:7, 4:5], + np.s_[:, 6:7], np.s_[3:9, 10:11]], + nrows=10, ncols=12, add_noise=False) - # Also make an image that doesn't have nonf data values, but has masked - # values at the same locations, to make sure they give the same results. + # Also make an image that doesn't have nonfinite data values, but has masked + # values at the same locations, to make sure those cases are equivalent mask = ~np.isfinite(image1) - dat = self.mk_img(nrows=10, ncols=12, add_noise=False) + dat = mk_img(nrows=10, ncols=12, add_noise=False) image2 = NDData(dat, mask=mask) for imgg in [image1, image2]: @@ -421,5 +508,5 @@ def test_mask_treatment_omit(self, peak_method, expected): mask = np.isfinite(y_bins) all_bin_fit = fitter(trace.trace_model, x_bins[mask], y_bins[mask]) all_bin_fit = all_bin_fit((np.arange(12))) - + np.testing.assert_allclose(trace.trace, all_bin_fit) diff --git a/specreduce/tracing.py b/specreduce/tracing.py index 018e62c..fdce151 100644 --- a/specreduce/tracing.py +++ b/specreduce/tracing.py @@ -37,9 +37,11 @@ def __post_init__(self): self.trace = np.ones_like(self.image[0]) * self.trace_pos # masking options not relevant for basic Trace - self.mask_treatment = None - self._valid_mask_treatment_methods = None - + self._mask_treatment = None + self._valid_mask_treatment_methods = [None] + + # eventually move this to common SpecreduceOperation base class + self.validate_masking_options() def __getitem__(self, i): return self.trace[i] @@ -48,6 +50,11 @@ def __getitem__(self, i): def shape(self): return self.trace.shape + def validate_masking_options(self): + if self.mask_treatment not in self.valid_mask_treatment_methods: + raise ValueError( + f'`mask_treatment` {self.mask_treatment} not one of {self.valid_mask_treatment_methods}') + def shift(self, delta): """ Shift the trace by delta pixels perpendicular to the axis being traced @@ -66,7 +73,7 @@ def _bound_trace(self): Mask trace positions that are outside the upper/lower bounds of the image. """ ny = self.image.shape[0] - self.trace = np.ma.masked_outside(self.trace, 0, ny-1) + self.trace = np.ma.masked_outside(self.trace, 0, ny - 1) def __add__(self, delta): """ @@ -84,6 +91,14 @@ def __sub__(self, delta): """ return self.__add__(-delta) + @property + def mask_treatment(self): + return self._mask_treatment + + @property + def valid_mask_treatment_methods(self): + return self._valid_mask_treatment_methods + @dataclass class FlatTrace(Trace, _ImageParser): @@ -106,9 +121,9 @@ def __post_init__(self): self.set_position(self.trace_pos) - # masking options not relevant for FlatTrace - self.mask_treatment = None - self._valid_mask_treatment_methods = None + # masking options not relevant for basic Trace + self._mask_treatment = None + self._valid_mask_treatment_methods = [None] def set_position(self, trace_pos): """ @@ -133,25 +148,32 @@ class ArrayTrace(Trace, _ImageParser): Parameters ---------- - trace : `numpy.ndarray` - Array containing trace positions + trace : `numpy.ndarray` or `numpy.ma.MaskedArray` + Array containing trace positions. """ trace: np.ndarray def __post_init__(self): - # masking options not relevant for FlatTrace. any non-finite or masked + # masking options not relevant for ArrayTrace. any non-finite or masked # data in `image` will not affect array trace - self.mask_treatment = None - self._valid_mask_treatment_methods = None + self._mask_treatment = None + self._valid_mask_treatment_methods = [None] + + # masked array will have a .data, regular array will not. + trace_data = getattr(self.trace, 'data', self.trace) - # but we do need to mask non-finite data in input trace array. - nonf_mask = ~np.isfinite(self.trace) + # but we do need to mask uncaught non-finite values in input trace array + # which should also be combined with any existing mask in the input `trace` + if hasattr(self.trace, 'mask'): + total_mask = np.logical_or(self.trace.mask, ~np.isfinite(trace_data)) + else: + total_mask = ~np.isfinite(trace_data) - # if there are any non-finite values in input `trace` array, or later on - # any padded values then trace.trace will have masked values - if np.any(nonf_mask): - self.trace = np.ma.MaskedArray(self.trace, nonf_mask) + # always work with masked array, even if there is no masked + # or nonfinite data, in case padding is needed. if not, mask will be + # dropped at the end and a regular array will be returned. + self.trace = np.ma.MaskedArray(trace_data, total_mask) self.image = self._parse_image(self.image) @@ -169,6 +191,14 @@ def __post_init__(self): self._bound_trace() + # warn if entire trace is masked + if np.all(self.trace.mask): + warnings.warn("Entire trace array is masked.") + + # and return plain array if nothing is masked + if not np.any(self.trace.mask): + self.trace = self.trace.data + @dataclass class FitTrace(Trace, _ImageParser): @@ -248,10 +278,10 @@ class FitTrace(Trace, _ImageParser): peak_method: str = 'max' _crossdisp_axis = 0 _disp_axis = 1 - mask_treatment : str = 'filter' + mask_treatment: str = 'filter' _valid_mask_treatment_methods = ('filter', 'omit', 'zero-fill') # for testing purposes only, save bin peaks if requested - _save_bin_peaks_testing : bool = False + _save_bin_peaks_testing: bool = False def __post_init__(self): @@ -360,7 +390,7 @@ def _fit_trace(self, img): # binned columns, summed along disp. axis. # or just a single, unbinned column if no bins - z_i = img[ilum2, x_bins[i]:x_bins[i+1]].sum(axis=self._disp_axis) + z_i = img[ilum2, x_bins[i]:x_bins[i + 1]].sum(axis=self._disp_axis) # if this bin is fully masked, set bin peak to NaN so it can be # filtered in the final fit to all bin peaks for the trace @@ -402,7 +432,7 @@ def _fit_trace(self, img): z_i_cumsum = np.cumsum(z_i) # find the interpolated index where the cumulative array reaches # half the total cumulative values - y_bins[i] = np.interp(z_i_cumsum[-1]/2., z_i_cumsum, ilum2) + y_bins[i] = np.interp(z_i_cumsum[-1] / 2., z_i_cumsum, ilum2) # NOTE this reflects current behavior, should eventually be changed # to set to nan by default (or zero fill / interpoate option once @@ -434,7 +464,7 @@ def _fit_trace(self, img): # for testing purposes only, save bin peaks if requested if self._save_bin_peaks_testing: - self._bin_peaks_testing = (x_bins, y_bins) + self._bin_peaks_testing = (x_bins, y_bins) # filter non-finite bin peaks before filtering to all bin peaks y_finite = np.where(np.isfinite(y_bins))[0]