From 9762b3f02a7fb2401a0c4804f41f9bd9e86d4733 Mon Sep 17 00:00:00 2001 From: Jae-Joon Lee Date: Mon, 13 May 2024 21:37:42 +0900 Subject: [PATCH] respect order_start/order_end in extraction --- igrins/igrins_recipes/recipe_extract.py | 2 +- igrins/procedures/apertures.py | 27 ++++++- igrins/procedures/slit_profile.py | 6 +- igrins/procedures/target_spec.py | 94 +++++++++++++++++++------ 4 files changed, 99 insertions(+), 30 deletions(-) diff --git a/igrins/igrins_recipes/recipe_extract.py b/igrins/igrins_recipes/recipe_extract.py index 2a66500..5bc1325 100644 --- a/igrins/igrins_recipes/recipe_extract.py +++ b/igrins/igrins_recipes/recipe_extract.py @@ -91,7 +91,7 @@ def estimate_slit_profile_extended(obsset, Step("Setup extraction parameters", setup_extraction_parameters, height_2dspec=0, - order_range="-1,-1", + order_range="", correct_flexure=False, mask_cosmics=False, ), diff --git a/igrins/procedures/apertures.py b/igrins/procedures/apertures.py index 5167af3..dc2763c 100644 --- a/igrins/procedures/apertures.py +++ b/igrins/procedures/apertures.py @@ -29,7 +29,7 @@ def set_order_minmax_to_extract(self, start_order, end_order): if (end_order is None) or (end_order < 0): end_order = self.orders[-1] - self.orders_to_extract = range(start_order, end_order+1) + self.orders_to_extract = list(range(start_order, end_order+1)) def __init__(self, orders, bottomup_solutions, basename="", order_minmax_to_extract=(-1, -1)): @@ -176,6 +176,9 @@ def extract_spectra_v2(self, data, f1=0., f2=1.): return s_list def extract_spectra_from_ordermap(self, data, order_map): + # FIXME please check is orders_to_extract is propertly respected and + # then remove the raise statement + raise RuntimeError("This method is not tested.") slices = ni.find_objects(order_map) s_list = [] for o in self.orders_to_extract: @@ -239,6 +242,15 @@ def extract_simple_from_shifted(self, ordermap, slices = ni.find_objects(ordermap) for o in self.orders_to_extract: + if o > len(slices) or slices[o-1] is None: + s = np.empty(2048, dtype=float) + v = np.empty(2048, dtype=float) + s.fill(np.nan) + v.fill(np.nan) + s_list.append(s) + v_list.append(v) + continue + sl = slices[o-1][0], slice(0, 2048) msk = (ordermap[sl] == o) @@ -284,6 +296,15 @@ def extract_stellar_from_shifted(self, ordermap_bpixed, # hl.append(pyfits.PrimaryHDU()) for o in self.orders_to_extract: + if o > len(slices) or slices[o-1] is None: + s = np.empty(2048, dtype=float) + v = np.empty(2048, dtype=float) + s.fill(np.nan) + v.fill(np.nan) + s_list.append(s) + v_list.append(v) + continue + sl = slices[o-1][0], slice(0, 2048) msk = (ordermap_bpixed[sl] == o) & msk1[sl] @@ -627,7 +648,9 @@ def make_synth_map(self, order_map, slitpos_map, xx = np.arange(2048) slices = ni.find_objects(order_map) - for o, s in zip(self.orders, s_list): + for o, s in zip(self.orders_to_extract, s_list): + if o > len(slices) or slices[o-1] is None: + continue sl = slices[o-1][0], slice(0, 2048) msk = (order_map[sl] == o) diff --git a/igrins/procedures/slit_profile.py b/igrins/procedures/slit_profile.py index 4022e28..b3567f2 100644 --- a/igrins/procedures/slit_profile.py +++ b/igrins/procedures/slit_profile.py @@ -124,9 +124,7 @@ def profile(order, x_pixel, y_slit_pos): postfix=obsset.basename_postfix).data data_minus_flattened = data_minus / orderflat - from .aperture_helper import get_aperture_from_obsset - orders = helper.get("orders") - ap = get_aperture_from_obsset(obsset, orders=orders) + ap = helper.get_aperture(obsset) ordermap = helper.get("ordermap") ordermap_bpixed = helper.get("ordermap_bpixed") @@ -145,7 +143,7 @@ def profile(order, x_pixel, y_slit_pos): profile_x, profile_y = _get_norm_profile(bins, hh0) # profile = get_profile_func(profile_x, profile_y) - slit_profile_dict = dict(orders=orders, + slit_profile_dict = dict(orders=ap.orders_to_extract, ab_mode=do_ab, slit_profile_list=slit_profile_list, profile_x=profile_x, diff --git a/igrins/procedures/target_spec.py b/igrins/procedures/target_spec.py index 3a5ff0c..07b3fb1 100644 --- a/igrins/procedures/target_spec.py +++ b/igrins/procedures/target_spec.py @@ -15,18 +15,19 @@ def _get_int_from_config(obsset, kind, default): return v -def setup_extraction_parameters(obsset, order_range="-1,-1", +def setup_extraction_parameters(obsset, order_range="", height_2dspec=0, correct_flexure=False, mask_cosmics=False): - _order_range_s = order_range - try: - order_start, order_end = map(int, _order_range_s.split(",")) - except Exception: - msg = "Failed to parse order range: {}".format(_order_range_s) - raise ValueError(msg) - - order_start = _get_int_from_config(obsset, "ORDER_START", order_start) - order_end = _get_int_from_config(obsset, "ORDER_END", order_end) + if order_range: + _order_range_s = order_range + try: + order_start, order_end = map(int, _order_range_s.split(",")) + except Exception: + msg = "Failed to parse order range: {}".format(_order_range_s) + raise ValueError(msg) + else: + order_start = _get_int_from_config(obsset, "ORDER_START", -1) + order_end = _get_int_from_config(obsset, "ORDER_END", -1) height_2dspec = _get_int_from_config(obsset, "HEIGHT_2DSPEC", height_2dspec) @@ -271,11 +272,42 @@ def convert_data(d): return header.copy(), hdu.data, convert_data +def _reorder_data_to_orders_to_extract(ap, d0): + d_map = dict(zip(ap.orders, d0)) + d0 = np.asarray(d0) + d = np.empty((len(ap.orders_to_extract),) + d0.shape[1:], dtype=d0.dtype) + d.fill(np.nan) + + for i, o in enumerate(ap.orders_to_extract): + _d = d_map.get(o, None) + if _d is None: + continue + d[i] = _d + + return d + +from astropy.io.fits.header import Header +def remove_wat(header): + cards = [card for card in header.cards if not card.keyword.startswith("WAT")] + header = Header(cards) + return header + def store_1dspec(obsset, v_list, s_list, sn_list=None): basename_postfix = obsset.basename_postfix - wvl_header, wvl_data, convert_data = get_wvl_header_data(obsset) + wvl_header, wvl_data0, convert_data = get_wvl_header_data(obsset) + + if (obsset.get_recipe_parameter("order_start"), + obsset.get_recipe_parameter("order_end")) != (-1, -1): + # If custom extraction orders are used, just drop WAT headers for simplicity. + # FIXME We should rewrite the WAT header instead of removing it. + wvl_header = remove_wat(wvl_header) + + helper = ResourceHelper(obsset) + ap = helper.get("aperture") + + wvl_data = _reorder_data_to_orders_to_extract(ap, wvl_data0) d = np.array(v_list) v_data = convert_data(d.astype("float32")) @@ -365,7 +397,12 @@ def store_2dspec(obsset, d0_shft_list, msk_shft_list = _ with np.errstate(invalid="ignore"): - d = np.array(d0_shft_list) / np.array(msk_shft_list) + d0 = np.array(d0_shft_list) / np.array(msk_shft_list) + + helper = ResourceHelper(obsset) + ap = helper.get("aperture") + + d = _reorder_data_to_orders_to_extract(ap, d0) hdul = obsset.get_hdul_to_write(([], convert_data(d.astype("float32")))) # wvl_header.update(hdul[0].header) @@ -383,7 +420,9 @@ def store_2dspec(obsset, d0_shft_list, msk_shft_list = _ with np.errstate(invalid="ignore"): - d = np.array(d0_shft_list) / np.array(msk_shft_list) + d0 = np.array(d0_shft_list) / np.array(msk_shft_list) + + d = _reorder_data_to_orders_to_extract(ap, d0) hdul = obsset.get_hdul_to_write(([], convert_data(d.astype("float32")))) # wvl_header.update(hdul[0].header) @@ -444,17 +483,21 @@ def extract_stellar_spec(obsset, extraction_mode="optimal", s_list, v_list, cr_mask, aux_images = _ # calculate S/N per resolution - wvl_solutions = helper.get("wvl_solutions") + wvl_solutions0 = helper.get("wvl_solutions") + wvl_solutions = _reorder_data_to_orders_to_extract(ap, wvl_solutions0) sn_list = [] for wvl, s, v in zip(wvl_solutions, s_list, v_list): - if pixel_per_res_element is None: - dw = np.gradient(wvl) - _pixel_per_res_element = (wvl/40000.)/dw + if wvl is None: + _pixel_per_res_element = np.nan else: - _pixel_per_res_element = float(pixel_per_res_element) + if pixel_per_res_element is None: + dw = np.gradient(wvl) + _pixel_per_res_element = (wvl/40000.)/dw + else: + _pixel_per_res_element = float(pixel_per_res_element) # print pixel_per_res_element[1024] # len(pixel_per_res_element) = 2047. But we ignore it. @@ -648,17 +691,22 @@ def extract_extended_spec(obsset, # calculate S/N per resolution helper = ResourceHelper(obsset) - wvl_solutions = helper.get("wvl_solutions") + ap = helper.get("aperture") + wvl_solutions0 = helper.get("wvl_solutions") + wvl_solutions = _reorder_data_to_orders_to_extract(ap, wvl_solutions0) sn_list = [] for wvl, s, v in zip(wvl_solutions, s_list, v_list): - if pixel_per_res_element is None: - dw = np.gradient(wvl) - _pixel_per_res_element = (wvl/40000.)/dw + if wvl is None: + _pixel_per_res_element = np.nan else: - _pixel_per_res_element = float(pixel_per_res_element) + if pixel_per_res_element is None: + dw = np.gradient(wvl) + _pixel_per_res_element = (wvl/40000.)/dw + else: + _pixel_per_res_element = float(pixel_per_res_element) # print pixel_per_res_element[1024] # len(pixel_per_res_element) = 2047. But we ignore it.