From 4221096588fdc0ed9a1177e72734bf3ddbe1ddae Mon Sep 17 00:00:00 2001 From: Curtis McCully Date: Tue, 12 Nov 2024 08:35:41 -0500 Subject: [PATCH 1/3] Fixes based on commissioning. We now trim the edges of the orders better where the background fitting algorithm introduces artifacts. --- banzai_floyds/extract.py | 35 +++++++++++++++++------------ banzai_floyds/orders.py | 2 +- banzai_floyds/settings.py | 2 +- banzai_floyds/tests/test_extract.py | 28 ++++++++++++++++++++++- banzai_floyds/tests/utils.py | 4 +++- banzai_floyds/trim.py | 32 +++++++++++++++++++------- 6 files changed, 77 insertions(+), 26 deletions(-) diff --git a/banzai_floyds/extract.py b/banzai_floyds/extract.py index 74c60ac..1334ccf 100644 --- a/banzai_floyds/extract.py +++ b/banzai_floyds/extract.py @@ -57,6 +57,9 @@ def extract(binned_data, bin_key='order_wavelength_bin', data_keyword='data', ba if np.max(data_to_sum[weights_key][data_to_sum['extraction_window']]) < 5e-3: continue + data_to_sum = data_to_sum[data_to_sum['mask'] == 0] + if len(data_to_sum) == 0: + continue wavelength_bin_width = data_to_sum[bin_key + '_width'][0] # This should be equivalent to Horne 1986 optimal extraction flux = data_to_sum[data_keyword] - data_to_sum[background_key] @@ -110,22 +113,26 @@ def do_stage(self, image): # Scale the background in the same way we scaled the data so we can still subtract it cleanly image.binned_data['flux_background'] = image.binned_data['background'] * image.binned_data['flux'] image.binned_data['flux_background'] /= image.binned_data['data'] - overlap_region = [max([domain[0] for domain in image.wavelengths.domains]), - min([domain[1] for domain in image.wavelengths.domains])] - order_2 = image.binned_data['order'] == 2 - order_1 = image.binned_data['order'] == 1 - in_overlap = np.logical_and(image.binned_data['wavelength'] > overlap_region[0], - image.binned_data['wavelength'] < overlap_region[1]) - - # Since we are now combining two orders, we need to divide the weights by 2 in the overlap region - image.binned_data['combined_weights'] = image.binned_data['weights'] - image.binned_data['combined_weights'][in_overlap] /= 2 + overlap_region = [max([domain[0] for domain in image.wavelengths.wavelength_domains]), + min([domain[1] for domain in image.wavelengths.wavelength_domains])] # Normalize the orders to make sure they overlap - normalization = np.median(image.binned_data['flux'][np.logical_and(order_1, in_overlap)]) - normalization /= np.median(image.binned_data['flux'][np.logical_and(order_2, in_overlap)]) + # order_1 onto order_2 because order_1 has higher resolution + extracted_order1 = image.extracted['order'] == 1 + extracted_order2 = image.extracted['order'] == 2 + in_extracted_overlap = np.logical_and(image.extracted['wavelength'] > overlap_region[0], + image.extracted['wavelength'] < overlap_region[1]) + waves_to_interp = image.extracted['wavelength'][np.logical_and(extracted_order1, in_extracted_overlap)] + overlap_order2 = np.logical_and(extracted_order2, in_extracted_overlap) + flux_to_ratio = np.interp(waves_to_interp, image.extracted['wavelength'][overlap_order2], + image.extracted['flux'][overlap_order2]) + order_ratio = image.extracted['flux'][np.logical_and(extracted_order1, in_extracted_overlap)] + order_ratio /= flux_to_ratio + normalization = np.median(order_ratio) + order_2 = image.binned_data['order'] == 2 for key in ['flux', 'fluxerror', 'flux_background']: image.binned_data[key][order_2] *= normalization image.spectrum = extract(image.binned_data, data_keyword='flux', bin_key='wavelength_bin', - background_key='flux_background', background_out_key='background', - uncertainty_key='fluxerror', weights_key='combined_weights', include_order=False) + background_key='flux_background', background_out_key='background', flux_keyword='flux', + flux_error_key='fluxerror', uncertainty_key='fluxerror', + weights_key='weights', include_order=False) return image diff --git a/banzai_floyds/orders.py b/banzai_floyds/orders.py index 31ea8d4..61ac865 100644 --- a/banzai_floyds/orders.py +++ b/banzai_floyds/orders.py @@ -406,7 +406,7 @@ class OrderSolver(Stage): POLYNOMIAL_ORDER = 3 ORDER_REGIONS = {'ogg': [(0, 1550), (500, 1835)], - 'coj': [(0, 1600), (615, 1920)]} + 'coj': [(55, 1600), (615, 1920)]} def do_stage(self, image): if image.orders is None: diff --git a/banzai_floyds/settings.py b/banzai_floyds/settings.py index e7f767c..5026990 100644 --- a/banzai_floyds/settings.py +++ b/banzai_floyds/settings.py @@ -16,7 +16,7 @@ 'banzai_floyds.profile.ProfileFitter', 'banzai_floyds.background.BackgroundFitter', 'banzai_floyds.extract.Extractor', - # 'banzai_floyds.trim.Trimmer', + 'banzai_floyds.trim.Trimmer', 'banzai_floyds.flux.StandardLoader', 'banzai_floyds.flux.FluxSensitivity', 'banzai_floyds.flux.FluxCalibrator', diff --git a/banzai_floyds/tests/test_extract.py b/banzai_floyds/tests/test_extract.py index 0bc5144..dd59285 100644 --- a/banzai_floyds/tests/test_extract.py +++ b/banzai_floyds/tests/test_extract.py @@ -1,7 +1,7 @@ import numpy as np from banzai import context from banzai_floyds.tests.utils import generate_fake_science_frame -from banzai_floyds.extract import Extractor, extract, set_extraction_region +from banzai_floyds.extract import Extractor, extract, set_extraction_region, CombinedExtractor from banzai_floyds.utils.binning_utils import bin_data from collections import namedtuple from astropy.table import Table @@ -68,3 +68,29 @@ def test_full_extraction_stage(): residuals = frame['EXTRACTED'].data['fluxraw'] - expected residuals /= frame['EXTRACTED'].data['fluxrawerr'] assert (np.abs(residuals) < 3).sum() > 0.99 * len(frame['EXTRACTED'].data['fluxraw']) + + +def test_combined_extraction(): + np.random.seed(125325) + input_context = context.Context({}) + frame = generate_fake_science_frame(flat_spectrum=False, include_sky=True) + frame.binned_data = bin_data(frame.data, frame.uncertainty, frame.wavelengths, frame.orders) + fake_profile_width_funcs = [Legendre(frame.input_profile_sigma,) for _ in frame.input_profile_centers] + frame.profile = frame.input_profile_centers, fake_profile_width_funcs, None + frame.binned_data['background'] = frame.input_sky[frame.binned_data['y'].astype(int), + frame.binned_data['x'].astype(int)] + frame.extraction_windows = [[-5.0, 5.0], [-5.0, 5.0]] + set_extraction_region(frame) + frame.sensitivity = Table({'wavelength': [0, 1e6, 0, 1e6], 'sensitivity': [1, 1, 1, 1], 'order': [1, 1, 2, 2]}) + frame.telluric = Table({'wavelength': [0, 1e6], 'telluric': [1, 1]}) + extracted_waves = np.arange(3000.0, 10000.0) + flux = np.ones(len(extracted_waves) * 2) + orders = np.hstack([np.ones(len(extracted_waves)), np.ones(len(extracted_waves)) * 2]) + frame.extracted = Table({'wavelength': np.hstack([extracted_waves, extracted_waves]), 'flux': flux, + 'order': orders}) + stage = CombinedExtractor(input_context) + frame = stage.do_stage(frame) + expected = np.interp(frame['SPECTRUM'].data['wavelength'], frame.input_spectrum_wavelengths, frame.input_spectrum) + residuals = frame['SPECTRUM'].data['flux'] - expected + residuals /= frame['SPECTRUM'].data['fluxerror'] + assert (np.abs(residuals) < 3).sum() > 0.99 * len(frame['SPECTRUM'].data['flux']) diff --git a/banzai_floyds/tests/utils.py b/banzai_floyds/tests/utils.py index 2870d35..59d98d0 100644 --- a/banzai_floyds/tests/utils.py +++ b/banzai_floyds/tests/utils.py @@ -168,7 +168,9 @@ def generate_fake_science_frame(include_sky=False, flat_spectrum=True, fringe=Fa frame = FLOYDSObservationFrame([CCDData(data, fits.Header({'DAY-OBS': '20230101', - 'DATE-OBS': '2023-01-01 12:41:56.11'}), + 'DATE-OBS': '2023-01-01 12:41:56.11', + 'HEIGHT': 0, + 'AIRMASS': 1.0}), uncertainty=errors)], 'foo.fits') frame.input_profile_centers = profile_centers diff --git a/banzai_floyds/trim.py b/banzai_floyds/trim.py index 03c3cc7..0f063ef 100644 --- a/banzai_floyds/trim.py +++ b/banzai_floyds/trim.py @@ -4,13 +4,29 @@ class Trimmer(Stage): def do_stage(self, image): - image.extracted = image.extracted[image.extracted['wavelength'] < 10500.0] - order_1 = image.extracted['order'] == 1 - # We cut off 50 points here in the red order due to chip edge effects - # 50 is a little arbirary, but works for early tests on standard star observations + # We need to trim ~10 pixels off each order to avoid artifacts with the background + # The tilt of the arc lines is abou 8 degrees, so this corresponds to about just cutting wavelengths + # The don't fully fall on the chip + cuts = {} + for order in [1, 2]: + order_wavelengths = image.extracted['wavelength'][image.extracted['order'] == order] + order_wavelengths = order_wavelengths[np.argsort(order_wavelengths)] + cuts[order] = order_wavelengths[10], order_wavelengths[-10] + keep_extracted = np.zeros_like(image.extracted, dtype=bool) + for order, cut in cuts.items(): + passes_cut = image.extracted['order'] == order + passes_cut = np.logical_and(passes_cut, image.extracted['wavelength'] > cut[0]) + passes_cut = np.logical_and(passes_cut, image.extracted['wavelength'] < cut[1]) + keep_extracted = np.logical_or(keep_extracted, passes_cut) + + keep_extracted = np.logical_and(keep_extracted, image.extracted['wavelength'] < 10500.0) + image.extracted = image.extracted[keep_extracted] + + for order, cut in cuts.items(): + in_order = image.binned_data['order'] == order + should_mask = np.logical_or(image.binned_data['wavelength'] < cuts[order][0], + image.binned_data['wavelength'] > cuts[order][1]) + should_mask = np.logical_and(should_mask, in_order) + image.binned_data['mask'][should_mask] = 1 - sorted_wavelengths = np.argsort(image.extracted[order_1]['wavelength']) - wavelength_cut = image.extracted[order_1][sorted_wavelengths][50]['wavelength'] - image.extracted = image.extracted[np.logical_or(image.extracted['order'] == 2, - image.extracted['wavelength'] >= wavelength_cut)] return image From 126f2980ecab5723d7fb22ef817c24ee38533b12 Mon Sep 17 00:00:00 2001 From: Curtis McCully Date: Tue, 12 Nov 2024 08:40:16 -0500 Subject: [PATCH 2/3] Added changelog. --- CHANGES.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGES.md b/CHANGES.md index 2bf2a4f..e35b560 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,3 +1,8 @@ +0.11.1 (2024-11-12) +------------------- +- Fixes to the quality of the reductions +- We now trim the edges of orders better to remove artifacts + 0.11.0 (2024-11-05) ------------------- - Added the ability to combine the extraction from both orders into a single spectrum From 63b09e76771058ae261f9180b22a09f5001a7e2e Mon Sep 17 00:00:00 2001 From: Curtis McCully Date: Tue, 12 Nov 2024 14:53:45 -0500 Subject: [PATCH 3/3] Code style fix. --- banzai_floyds/tests/test_extract.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/banzai_floyds/tests/test_extract.py b/banzai_floyds/tests/test_extract.py index dd59285..025732a 100644 --- a/banzai_floyds/tests/test_extract.py +++ b/banzai_floyds/tests/test_extract.py @@ -86,7 +86,7 @@ def test_combined_extraction(): extracted_waves = np.arange(3000.0, 10000.0) flux = np.ones(len(extracted_waves) * 2) orders = np.hstack([np.ones(len(extracted_waves)), np.ones(len(extracted_waves)) * 2]) - frame.extracted = Table({'wavelength': np.hstack([extracted_waves, extracted_waves]), 'flux': flux, + frame.extracted = Table({'wavelength': np.hstack([extracted_waves, extracted_waves]), 'flux': flux, 'order': orders}) stage = CombinedExtractor(input_context) frame = stage.do_stage(frame)