diff --git a/nirc2_reduce/data/header_kw_dicts/nirc2.yaml b/nirc2_reduce/data/header_kw_dicts/nirc2.yaml index 5329670..54ab635 100644 --- a/nirc2_reduce/data/header_kw_dicts/nirc2.yaml +++ b/nirc2_reduce/data/header_kw_dicts/nirc2.yaml @@ -1,7 +1,9 @@ --- object: OBJECT instrument: INSTRUME +observer: OBSERVER date: DATE-OBS +time: EXPSTART filter: FILTER wl: EFFWAVE subc: NAXIS1 @@ -15,6 +17,10 @@ lamps: kw: FLSPECTR lampon: "on" lampoff: "off" +shutter: + kw: SHRNAME + shutopen: "open" + shutclose: "closed" warpx_file: nirc2_distort_X_post20150413_v1.fits warpy_file: nirc2_distort_Y_post20150413_v1.fits -dfits_kws: ["OBJECT", "DATE-OBS", "FILTER", "EFFWAVE", "NAXIS1", "ITIME", "COADDS", "AXESTAT", "FLSPECTR"] \ No newline at end of file +dfits_kws: ["OBJECT", "DATE-OBS", "FILTER", "EFFWAVE", "NAXIS1", "ITIME", "COADDS", "AXESTAT", "FLSPECTR", "SHRNAME", "EXPSTART", "OBSERVER"] \ No newline at end of file diff --git a/nirc2_reduce/data/header_kw_dicts/osiris.yaml b/nirc2_reduce/data/header_kw_dicts/osiris.yaml index 4766d64..d2b7dfb 100644 --- a/nirc2_reduce/data/header_kw_dicts/osiris.yaml +++ b/nirc2_reduce/data/header_kw_dicts/osiris.yaml @@ -1,7 +1,9 @@ --- object: OBJECT instrument: INSTRUME +observer: date: DATE-OBS +time: filter: IFILTER wl: EFFWAVE subc: DETENDX @@ -15,6 +17,10 @@ lamps: kw: FLSPECTR lampon: "on" lampoff: "off" +shutter: + kw: + shutopen: "" + shutclose: "" warpx_file: nirc2_distort_X_post20150413_v1.fits warpy_file: nirc2_distort_Y_post20150413_v1.fits -dfits_kws: ["OBJECT", "INSTRUME", "DATE-OBS", "IFILTER", "EFFWAVE", "NAXIS1", "TRUITIME", "COADDS", "AXESTAT", "FLSPECTR"] \ No newline at end of file +dfits_kws: ["OBJECT", "INSTRUME", "DATE-OBS", "IFILTER", "EFFWAVE", "NAXIS1", "TRUITIME", "COADDS", "AXESTAT", "FLSPECTR", ] \ No newline at end of file diff --git a/nirc2_reduce/flats.py b/nirc2_reduce/flats.py index e2ecfce..38089b3 100755 --- a/nirc2_reduce/flats.py +++ b/nirc2_reduce/flats.py @@ -30,10 +30,17 @@ def __init__(self, fnames_off, fnames_on): self.frames_off = np.asarray([Image(f).data for f in fnames_off]) self.frames_on = np.asarray([Image(f).data for f in fnames_on]) - off = np.median(self.frames_off, axis=0) - on = np.median(self.frames_on, axis=0) - flat = on - off - self.flat = flat / np.median(flat) + off = np.nanmedian(self.frames_off, axis=0) + on = np.nanmedian(self.frames_on, axis=0) + if (np.nanmean(on) - np.nanmean(off))/np.nanmean(on) < 0.01: + warnings.warn(f'brightness difference between domeflaton and domeflatoff is near zero \ + for flat including fname {fnames_on[0]}. setting to off. \ + if this is a thermal filter (lp, ms) then this is probably ok') + flat = off + else: + flat = on - off + + self.flat = flat / np.nanmedian(flat) def write(self, outfile): """ @@ -76,10 +83,10 @@ def make_badpx_map(self, outfile, tol=0.07, blocksize=6): warnings.filterwarnings('ignore', category=RuntimeWarning) med = np.median(flatblock) - # if not within tolerance, set to NaN - mapblock[np.where(flatblock / med > 1 + tol)] = 0 - mapblock[np.where(flatblock / med < 1 - tol)] = 0 - badpx_map[i : i + blocksize, j : j + blocksize] = mapblock + # if not within tolerance, set to NaN + mapblock[np.where(flatblock / med > 1 + tol)] = 0 + mapblock[np.where(flatblock / med < 1 - tol)] = 0 + badpx_map[i : i + blocksize, j : j + blocksize] = mapblock self.badpx_map = badpx_map # change some header info and write to .fits hdulist_out = self.dummy_fits.hdulist diff --git a/nirc2_reduce/multi_reduce.py b/nirc2_reduce/multi_reduce.py index de23366..64fb9f3 100755 --- a/nirc2_reduce/multi_reduce.py +++ b/nirc2_reduce/multi_reduce.py @@ -74,6 +74,7 @@ def process_flats( """ date_kw = self.header_kw_dict['date'] filter_kw = self.header_kw_dict['filter'] + subc_kw = self.header_kw_dict['subc'] wl_kw = self.header_kw_dict['wl'] standard_filts = np.array(list(standard_wleff.keys())) standard_wls = np.array([float(standard_wleff[key]) for key in standard_filts]) @@ -84,34 +85,37 @@ def process_flats( counter = 0 for i, tab in enumerate(tabs_by_filt): - flatoff, flaton = sort_rawfiles.get_flats(self.tab, self.instrument) - date = self.tab[date_kw].data[0] - - # if there are new flats in that filter, make master flat and badpx map - # and put them into flatdir - if (len(flatoff) > 0) and (len(flaton) > 0): - counter += 1 - - wl_eff = float(tab[wl_kw].data[0]) - standard_idx, _ = find_nearest(standard_wls, wl_eff) - nearest_stand_filt = standard_filts[standard_idx] - short_name = nearest_stand_filt - badpx_outpath = os.path.join( - flatdir, f"{date}_badpx_map_{short_name}.fits" - ) - flat_outpath = os.path.join( - flatdir, f"{date}_flat_master_{short_name}.fits" - ) - flatobj = flats.Flats(flatoff, flaton) - flatobj.write(flat_outpath) - flatobj.make_badpx_map(badpx_outpath, **badpx_kwargs) - print(f"wrote files {flat_outpath} and {badpx_outpath}") + + # sort by subarray + subcs, tabs_by_subc = sort_rawfiles.split_by_kw(tab, subc_kw) + for j, tab in enumerate(tabs_by_subc): + + flatoff, flaton = sort_rawfiles.get_flats(tab, self.instrument) + # if there are new flats in that filter, make master flat and badpx map + # and put them into flatdir + if (len(flatoff) > 0) and (len(flaton) > 0): + counter += 1 + date = tab[date_kw].data[0] + wl_eff = float(tab[wl_kw].data[0]) + standard_idx, _ = find_nearest(standard_wls, wl_eff) + nearest_stand_filt = standard_filts[standard_idx] + short_name = nearest_stand_filt + badpx_outpath = os.path.join( + flatdir, f"{date}_badpx_map_{subcs[j]}_{short_name}.fits" + ) + flat_outpath = os.path.join( + flatdir, f"{date}_flat_master_{subcs[j]}_{short_name}.fits" + ) + flatobj = flats.Flats(flatoff, flaton) + flatobj.write(flat_outpath) + flatobj.make_badpx_map(badpx_outpath, **badpx_kwargs) + print(f"wrote files {flat_outpath} and {badpx_outpath}") if counter == 0: warnings.warn("No flats processed! (none found by sort_rawfiles.get_flats)") return - def _find_flats(self, flatdir, date, filt): + def _find_flats(self, flatdir, date, filt, subc): """ find nearest-in-time flat and badpx map searches flatdir for nearest-in-time flats @@ -125,6 +129,7 @@ def _find_flats(self, flatdir, date, filt): directory to search for flats date : str, required. filt : str, required. + subc : str, required. Returns ------- @@ -135,15 +140,15 @@ def _find_flats(self, flatdir, date, filt): """ # find all the flats and bad pixel maps - all_flats = glob.glob(f"{flatdir}/*flat*{filt}*.fits") - all_badpx = glob.glob(f"{flatdir}/*badpx*{filt}*.fits") + all_flats = glob.glob(f"{flatdir}/*flat*{subc}*{filt}*.fits") + all_badpx = glob.glob(f"{flatdir}/*badpx*{subc}*{filt}*.fits") if len(all_flats) < 1: raise FileNotFoundError( - f"No flats in directory {flatdir} match wildcard *flat*{filt}*.fits" + f"No flats in directory {flatdir} match wildcard *flat*{subc}*{filt}*.fits" ) if len(all_badpx) < 1: raise FileNotFoundError( - f"No badpx maps in directory {flatdir} match wildcard *badpx*{filt}*.fits" + f"No badpx maps in directory {flatdir} match wildcard *badpx*{subc}*{filt}*.fits" ) flat_dates = [f.split("/")[-1].split("_")[0] for f in all_flats] badpx_dates = [f.split("/")[-1].split("_")[0] for f in all_badpx] @@ -182,7 +187,8 @@ def run( self, outdir, flatdir, - filts_want=None): + filts_want=None, + show=False): """ Parameters ---------- @@ -193,6 +199,8 @@ def run( filts_want : list-like or None, optional. default None. list of filter names to run if None, runs all filters found in rawdir + show : bool, optional. default False + show final quicklook images while running? Writes ------ @@ -200,10 +208,13 @@ def run( """ object_kw = self.header_kw_dict['object'] date_kw = self.header_kw_dict['date'] + time_kw = self.header_kw_dict['time'] wl_kw = self.header_kw_dict['wl'] + subc_kw = self.header_kw_dict['subc'] filter_kw = self.header_kw_dict['filter'] flatposkw = self.header_kw_dict['isdome']['kw'] flatposarg = self.header_kw_dict['isdome']['yesdome'] + trackingarg = self.header_kw_dict['isdome']['nodome'] if not os.path.exists(outdir): os.mkdir(outdir) @@ -212,11 +223,9 @@ def run( standard_filts = np.array(list(standard_wleff.keys())) standard_wls = np.array([float(standard_wleff[key]) for key in standard_filts]) - # ignore all files with dome in flat position + # select only files where scope is tracking, i.e., the science frames isinflatpos, flatpostabs = sort_rawfiles.split_by_kw(self.tab, flatposkw) - tab = flatpostabs[ - np.argwhere([not (s.startswith(flatposarg)) for s in isinflatpos])[0, 0] - ] + tab = flatpostabs[np.argwhere(isinflatpos == trackingarg)[0,0]] date = tab[date_kw].data[0] # split by object @@ -229,28 +238,32 @@ def run( filts, filt_tabs = sort_rawfiles.split_by_kw(targ_tab, filter_kw) # loop over all filters + failures = 0 for i, filt_name in enumerate(filts): filt_str = filt_name.replace(" ", "") filt_str = filt_str.replace("+", "") + filt_str = filt_str.replace("clear", "") targ_str = targ.split(" ")[0] if (filts_want is not None) and (filt_name not in filts_want): continue print(f"Starting filter {filt_name} ({i+1} of {len(filts)})") - failures = 0 filt_tab = filt_tabs[i] - fnames = np.sort(filt_tab["FILENAME"].data) + fnames_i = np.argsort(filt_tab["FILENAME"].data) + fnames = filt_tab["FILENAME"].data[fnames_i] + time_strs = filt_tab[time_kw].data[fnames_i] + time_strs = [s[:5].replace(":", "")+'UT' for s in time_strs] # check right number of frames for bxy3 if len(fnames) < 3: warnings.warn( - f"Fewer than 3 frames found for object {targ}, filter {filt_name}; skipping filter {filt_name}!" + f"Fewer than 3 frames found for object {targ}, filter {filt_name}; skipping filter {filt_name}!", stacklevel=2 ) failures += 1 continue if len(fnames) > 3: warnings.warn( - f"More than 3 frames found for object {targ}, filter {filt_name}; using last three!" + f"More than 3 frames found for object {targ}, filter {filt_name}; using last three!", stacklevel=2 ) fnames = fnames[-3:] @@ -258,11 +271,13 @@ def run( wl_eff = float(tab[wl_kw].data[0]) standard_idx, _ = find_nearest(standard_wls, wl_eff) flat_filt = standard_filts[standard_idx] + subc = filt_tab[subc_kw].data[0] try: - flat_fname, badpx_fname = self._find_flats(flatdir, date, flat_filt) + flat_fname, badpx_fname = self._find_flats(flatdir, date, flat_filt, subc) + print(f'Applying flat {flat_fname}, badpx {badpx_fname}') except (ValueError, FileNotFoundError): warnings.warn( - f"Could not find flat for filter {flat_filt} as requested by {date} {targ} {filt_name}; skipping filter {filt_name}!" + f"Could not find flat for filter {flat_filt} as requested by {date} {targ} {filt_name} {subc}; skipping filter {filt_name}!", stacklevel=2 ) failures += 1 continue @@ -296,16 +311,16 @@ def run( ) obs.trim() obs.stack() - obs.crop_final(50) - obs.plot_final( - show=False, - png_file=os.path.join(outdir, f"{date}_{targ_str}_{filt_str}.png"), - ) obs.write_final( os.path.join( outdir, f"{date}_{targ_str}_stacked_nophot_{filt_str}.fits" ) ) + obs.com_crop_final(100) + obs.plot_final( + show=show, + png_file=os.path.join(outdir, f"{targ_str}_{filt_str}_{time_strs[0]}.png"), + ) print(f"Object {targ_str} finished with {failures} failed filters") return diff --git a/nirc2_reduce/observation.py b/nirc2_reduce/observation.py index 3277748..132dddb 100755 --- a/nirc2_reduce/observation.py +++ b/nirc2_reduce/observation.py @@ -6,7 +6,7 @@ from .prettycolors import make_colormap, get_colormap from scipy.signal import medfilt from scipy.interpolate import RectBivariateSpline -from scipy.ndimage import rotate +from scipy.ndimage import rotate, center_of_mass import astroscrappy from image_registration.chi2_shifts import chi2_shift from image_registration.fft_tools.shift import shiftnd, shift2d @@ -83,6 +83,7 @@ def __init__(self, fnames, instrument): self.frames = np.asarray([Image(f).data for f in fnames]) self.subc = int(self.dummy_fits.header[self.header_kw_dict['subc']]) self.target = self.dummy_fits.header[self.header_kw_dict['object']] + self.filter = self.dummy_fits.header[self.header_kw_dict['filter']] if self.frames[0].shape[0] != self.subc or self.frames[0].shape[1] != self.subc: # subarrays smaller than 512x512 on Keck are not square. chop out center @@ -339,6 +340,48 @@ def crop_final(self, bw): """ szx, szy = self.final.shape[0], self.final.shape[1] self.final = self.final[bw : szx - bw, bw : szy - bw] + + def com_crop_final(self, wx, wy=None, cutoff=None): + """ + Attempts to crop final image around bright source in frame + using a center-of-mass approach + + Parameters + ---------- + wx : int, required. + distance from com to keep in +x and -x direction + i.e., half-width of output image + wy : int, optional. Default wx + distance from com to keep in +y and -y direction + cutoff : float, optional. Default None + the cutoff below which the data are considered noise + if None, code attempts to figure out rms noise based + on mean and stddev in the corners, and then + takes 5x that value as the cutoff + """ + if wy is None: + wy = wx + + # simple COM fails for large images because the few + # on-source data points do not do enough + # need to set the background to zero first + if cutoff is None: + a = 16 + a0 = int(self.final.shape[0]/a) + a1 = int(self.final.shape[1]/a) + corners = np.concatenate([self.final[:a0, :a1].flatten(), + self.final[-a0:, :a1].flatten(), + self.final[a0:, -a1:].flatten(), + self.final[-a0:, -a1:].flatten()]) + mn = np.mean(corners) + std = np.std(corners) + cutoff = np.abs(mn + std) + + data_to_calc = np.copy(self.final) + data_to_calc[data_to_calc < 10*cutoff] = 0.0 + com = center_of_mass(data_to_calc) + + self.final = self.final[int(com[0]-wy):int(com[0]+wy), int(com[1]-wx):int(com[1]+wx)] def plot_frames(self, png_file=None, figsz=3): """ @@ -381,9 +424,11 @@ def plot_final(self, show=True, png_file=None): try: cmap = get_colormap(self.target.split(" ")[0]) except: - print("No custom colormap defined for target, setting to default") + print(f"No custom colormap defined for target {self.target}, setting to default") cmap = cm.viridis ax.imshow(self.final, cmap=cmap, origin="lower") + ax.set_xticks([]) + ax.set_yticks([]) plt.tight_layout() if png_file is not None: fig.savefig(png_file, dpi=300) @@ -565,7 +610,7 @@ def apply_sky(self, fname): frames_skysub.append(frame - sky_norm) self.frames = np.asarray(frames_skysub) else: - print("Sky subtraction not applied (counts too low for good stats)") + print(f"Sky subtraction not applied to {self.target} {self.filter} {self.subc} (counts too low for good stats)") def trim(self): """ diff --git a/nirc2_reduce/sort_rawfiles.py b/nirc2_reduce/sort_rawfiles.py index 1c88969..a6337fb 100755 --- a/nirc2_reduce/sort_rawfiles.py +++ b/nirc2_reduce/sort_rawfiles.py @@ -82,8 +82,8 @@ def split_by_kw(tab, kw): def get_flats( tab, instrument, - ignore_objects=["FLAT", "BADPX", "FLAT_MASTER", "DOME_FLAT_MASTER", "BADPX_MAP"], -): + ignore_objects=['',] + ): """ scrub table from dfits_fitsort to find domeflaton, domeflatoff filenames default params are for NIRC2 @@ -97,8 +97,8 @@ def get_flats( name of instrument used, e.g. nirc2. will look for file data/{instrument}.yaml in order to scrub header keywords - ignore_objects : list, optional. - special values in header[object] to ignore + ignore_objects : list, optional, default empty list. + special values in header[object] to ignore. Returns ------- @@ -119,12 +119,30 @@ def get_flats( # find dome position targnames, tabs = split_by_kw(tab, header_kw_dict['isdome']['kw']) - dometab = tabs[np.argwhere([s.startswith(header_kw_dict['isdome']['yesdome'].strip()) for s in targnames])[0, 0]] + domebool = [s.startswith(header_kw_dict['isdome']['yesdome'].strip()) for s in targnames] + if not np.any(np.array(domebool)): + return [], [] + dometab = tabs[np.argwhere(domebool)[0, 0]] + + # remove darks using shutter on/off kw + targnames, darktabs = split_by_kw(dometab, header_kw_dict['shutter']['kw']) + targnames = np.array([s.strip() for s in targnames]) + shutbool = (targnames == header_kw_dict['shutter']['shutopen']) + if not np.any(np.array(shutbool)): + return [], [] + flattab = darktabs[np.argwhere(shutbool)[0, 0]] # find ons and offs - onoff, dometabs = split_by_kw(dometab, header_kw_dict['lamps']['kw']) - ons = dometabs[np.argwhere(onoff == header_kw_dict['lamps']['lampon'])[0, 0]] - offs = dometabs[np.argwhere(onoff == header_kw_dict['lamps']['lampoff'])[0, 0]] + onoff, flattabs = split_by_kw(flattab, header_kw_dict['lamps']['kw']) + onbool = [s.startswith(header_kw_dict['lamps']['lampon'].strip()) for s in onoff] + if not np.any(np.array(onbool)): + print('Found some frames taken in dome flat position, but no frames with lamps on!') + print('Might be AO calibration; simply skipping make dome flats step') + print('But if you were expecting dome flats, something went wrong!') + return [], [] + + ons = flattabs[np.argwhere(onoff == header_kw_dict['lamps']['lampon'])[0, 0]] + offs = flattabs[np.argwhere(onoff == header_kw_dict['lamps']['lampoff'])[0, 0]] flatoff = offs["FILENAME"].data flaton = ons["FILENAME"].data diff --git a/nirc2_reduce/tests/test_multi_reduce.py b/nirc2_reduce/tests/test_multi_reduce.py index d8a1f31..508ff86 100644 --- a/nirc2_reduce/tests/test_multi_reduce.py +++ b/nirc2_reduce/tests/test_multi_reduce.py @@ -52,13 +52,13 @@ def test_MultiBxy3(datadir, rawdir, reddir, flatdir): ].data assert np.allclose(stack_test, stack_expected, rtol=1e-3) - flat_test = fits.open(os.path.join(flatdir, "2017-07-25_flat_master_h.fits"))[ + flat_test = fits.open(os.path.join(flatdir, "2017-07-25_flat_master_1024_h.fits"))[ 0 ].data flat_expected = fits.open(os.path.join(datadir, "flat_expected.fits"))[0].data assert np.allclose(flat_test, flat_expected, rtol=1e-3) - badpx_test = fits.open(os.path.join(flatdir, "2017-07-25_badpx_map_h.fits"))[0].data + badpx_test = fits.open(os.path.join(flatdir, "2017-07-25_badpx_map_1024_h.fits"))[0].data badpx_expected = fits.open(os.path.join(datadir, "badpx_map_expected.fits"))[0].data assert np.allclose(badpx_test, badpx_expected, rtol=1e-3) diff --git a/nirc2_reduce/tests/test_sort_rawfiles.py b/nirc2_reduce/tests/test_sort_rawfiles.py index 682f2e8..e63ee1c 100644 --- a/nirc2_reduce/tests/test_sort_rawfiles.py +++ b/nirc2_reduce/tests/test_sort_rawfiles.py @@ -25,7 +25,7 @@ def test_dfits_fitsort(rawdir): warnings.simplefilter("ignore") input_wildcard = os.path.join(rawdir, "o*.fits") tab = sort_rawfiles.dfits_fitsort(input_wildcard, - ["OBJECT", "DATE-OBS", "FILTER", "AXESTAT", "FLSPECTR"]) + ["OBJECT", "DATE-OBS", "FILTER", "AXESTAT", "FLSPECTR", "SHRNAME"]) flatoff, flaton = sort_rawfiles.get_flats(tab, 'nirc2') offs_expected = np.array([os.path.join(rawdir, f"off{i}.fits") for i in range(5)])