Skip to content

respect order_start/order_end in extraction #35

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
leejjoon wants to merge 1 commit into
base: master
Choose a base branch
from
Open

respect order_start/order_end in extraction #35

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion igrins/igrins_recipes/recipe_extract.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
),
Expand Down
27 changes: 25 additions & 2 deletions igrins/procedures/apertures.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)):
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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]

Expand Down Expand Up @@ -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)

Expand Down
6 changes: 2 additions & 4 deletions igrins/procedures/slit_profile.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand All @@ -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,
Expand Down
94 changes: 71 additions & 23 deletions igrins/procedures/target_spec.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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"))
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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.
Expand Down