diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 80c8730..ffb9fc2 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -10,6 +10,16 @@ on: branches: [ main, master ] jobs: + + pre-commit: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@08c6903cd8c0fde910a37f88322edcfb5dd907a8 # v5.0.0 + - uses: actions/setup-python@a26af69be951a213d495a4c3e4e4022e16d87065 # v5.6.0 + with: + python-version: '3.12' + - uses: pre-commit/action@2c7b3805fd2a0fd8c1884dcaebf91fc102a13ecd # v3.0.1 + build: runs-on: ubuntu-latest diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..63eb5be --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,22 @@ +repos: + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v5.0.0 + hooks: + - id: check-added-large-files + - id: check-ast + - id: check-case-conflict + - id: check-yaml + args: ["--unsafe"] + - id: check-toml + - id: check-merge-conflict + - id: check-symlinks + - id: debug-statements + - id: detect-private-key + # - id: end-of-file-fixer + # - id: trailing-whitespace + - repo: https://github.com/astral-sh/ruff-pre-commit + rev: 'v0.12.3' + hooks: + - id: ruff + args: ["--fix"] + - id: ruff-format \ No newline at end of file diff --git a/.ruff.toml b/.ruff.toml new file mode 100644 index 0000000..1314cef --- /dev/null +++ b/.ruff.toml @@ -0,0 +1,66 @@ +extend = "pyproject.toml" + +exclude = [ + "docs", + ".tox", + ".eggs", + "build", + "*.ipynb", + "**/tests/**/__init__.py", +] +line-length = 100 + +target-version = "py312" + +[format] +quote-style = "double" +indent-style = "space" +docstring-code-format = true + +[lint] +select = [ + "F", # Pyflakes (part of default flake8) + "E", # pycodestyle (part of default flake8) + "W", # pycodestyle (part of default flake8) + "D", # docstrings, see also numpydoc pre-commit action + "N", # pep8-naming (naming conventions) + "A", # flake8-builtins (prevent shadowing of builtins) + "ARG", # flake8-unused-arguments (prevent unused arguments) + "B", # flake8-bugbear (miscellaneous best practices to avoid bugs) + "C4", # flake8-comprehensions (best practices for comprehensions) + "I", # isort + "ICN", # flake8-import-conventions (enforce import conventions) + "INP", # flake8-no-pep420 (prevent use of PEP420, i.e. implicit name spaces) + "ISC", # flake8-implicit-str-concat (conventions for concatenating long strings) + "LOG", # flake8-logging + "NPY", # numpy-specific rules + "PGH", # pygrep-hooks (ensure appropriate usage of noqa and type-ignore) + # "PTH", # flake8-use-pathlib (enforce using Pathlib instead of os) + "S", # flake8-bandit (security checks) + "SLF", # flake8-self (prevent using private class members outside class) + "SLOT", # flake8-slots (require __slots__ for immutable classes) + # "T20", # flake8-print (prevent print statements in code) + "TRY", # tryceratops (best practices for try/except blocks) + "UP", # pyupgrade (simplified syntax allowed by newer Python versions) + "YTT", # flake8-2020 (prevent some specific gotchas from sys.version) + "TID252" # Checks for relative imports +] +ignore = ["D100", "D104", "TRY003"] + +[lint.per-file-ignores] +"**/tests/test_*.py" = [ + "D", + "S", +] + +[lint.pydocstyle] +convention = "numpy" + +[lint.flake8-annotations] +ignore-fully-untyped = true # Turn off annotation checking for fully untyped code + +[lint.flake8-tidy-imports] +ban-relative-imports = "all" # Disallow all relative imports. + +[lint.isort] +known-first-party = ["nirc2_reduce"] \ No newline at end of file diff --git a/nirc2_reduce/_version.py b/nirc2_reduce/_version.py index 4f813b9..f826ffa 100644 --- a/nirc2_reduce/_version.py +++ b/nirc2_reduce/_version.py @@ -12,11 +12,8 @@ TYPE_CHECKING = False if TYPE_CHECKING: - from typing import Tuple - from typing import Union - - VERSION_TUPLE = Tuple[Union[int, str], ...] - COMMIT_ID = Union[str, None] + VERSION_TUPLE = tuple[int | str, ...] + COMMIT_ID = str | None else: VERSION_TUPLE = object COMMIT_ID = object @@ -28,7 +25,7 @@ commit_id: COMMIT_ID __commit_id__: COMMIT_ID -__version__ = version = '0.1.dev53+g310c26299.d20250816' -__version_tuple__ = version_tuple = (0, 1, 'dev53', 'g310c26299.d20250816') +__version__ = version = "0.1.dev55+g3741c605e.d20250817" +__version_tuple__ = version_tuple = (0, 1, "dev55", "g3741c605e.d20250817") -__commit_id__ = commit_id = 'g310c26299' +__commit_id__ = commit_id = "g3741c605e" diff --git a/nirc2_reduce/flats.py b/nirc2_reduce/flats.py index 8c43344..dcac8cb 100755 --- a/nirc2_reduce/flats.py +++ b/nirc2_reduce/flats.py @@ -1,24 +1,25 @@ -import numpy as np -import matplotlib.pyplot as plt -from astropy.io import fits -from .image import Image import warnings +import matplotlib.pyplot as plt +import numpy as np + +from nirc2_reduce.image import Image + + class Flats: - """ - class for dome flats - """ + """Class for dome flats.""" def __init__(self, fnames_off, fnames_on): """ - subtract lights off from lights on, - divide by median value to center it around 1 - + Subtract lights off from lights on. + + Divide by median value to center it around 1 + Parameters ---------- fnames_off : list, required. fits files for dome flat OFF fnames_on : list, required. fits files for dome flat ON - + Attributes ---------- dummy_fits : Image object used to hijack header info @@ -32,16 +33,23 @@ def __init__(self, fnames_off, fnames_on): 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', stacklevel=2) + if (np.nanmean(on) - np.nanmean(off)) / np.nanmean(on) < 0.01: + warnings.warn( + "brightness difference between domeflaton and domeflatoff is near zero " + f"for flat including fname {fnames_on[0]}. setting to off. " + "if this is a thermal filter (lp, ms) then this is probably ok", + stacklevel=2, + ) flat = off else: flat = on - off - + self.flat = flat / np.nanmedian(flat) def write(self, outfile): """ + Write the flat to fits. + Parameters ---------- outfile : str, required. output fits filename @@ -53,23 +61,25 @@ def write(self, outfile): hdulist_out[0].writeto(outfile, overwrite=True) def plot(self): + """Make a simple plot of the flat.""" plt.imshow(self.flat, origin="lower") plt.show() def make_badpx_map(self, outfile, tol=0.07, blocksize=6): """ - Find pixels whose values are very far from the average of their neighbors - Bad pixel is defined as + Find pixels whose values are very far from the average of their neighbors. + + Bad pixel is defined as abs(pixel value / median of nearby pixels - 1) > tol - + Parameters ---------- - outfile : str, required. + outfile : str, required. fits filename to write to tol : float, optional. default 0.07 - fractional tolerance. + fractional tolerance. blocksize : int, optional. default 6 - number of pixels over which to average + number of pixels over which to average in each direction """ badpx_map = np.ones(self.flat.shape) @@ -78,7 +88,7 @@ def make_badpx_map(self, outfile, tol=0.07, blocksize=6): flatblock = self.flat[i : i + blocksize, j : j + blocksize] mapblock = badpx_map[i : i + blocksize, j : j + blocksize] with warnings.catch_warnings(): - warnings.filterwarnings('ignore', category=RuntimeWarning) + warnings.filterwarnings("ignore", category=RuntimeWarning) med = np.median(flatblock) # if not within tolerance, set to NaN diff --git a/nirc2_reduce/image.py b/nirc2_reduce/image.py index 578a472..0fd4536 100755 --- a/nirc2_reduce/image.py +++ b/nirc2_reduce/image.py @@ -1,14 +1,13 @@ -import numpy as np +import warnings + import matplotlib.pyplot as plt from astropy.io import fits -import warnings class Image: """ - Thin wrapper to astropy.io.fits providing convenience functions - and suppression of warnings in fits.open - + Thin wrapper to astropy.io.fits. + To do ----- make agnostic to specific .fits header keywords @@ -17,10 +16,12 @@ class Image: def __init__(self, fname): """ + Provide convenience function for opening fits files. + Parameters ---------- fname : str, required. input .fits filename - + Attributes ---------- hdulist: .fits format hdulist structure @@ -31,7 +32,8 @@ def __init__(self, fname): """ warnings.filterwarnings( "ignore", - "The following header keyword is invalid or follows an unrecognized non-standard convention", + "The following header keyword is invalid or follows an unrecognized " + "non-standard convention", ) warnings.filterwarnings( "ignore", @@ -39,7 +41,8 @@ def __init__(self, fname): ) warnings.filterwarnings( "ignore", - "Header block contains null bytes instead of spaces for padding, and is not FITS-compliant", + "Header block contains null bytes instead of spaces for padding, " + "and is not FITS-compliant", ) self.hdulist = fits.open(fname, ignore_missing_end=True) self.header = self.hdulist[0].header @@ -51,7 +54,7 @@ def __init__(self, fname): try: targ = self.header["OBJECT"] self.target = targ.split()[0].strip(", \n").capitalize() - except: + except KeyError: self.target = "Unknown" try: fwi = self.header["FWINAME"].strip(", \n") @@ -60,14 +63,25 @@ def __init__(self, fname): self.filt = fwi elif fwi == "clear" or fwi == "PK50_1.5": self.filt = fwo - except: + else: + self.filt = "Unknown" + except KeyError: self.filt = "Unknown" def plot(self): + """Make a simple plot of the data.""" plt.imshow(self.data, origin="lower left") plt.show() def write(self, fname): + """ + Write the image to fits. + + Parameters + ---------- + fname : str + Filename to write to. + """ self.hdulist[0].header = self.header self.hdulist[0].data = self.data self.hdulist.writeto(fname, overwrite=True) diff --git a/nirc2_reduce/multi_reduce.py b/nirc2_reduce/multi_reduce.py index 9b59924..37fc942 100755 --- a/nirc2_reduce/multi_reduce.py +++ b/nirc2_reduce/multi_reduce.py @@ -1,42 +1,49 @@ -from . import sort_rawfiles, observation, flats, image -import matplotlib.pyplot as plt -import numpy as np import glob +import importlib import os import warnings -from astropy.time import Time -import importlib + +import numpy as np import yaml +from astropy.time import Time + import nirc2_reduce.data.header_kw_dicts as inst_info +from nirc2_reduce import flats, observation, sort_rawfiles -''' +""" To do: instead of nearest_stand_filt, use header's effective wavelength! -''' -standard_wleff = {'j': 1.248, 'h': 1.633, 'kp': 2.124, 'lp':3.776, 'ms':4.674,} +""" +standard_wleff = { + "j": 1.248, + "h": 1.633, + "kp": 2.124, + "lp": 3.776, + "ms": 4.674, +} + -def find_nearest(array, value): +def _find_nearest(array, value): array = np.asarray(array) idx = (np.abs(array - value)).argmin() return idx, array[idx] + class MultiReduce: - """ - - - superclass to multiBxy3 and multiNod - """ + """Superclass to multiBxy3 and multiNod.""" - def __init__(self,rawdir,instrument): + def __init__(self, rawdir, instrument): """ + Create a MultiReduce object. + Parameters ---------- rawdir : str, required. directory containing the raw fits files instrument : str, required. - name of instrument used, e.g. nirc2. + name of instrument used, e.g. nirc2. will look for file data/{instrument}.yaml in order to scrub header keywords - + Attributes ---------- rawdir : str @@ -50,46 +57,42 @@ def __init__(self,rawdir,instrument): yaml_bytes = file.read() self.header_kw_dict = yaml.safe_load(yaml_bytes) self.tab = sort_rawfiles.dfits_fitsort( - os.path.join(self.rawdir, "*.fits"), fits_kws=self.header_kw_dict['dfits_kws']) - + os.path.join(self.rawdir, "*.fits"), fits_kws=self.header_kw_dict["dfits_kws"] + ) - def process_flats( - self, - flatdir, - **badpx_kwargs): + def process_flats(self, flatdir, **badpx_kwargs): """ - check for new dome flats in rawdir, + Process the flatfield files. + + Check for new dome flats in rawdir, make new flats for each filter present and put into flatdir - + Parameters ---------- flatdir : str, required. directory to put into and/or retrieve flats if retrieving flats without making them, ensure - filenames match defaults; + filenames match defaults; {date}_badpx_map_{filt}.fits for bad pixel maps in given filter {date}_flat_master_{filt}.fits for flats in given filter **badpx_kwargs : dict, optional. keyword arguments to flats.Flats.make_badpx_map() """ - 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'] + 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]) - - + # sort by filter - filts, tabs_by_filt = sort_rawfiles.split_by_kw(self.tab, filter_kw) + _filts, tabs_by_filt = sort_rawfiles.split_by_kw(self.tab, filter_kw) counter = 0 - for i, tab in enumerate(tabs_by_filt): - + for tab in tabs_by_filt: # 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 @@ -97,7 +100,7 @@ def process_flats( counter += 1 date = tab[date_kw].data[0] wl_eff = float(tab[wl_kw].data[0]) - standard_idx, _ = find_nearest(standard_wls, wl_eff) + 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( @@ -111,18 +114,21 @@ def process_flats( 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)") + warnings.warn( + "No flats processed! (none found by sort_rawfiles.get_flats)", stacklevel=2 + ) return def _find_flats(self, flatdir, date, filt, subc): """ - find nearest-in-time flat and badpx map - searches flatdir for nearest-in-time flats + Find nearest-in-time flat and badpx map. + + Searches flatdir for nearest-in-time flats and badpx maps. assumes filenames in flats/ match the wildcards - DATE-OBS*flat*{filt}*.fits and + DATE-OBS*flat*{filt}*.fits and DATE-OBS*badpx*{filt}*.fits - + Parameters ---------- flatdir : str, required. @@ -130,7 +136,7 @@ def _find_flats(self, flatdir, date, filt, subc): date : str, required. filt : str, required. subc : str, required. - + Returns ------- str @@ -138,7 +144,6 @@ def _find_flats(self, flatdir, date, filt, subc): str filename of nearest-in-time badpx map """ - # find all the flats and bad pixel maps all_flats = glob.glob(f"{flatdir}/*flat*{subc}*{filt}*.fits") all_badpx = glob.glob(f"{flatdir}/*badpx*{subc}*{filt}*.fits") @@ -162,7 +167,8 @@ def _find_flats(self, flatdir, date, filt, subc): # ensure that flat has an associated bad pixel map if nearest not in badpx_dates: raise ValueError( - f"Nearest-in-time flat ({nearest}) has no associated bad pixel map (unmatched wildcard {nearest}*badpx*{filt}.fits)" + f"Nearest-in-time flat ({nearest}) has no associated bad pixel map " + f"(unmatched wildcard {nearest}*badpx*{filt}.fits)" ) return all_flats[nearest_idx], all_badpx[nearest_idx] @@ -170,8 +176,8 @@ def _find_flats(self, flatdir, date, filt, subc): class MultiBxy3(MultiReduce): """ - run multiple filters and objects in a single date of observing - + Run multiple filters and objects in a single date of observing. + Example ------- obs = MultiBxy3(rawdir) @@ -180,16 +186,12 @@ class MultiBxy3(MultiReduce): """ def __init__(self, rawdir, instrument): - super().__init__(rawdir, instrument) - def run( - self, - outdir, - flatdir, - filts_want=None, - show=False): + def run(self, outdir, flatdir, filts_want=None, show=False): """ + Run the multiple reductions. + Parameters ---------- outdir : str, required. @@ -201,31 +203,31 @@ def run( if None, runs all filters found in rawdir show : bool, optional. default False show final quicklook images while running? - + Writes ------ - + """ - 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'] - + 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) - - # read in standard filters + + # read in standard filters standard_filts = np.array(list(standard_wleff.keys())) standard_wls = np.array([float(standard_wleff[key]) for key in standard_filts]) # 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(isinflatpos == trackingarg)[0,0]] + tab = flatpostabs[np.argwhere(isinflatpos == trackingarg)[0, 0]] date = tab[date_kw].data[0] # split by object @@ -233,7 +235,7 @@ def run( # loop over all targets for j, targ in enumerate(targets): - print(f"Starting object {targ} ({j+1} of {len(targets)})") + print(f"Starting object {targ} ({j + 1} of {len(targets)})") targ_tab = target_tabs[j] filts, filt_tabs = sort_rawfiles.split_by_kw(targ_tab, filter_kw) @@ -247,50 +249,52 @@ def run( 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)})") + print(f"Starting filter {filt_name} ({i + 1} of {len(filts)})") filt_tab = filt_tabs[i] 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.split('T')[-1] for s in time_strs] - time_strs = [s[:5].replace(":", "")+'UT' for s in time_strs] + time_strs = [s.split("T")[-1] for s in time_strs] + 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}!", stacklevel=2 + 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!", stacklevel=2 + f"More than 3 frames found for object {targ}, " + f"filter {filt_name}; using last three!", + stacklevel=2, ) fnames = fnames[-3:] # find corresponding wideband flatfield and badpx map wl_eff = float(filt_tab[wl_kw].data[0]) - standard_idx, _ = find_nearest(standard_wls, wl_eff) + 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, subc) - print(f'Applying flat {flat_fname}, badpx {badpx_fname}') + 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} {subc}; skipping filter {filt_name}!", stacklevel=2 + f"Could not find flat for filter {flat_filt} as requested by " + f"{date} {targ} {filt_name} {subc}; skipping filter {filt_name}!", + stacklevel=2, ) failures += 1 continue # run all steps of bxy3 obs = observation.Bxy3(fnames, self.instrument) - obs.make_sky( - os.path.join(outdir, f"{date}_{targ_str}_sky_{filt_str}.fits") - ) - obs.apply_sky( - os.path.join(outdir, f"{date}_{targ_str}_sky_{filt_str}.fits") - ) + obs.make_sky(os.path.join(outdir, f"{date}_{targ_str}_sky_{filt_str}.fits")) + obs.apply_sky(os.path.join(outdir, f"{date}_{targ_str}_sky_{filt_str}.fits")) obs.apply_flat(flat_fname) obs.apply_badpx_map(badpx_fname) obs.remove_cosmic_rays() @@ -300,22 +304,14 @@ def run( obs.trim() obs.write_frames( [ - os.path.join( - outdir, f"{date}_{targ_str}_0_nophot_{filt_str}.fits" - ), - os.path.join( - outdir, f"{date}_{targ_str}_1_nophot_{filt_str}.fits" - ), - os.path.join( - outdir, f"{date}_{targ_str}_2_nophot_{filt_str}.fits" - ), + os.path.join(outdir, f"{date}_{targ_str}_0_nophot_{filt_str}.fits"), + os.path.join(outdir, f"{date}_{targ_str}_1_nophot_{filt_str}.fits"), + os.path.join(outdir, f"{date}_{targ_str}_2_nophot_{filt_str}.fits"), ], ) obs.stack() obs.write_final( - os.path.join( - outdir, f"{date}_{targ_str}_stacked_nophot_{filt_str}.fits" - ) + os.path.join(outdir, f"{date}_{targ_str}_stacked_nophot_{filt_str}.fits") ) obs.com_crop_final(100) obs.plot_final( @@ -327,13 +323,11 @@ def run( return -class MultiNod(MultiReduce): - def __init__(self, rawdir): - - raise NotImplementedError() +# class MultiNod(MultiReduce): +# def __init__(self, rawdir): +# raise NotImplementedError() - super().__init__(rawdir) +# super().__init__(rawdir) - def run(self): - - return +# def run(self): +# return diff --git a/nirc2_reduce/observation.py b/nirc2_reduce/observation.py index 634e772..d0e0fa5 100755 --- a/nirc2_reduce/observation.py +++ b/nirc2_reduce/observation.py @@ -1,36 +1,39 @@ -import numpy as np +import importlib.resources +import warnings +from datetime import datetime + +import astroscrappy import matplotlib.pyplot as plt -from matplotlib import cm +import numpy as np +import yaml from astropy.io import fits -from .image import Image -from .prettycolors import get_colormap -from scipy.signal import medfilt -from scipy.interpolate import RectBivariateSpline -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 shift2d -import importlib.resources -import warnings -import yaml -import nirc2_reduce.data.header_kw_dicts as inst_info -from datetime import datetime +from matplotlib import cm +from scipy.interpolate import RectBivariateSpline +from scipy.ndimage import center_of_mass, rotate +from scipy.signal import medfilt +import nirc2_reduce.data.header_kw_dicts as inst_info +from nirc2_reduce.image import Image +from nirc2_reduce.prettycolors import get_colormap -''' +""" To do ----- weight stacks by integration time (?) -''' +""" def crop_center(frame, subc): """ + Crop an array around its center. + Parameters ---------- frame : np.array, required. subc : int, required. - + Returns ------- np.array @@ -45,23 +48,26 @@ def crop_center(frame, subc): class Observation: """ - Generic class containing dithers and nods + Generic class containing dithers and nods. + To define a custom dither pattern, inherit from this class See Bxy3() and Nod() objects for usage """ def __init__(self, fnames, instrument): """ + Make a generic infrared imaging observation. + Parameters ---------- - fnames : list or str, required. - filenames of data frames. + fnames : list or str, required. + filenames of data frames. if str, assume passing single data frame instrument : str, required. - name of instrument used, e.g. nirc2. + name of instrument used, e.g. nirc2. will look for file data/{instrument}.yaml in order to scrub header keywords - + Attributes ---------- dummy_fits : nirc2_reduce.image.Image of the zeroth frame @@ -71,7 +77,9 @@ def __init__(self, fnames, instrument): target : str. the object you observed, scrubbed from the header. """ if isinstance(fnames, str): - fnames = [fnames,] + fnames = [ + fnames, + ] self.dummy_fits = Image(fnames[0]) # used to hijack header info date = self.dummy_fits.date self.instrument = instrument.lower() @@ -80,7 +88,7 @@ def __init__(self, fnames, instrument): "Instrument is nirc2, but DATE-OBS could not be read from header." "Either choose instrument nirc2_pre_oct23 or nirc2_post_oct23, " "or ensure the date can be read from the header." - ) + ) if instrument == "nirc2": date = datetime.strptime(date, "%Y-%m-%d") controller_update_date = datetime.strptime("2023-10-15", "%Y-%m-%d") @@ -93,9 +101,9 @@ def __init__(self, fnames, instrument): self.header_kw_dict = yaml.safe_load(yaml_bytes) 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']] + 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 @@ -106,11 +114,11 @@ def __init__(self, fnames, instrument): def apply_flat(self, fname): """ - applies flatfield correction - + Apply the flatfield correction. + Parameters ---------- - fname : str, required. + fname : str, required. filename of flat """ flat = Image(fname).data @@ -123,17 +131,16 @@ def apply_flat(self, fname): with np.errstate(divide="ignore", invalid="ignore"): frame_flat = frame / flat # solve NaNs, but first ensure we didn't get a huge fraction of NaNs - if np.sum(np.isnan(frame_flat))/(frame_flat.size) < 0.1: + if np.sum(np.isnan(frame_flat)) / (frame_flat.size) < 0.1: frame_flat[np.isnan(frame_flat)] = 0.0 frames_flat.append(frame_flat) - + self.frames = np.asarray(frames_flat) def apply_badpx_map(self, fname, kernel_size=7): """ - Replaces all bad pixels in input map with the median of the pixels - around it. - + Replace all bad pixels in input map with the median of the pixels around it. + Parameters ---------- fname : str, required. @@ -146,7 +153,8 @@ def apply_badpx_map(self, fname, kernel_size=7): if badpx_map.shape[0] != self.subc: badpx_map = crop_center(badpx_map, self.subc) bad_indices = np.where(badpx_map == 0) - ## test case - kernel of 7 shown to remove all bad pixels including weird stripe in bottom left of detector + ## test case - kernel of 7 shown to remove all bad pixels + ## including weird stripe in bottom left of detector # smoothed_badpx_map = medfilt(badpx_map,kernel_size = 7) # fig, (ax0,ax1) = plt.subplots(1,2, figsize=(10,6)) # ax0.imshow(badpx_map, origin = 'lower') @@ -161,34 +169,36 @@ def apply_badpx_map(self, fname, kernel_size=7): def dewarp(self): """ - Apply nirc2 distortion correction using updated maps + Apply nirc2 distortion correction. + + Uses updated maps from Ghez, Lu galactic center group (Service et al. 2016) doi:10.1088/1538-3873/128/967/095004 - Spline interpolation used to interpolate original image. + Spline interpolation used to interpolate original image. The spline is then evaluated at new pixel locations computed from the map. - Custom distortion solutions should be placed in the + Custom distortion solutions should be placed in the nirc2_reduce/nirc2_reduce/data/ folder - + Parameters ---------- warpx_file : str, optional. distortion map to apply in x warpy_file : str, optional. distortion map to apply in y - + Notes ----- - From Service+16 "These are lookup tables generated by evaluating the fits - from the previous section at the center of every pixel on the - NIRC2 detector and are the values that should be added to + From Service+16 "These are lookup tables generated by evaluating the fits + from the previous section at the center of every pixel on the + NIRC2 detector and are the values that should be added to raw NIRC2 positions to shift them to a distortion-free frame" - So the position offsets in the fits files should be added. - That this works right was checked by looking at difference maps + So the position offsets in the fits files should be added. + That this works right was checked by looking at difference maps in a subdirectory of tests/ and comparing with the expected vectors in Service+16 """ - warpx_file = self.header_kw_dict['warpx_file'] - warpy_file = self.header_kw_dict['warpy_file'] + warpx_file = self.header_kw_dict["warpx_file"] + warpy_file = self.header_kw_dict["warpy_file"] distortion_source_x = importlib.resources.open_binary( "nirc2_reduce.data.distortion", f"{warpx_file}" ) @@ -202,7 +212,7 @@ def dewarp(self): if ( self.subc != warpx.shape[0] ): # hardcode because dewarp arrays will always be full size of detector - ctr = warpx.shape[0]/2 + ctr = warpx.shape[0] / 2 ll = int(ctr - self.subc / 2) ul = int(ctr + self.subc / 2) warpx = warpx[ll:ul, ll:ul] @@ -218,7 +228,7 @@ def dewarp(self): flatx = mapx.flatten() flaty = mapy.flatten() frames_dewarp = [] - + for frame in self.frames: frame_spline = RectBivariateSpline(xx, yy, frame) frame_dw = frame_spline.ev(flatx, flaty).reshape(frame.shape).T @@ -229,12 +239,13 @@ def dewarp(self): def rotate(self, custom_angle=0.0, beta=0.252): """ - Rotate frame to North up plus an additional custom angle + Rotate frame to North up plus an additional custom angle. + if header contains 'ROTPOSN' and 'INSTANGL' keywords (i.e., nirc2-like), then rotation is custom_angle + (ROTPOSN-INSTANGL) - beta, counterclockwise. otherwise, rotation is simply custom_angle - beta, counterclockwise. see https://github.com/jluastro/nirc2_distortion/wiki - + Parameters ---------- custom_angle : float, optional. default 0.0. units degrees. @@ -244,7 +255,7 @@ def rotate(self, custom_angle=0.0, beta=0.252): rotation required to get North up according to astrometry solution clockwise is positive to agree with definition of beta given by the nirc2 distortion wiki - + Notes ----- This should be generalized to include scopes other than nirc2 more easily @@ -255,7 +266,9 @@ def rotate(self, custom_angle=0.0, beta=0.252): total_rotation_ccw = custom_angle + hdr["ROTPOSN"] - hdr["INSTANGL"] - beta except KeyError: warnings.warn( - "header keywords ROTPOSN and/or INSTANGL not found; setting rotation to custom_angle - beta" + "header keywords ROTPOSN and/or INSTANGL not found; " + "setting rotation to custom_angle - beta", + stacklevel=2, ) self.frames = np.array( [rotate(frame, total_rotation_ccw, reshape=False) for frame in self.frames] @@ -263,8 +276,8 @@ def rotate(self, custom_angle=0.0, beta=0.252): def remove_cosmic_rays(self, **kwargs): """ - Detects cosmic rays using the astroscrappy package - + Detect cosmic rays using the astroscrappy package. + Parameters ---------- **kwargs : dict, optional @@ -274,7 +287,7 @@ def remove_cosmic_rays(self, **kwargs): for frame in self.frames: if np.all(np.isnan(frame)): # catches seg fault when passing all nans to detect_cosmics - raise ValueError('Tried to pass all NaNs to astroscrappy.detect_cosmics') + raise ValueError("Tried to pass all NaNs to astroscrappy.detect_cosmics") crmask, cleanarr = astroscrappy.detect_cosmics(frame, cleantype="medmask", **kwargs) # fig, (ax0,ax1,ax2) = plt.subplots(1,3, figsize=(15,6)) # ax0.imshow(frame,origin = 'lower') @@ -286,9 +299,8 @@ def remove_cosmic_rays(self, **kwargs): def per_second(self): """ - Changes units to counts/second by dividing by - ITIME x COADDS - + Change units to counts/second by dividing by ITIME x COADDS. + Parameters ---------- itime_kw : str, optional. @@ -297,8 +309,8 @@ def per_second(self): keyword to scrub fits header for coadds """ header = self.dummy_fits.header - tint = header[self.header_kw_dict['itime']] - coadd = header[self.header_kw_dict['coadds']] + tint = header[self.header_kw_dict["itime"]] + coadd = header[self.header_kw_dict["coadds"]] persec_frames = [] for frame in self.frames: persec_frames.append(frame / (tint * coadd)) @@ -306,8 +318,8 @@ def per_second(self): def apply_photometry_frames(self, flux_per): """ - Simply multiplies each frame by flux density / (cts/s) - + Multiply each frame by flux density / (cts/s). + Parameters ---------- flux_per : float, required. units [flux density / (cts/s)] @@ -317,25 +329,27 @@ def apply_photometry_frames(self, flux_per): def apply_photometry_final(self, flux_per): """ - Simply multiplies each frame by flux density / (cts/s) - + Multiply final combined product by flux density / (cts/s). + Parameters ---------- flux_per : float, required. units [flux density / (cts/s)] - conversion factor + conversion factor """ self.final = self.final * flux_per def stack(self): """ Cross-correlate the images applying sub-pixel shift. + Shift found using DFT upsampling method from image_registration. Stack them on top of each other to increase SNR. """ shifted_data = [self.frames[0]] for frame in self.frames[1:]: - [dx, dy, dxerr, dyerr] = chi2_shift(self.frames[0], frame) - # error is nonzero only if you include per-pixel error of each image as an input. Should eventually do that, but no need for now. + # error is nonzero only if you include per-pixel error of each image as an input. + # Should eventually do that, but no need for now. + [dx, dy, _dxerr, _dyerr] = chi2_shift(self.frames[0], frame) shifted = shift2d(frame, -1 * dx, -1 * dy) shifted_data.append(shifted) @@ -343,21 +357,22 @@ def stack(self): def crop_final(self, bw): """ - Applies crop to the final image. - + Apply crop to the final image. + Parameters ---------- - bw : int, required. + bw : int, required. border width """ 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 - + Apply a crop to the final image. + + Attempts to find bright source in frame using a center-of-mass approach. + Parameters ---------- wx : int, required. @@ -367,42 +382,48 @@ def com_crop_final(self, wx, wy=None, cutoff=None): 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 + 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()]) + 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 + data_to_calc[data_to_calc < 10 * cutoff] = 0.0 if np.sum(data_to_calc) == 0.0: # prevents failures when low SNR - com = (self.final.shape[0]/2, self.final.shape[1]/2) + com = (self.final.shape[0] / 2, self.final.shape[1] / 2) else: 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)] + + 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): """ - Plot the individual frames any step in the process - + Plot the individual frames any step in the process. + Parameters ---------- png_file : str or None, optional. @@ -410,7 +431,7 @@ def plot_frames(self, png_file=None, figsz=3): if None, image not saved """ n = len(self.frames) - fig, axes = plt.subplots(1, n, figsize=(figsz * n, figsz+1)) + fig, axes = plt.subplots(1, n, figsize=(figsz * n, figsz + 1)) for i in range(n): plotframe = self.frames[i] vmax = np.nanmax(medfilt(plotframe, kernel_size=7)) @@ -419,7 +440,7 @@ def plot_frames(self, png_file=None, figsz=3): else: ax = axes ax.imshow(plotframe, origin="lower", vmin=0, vmax=vmax) - ax.set_title("Frame %d" % i) + ax.set_title(f"Frame {i}") ax.set_xticks([]) ax.set_yticks([]) if png_file is not None: @@ -428,6 +449,8 @@ def plot_frames(self, png_file=None, figsz=3): def plot_final(self, show=True, png_file=None): """ + Plot the final product. + Parameters ---------- show : bool, optional. @@ -439,7 +462,7 @@ def plot_final(self, show=True, png_file=None): fig, ax = plt.subplots(1, 1, figsize=(8, 8)) try: cmap = get_colormap(self.target.split(" ")[0]) - except: + except: # noqa: E722 print(f"No custom colormap defined for target {self.target}, setting to default") cmap = cm.viridis ax.imshow(self.final, cmap=cmap, origin="lower") @@ -454,11 +477,11 @@ def plot_final(self, show=True, png_file=None): def write_frames(self, outfiles): """ - Write each individual frame to .fits at any step in the process - + Write each individual frame to .fits at any step in the process. + Parameters ---------- - outfiles : list, required. + outfiles : list, required. filenames to write """ for i in range(len(self.frames)): @@ -470,11 +493,11 @@ def write_frames(self, outfiles): def write_final(self, outfile): """ - Writes the final stacked frame to .fits - + Write the final stacked frame to .fits. + Parameters ---------- - outfile : str, required. + outfile : str, required. filename to write png_file : str or None, optional. filename to which image is saved @@ -487,20 +510,19 @@ def write_final(self, outfile): class Nod(Observation): - """ - simple on-off nod dither pattern - """ + """Simple on-off nod dither pattern.""" def __init__(self, data_fname, sky_fname, instrument): """ + Initialize a nod-type observation. + Parameters ---------- - data_fname : str, required. + data_fname : str, required. filename of input .fits data frame - sky_fname : str, required. + sky_fname : str, required. filename of input .fits sky frame """ - super().__init__(data_fname, instrument) self.sky = Image(sky_fname).data @@ -510,52 +532,55 @@ def __init__(self, data_fname, sky_fname, instrument): self.sky = crop_center(self.sky, self.subc) def apply_sky(self): - """ - simply subtract the sky frame from the data - """ + """Simply subtract the sky frame from the data.""" self.frames = np.array([data - self.sky for data in self.frames]) def uranus_crop(self, bw): """ - Custom crop for Uranus data. We use the right side of the + Apply a custom crop for Uranus data. + + We use the right side of the NIRC2 detector to avoid the bad pixels in the lower left corner so just cutting off the left bit here - + Parameters ---------- bw : int, required. width of the crop """ - szx, szy = self.frames[0].shape[0], self.frames[0].shape[1] + szx, _szy = self.frames[0].shape[0], self.frames[0].shape[1] self.frames = np.array([data[bw : szx - bw, 2 * bw :] for data in self.frames]) class Bxy3(Observation): """ - The famous bxy3 dither at Keck, - which avoids the noisier lower left quadrant of the detector - input fnames MUST be in the conventional order + The famous bxy3 dither at Keck. + + bxy3 avoids the noisier lower left quadrant of the detector + input fnames MUST be in the conventional order where target is in positions [top left, bottom right, top right] """ def __init__(self, fnames, instrument): - ''' + """ + Initialize a bxy3 observation. + Parameters ---------- - fnames : list, required. + fnames : list, required. fits files representing a Keck bxy3 dither instrument : str, required. - ''' - + """ super().__init__(fnames, instrument) def make_sky(self, outfile): """ Make a sky frame via simple median-average of the frames. + Writes fits file containing the sky frame with header info same as zeroth input bxy3 file - + Parameters ---------- outfile : str, required. @@ -565,29 +590,29 @@ def make_sky(self, outfile): sky = np.median(self.frames, axis=0) # change some header info and write to .fits hdulist_out = self.dummy_fits.hdulist - hdulist_out[0].header[self.header_kw_dict['object']] = "SKY_FULL" + hdulist_out[0].header[self.header_kw_dict["object"]] = "SKY_FULL" hdulist_out[0].data = sky hdulist_out[0].writeto(outfile, overwrite=True) def apply_sky(self, fname): """ - identify patches of "normal" sky in each of the three bxy3 frames. - comes up with a median average value of background for each frame. - Then you normalize the "full" sky frame we got from median of bxy3 + Identify patches of "normal" sky in each of the three bxy3 frames. + + Comes up with a median average value of background for each frame. + Then you normalize the "full" sky frame we got from median of bxy3 to the value of the sky brightness in the single frame. Also normalizes everything to be around 1, so can directly multiply by data This ONLY works if bxy3 frames are input in the correct (usual) order, top left | bottom right | top right - + Parameters ---------- - fname : str, required. + fname : str, required. filename pointing to sky frame """ - - full_sky = Image( - fname - ).data # could just as easily use self.sky, but this way could theoretically load a different sky background if desired + # could just as easily use self.sky, but this way could theoretically + # load a different sky background if desired + full_sky = Image(fname).data if full_sky.shape[0] != self.subc: full_sky = crop_center(full_sky, self.subc) @@ -626,11 +651,15 @@ def apply_sky(self, fname): frames_skysub.append(frame - sky_norm) self.frames = np.asarray(frames_skysub) else: - print(f"Sky subtraction not applied to {self.target} {self.filter} {self.subc} (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): """ - Clips each image in bxy3 to its own quadrant. + Clip each image in bxy3 to its own quadrant. + Relies on frames input being in correct order. Should be applied just before stacking """ @@ -648,27 +677,31 @@ def trim(self): class DitherN(Observation): - ''' - Generic N-position dither; fnames can have arbitrary length + """ + Generic N-position dither. + + fnames can have arbitrary length; will make sky by simple median average - ''' - + """ + def __init__(self, fnames, instrument): - ''' + """ + Make a DitherN observation. + Parameters ---------- - fnames : list, required. + fnames : list, required. fits files representing a Keck bxy3 dither instrument : str, required. - ''' - + """ super().__init__(fnames, instrument) - + def make_sky(self, outfile): """ - Make a sky frame via simple median-average of the frames + Make a sky frame via simple median-average of the frames. + Identical to Bxy3.make_sky() - + Parameters ---------- outfile : str, required. @@ -678,20 +711,19 @@ def make_sky(self, outfile): sky = np.median(self.frames, axis=0) # change some header info and write to .fits hdulist_out = self.dummy_fits.hdulist - hdulist_out[0].header[self.header_kw_dict['object']] = "SKY_FULL" + hdulist_out[0].header[self.header_kw_dict["object"]] = "SKY_FULL" hdulist_out[0].data = sky hdulist_out[0].writeto(outfile, overwrite=True) - - + def apply_sky(self, fname): """ - simply subtract the sky frame from the data - + Simply subtract the sky frame from the data. + Parameters ---------- fname : str, required. filename representing sky frame fits """ - sky = Image(fname).data - frames = np.array([data - sky for data in self.frames]) - \ No newline at end of file + pass + # sky = Image(fname).data + # frames = np.array([data - sky for data in self.frames]) diff --git a/nirc2_reduce/prettycolors.py b/nirc2_reduce/prettycolors.py index 15ff50d..03c21db 100755 --- a/nirc2_reduce/prettycolors.py +++ b/nirc2_reduce/prettycolors.py @@ -3,12 +3,14 @@ def make_colormap(seq): """ + Make a custom colormap. + Parameters ---------- - seq : + seq : a sequence of floats and RGB-tuples. The floats should be increasing and in the interval (0,1). - + Returns ------- LinearSegmentedColormap @@ -27,8 +29,8 @@ def make_colormap(seq): def get_colormap(objname): """ - Custom colormaps for making pretty pictures for Twilight Zone website - + Create custom colormaps for pretty pictures for Twilight Zone website. + Parameters ---------- objname : str, required. name of target. currently supported: @@ -38,9 +40,7 @@ def get_colormap(objname): objmap = "gist_heat" elif objname.lower() == "neptune": c = mcolors.ColorConverter().to_rgb - rvb = make_colormap( - [c("black"), c("blue"), 0.33, c("blue"), c("white"), 0.66, c("white")] - ) + rvb = make_colormap([c("black"), c("blue"), 0.33, c("blue"), c("white"), 0.66, c("white")]) objmap = rvb elif objname.lower() == "titan": c = mcolors.ColorConverter().to_rgb diff --git a/nirc2_reduce/sort_rawfiles.py b/nirc2_reduce/sort_rawfiles.py index a6337fb..5d55953 100755 --- a/nirc2_reduce/sort_rawfiles.py +++ b/nirc2_reduce/sort_rawfiles.py @@ -1,11 +1,13 @@ -import warnings import glob +import importlib +import warnings + import numpy as np +import yaml from astropy.io import fits from astropy.table import Table + import nirc2_reduce.data.header_kw_dicts as inst_info -import importlib -import yaml warnings.filterwarnings( "ignore", @@ -17,16 +19,14 @@ ) -def dfits_fitsort( - input_wildcard, - fits_kws=[ "OBJECT", "DATE-OBS", "FILTER"] - ): - """ - Python implementation of the dfits | fitsort bash script workflow +def dfits_fitsort(input_wildcard, fits_kws=None): + r""" + Python implementation of the dfits | fitsort bash script workflow. + Searches headers of all fits in input_dir to find the requested keywords makes a table of filename, kw0, kw1, ... - + Parameters ---------- input_wildcard : str, required. path to filenames @@ -35,23 +35,33 @@ def dfits_fitsort( example 'raw/2023jul14/\*.fits' fits_kws : list, optional. default ['OBJECT', 'DATE-OBS', 'FILTER'] which header keywords to include in the output table - + Returns ------- astropy.table.Table filenames and values corresponding to fits_kws """ + if fits_kws is None: + fits_kws = ["OBJECT", "DATE-OBS", "FILTER"] fnames = glob.glob(input_wildcard) fnames = np.sort(fnames) with warnings.catch_warnings(): warnings.simplefilter("ignore") hdrs = [fits.getheader(f, 0, ignore_missing_end=True) for f in fnames] vals = [ - [fnames[i],] + [hdr[key] for key in fits_kws] for i, hdr in enumerate(hdrs) + [ + fnames[i], + ] + + [hdr[key] for key in fits_kws] + for i, hdr in enumerate(hdrs) ] - names = ["FILENAME",] + fits_kws - dtype = [str,] * len(names) + names = [ + "FILENAME", + ] + fits_kws + dtype = [ + str, + ] * len(names) tab = Table(np.array(vals), names=names, dtype=dtype) return tab @@ -59,13 +69,15 @@ def dfits_fitsort( def split_by_kw(tab, kw): """ + Split a single table into multiple tables according to unique values of kwd. + Parameters ---------- tab : astropy Table, required. table of fits header params output by dfits_fitsort() kw : str, required. table header column to split fnames by - + Returns ------- list @@ -79,70 +91,71 @@ def split_by_kw(tab, kw): return unq.data, split_tables -def get_flats( - tab, - instrument, - ignore_objects=['',] - ): +def get_flats(tab, instrument, ignore_objects=None): """ - scrub table from dfits_fitsort to find domeflaton, domeflatoff filenames - default params are for NIRC2 - + Scrub table from dfits_fitsort to find domeflaton, domeflatoff filenames. + Parameters ---------- tab : astropy Table, required. table of fits header params output by dfits_fitsort() must have FILENAME kw instrument : str, required. - name of instrument used, e.g. nirc2. + name of instrument used, e.g. nirc2. will look for file data/{instrument}.yaml in order to scrub header keywords ignore_objects : list, optional, default empty list. special values in header[object] to ignore. - + Returns ------- - list + list filenames for dome flat off list filenames for dome flat on """ + if ignore_objects is None: + ignore_objects = [ + "", + ] instrument = instrument.lower() with importlib.resources.open_binary(inst_info, f"{instrument}.yaml") as file: yaml_bytes = file.read() header_kw_dict = yaml.safe_load(yaml_bytes) # ignore pre-existing bad pixel maps and master flats - rm_bool = np.sum(np.array([tab[header_kw_dict['object']] == arg for arg in ignore_objects]), axis=0) + rm_bool = np.sum( + np.array([tab[header_kw_dict["object"]] == arg for arg in ignore_objects]), axis=0 + ) rm_i = list(np.argwhere(rm_bool).flatten()) tab.remove_rows(rm_i) # find dome position - targnames, tabs = split_by_kw(tab, header_kw_dict['isdome']['kw']) - domebool = [s.startswith(header_kw_dict['isdome']['yesdome'].strip()) for s in targnames] + targnames, tabs = split_by_kw(tab, header_kw_dict["isdome"]["kw"]) + 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, 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']) + 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, flattabs = split_by_kw(flattab, header_kw_dict['lamps']['kw']) - onbool = [s.startswith(header_kw_dict['lamps']['lampon'].strip()) for s in onoff] + 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!') + 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]] + + 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_flats.py b/nirc2_reduce/tests/test_flats.py index 550edda..fb65338 100644 --- a/nirc2_reduce/tests/test_flats.py +++ b/nirc2_reduce/tests/test_flats.py @@ -1,48 +1,48 @@ -from pytest import fixture -from .. import flats -import numpy as np import os + +import numpy as np from astropy.io import fits +from pytest import fixture + +from nirc2_reduce import flats # until packaging is better and rootdir is found, can use # pytest --rootdir=/Users/emolter/Python/nirc2_reduce test_flats.py @fixture -def datadir(request, tmpdir): +def datadir(request): rootdir = request.config.rootdir path = os.path.join(rootdir, "nirc2_reduce", "tests", "data") return path @fixture -def rawdir(request, tmpdir): +def rawdir(request): rootdir = request.config.rootdir path = os.path.join(rootdir, "nirc2_reduce", "tests", "data", "raw") return path def test_datafiles_working(rawdir, datadir): - assert os.path.isfile(os.path.join(rawdir, "off0.fits")) assert os.path.isfile(os.path.join(datadir, "flat_expected.fits")) def test_flats(datadir, rawdir): """ - fnames required in data/raw: + Fnames required in data/raw: f'off{i}.fits') for i in range(5) f'on{i}.fits') for i in range(5) fnames required in data/ flat_expected.fits badpx_map_expected.fits """ - # make the flat domeflatoff = [os.path.join(rawdir, f"off{i}.fits") for i in range(5)] domeflaton = [os.path.join(rawdir, f"on{i}.fits") for i in range(5)] flat = flats.Flats(domeflatoff, domeflaton) - assert (not np.any(np.isnan(flat.flat))) + assert not np.any(np.isnan(flat.flat)) flat.write(os.path.join(datadir, "flat_test.fits")) # make the badpx map diff --git a/nirc2_reduce/tests/test_multi_reduce.py b/nirc2_reduce/tests/test_multi_reduce.py index db076fc..55aa63f 100644 --- a/nirc2_reduce/tests/test_multi_reduce.py +++ b/nirc2_reduce/tests/test_multi_reduce.py @@ -1,60 +1,52 @@ -import pytest -from pytest import fixture -import importlib -from .. import multi_reduce -import numpy as np -import os import glob +import os + +import numpy as np from astropy.io import fits -import warnings +from pytest import fixture + +from nirc2_reduce import multi_reduce @fixture -def datadir(request, tmpdir): +def datadir(request): rootdir = request.config.rootdir path = os.path.join(rootdir, "nirc2_reduce", "tests", "data") return path @fixture -def rawdir(request, tmpdir): +def rawdir(request): rootdir = request.config.rootdir path = os.path.join(rootdir, "nirc2_reduce", "tests", "data", "raw") return path @fixture -def reddir(request, tmpdir): +def reddir(request): rootdir = request.config.rootdir path = os.path.join(rootdir, "nirc2_reduce", "tests", "data", "reduced") return path @fixture -def flatdir(request, tmpdir): +def flatdir(request): rootdir = request.config.rootdir path = os.path.join(rootdir, "nirc2_reduce", "tests", "data", "flats") return path -def test_MultiBxy3(datadir, rawdir, reddir, flatdir): - +def test_multi_bxy3(datadir, rawdir, reddir, flatdir): obs = multi_reduce.MultiBxy3(rawdir, "nirc2_pre_oct23") obs.process_flats(flatdir) # looks for raw flat frames in rawdir, puts into flatdir obs.run(reddir, flatdir) # test major output files were written correctly - stack_test = fits.open( - os.path.join(reddir, "2017-07-25_Neptune_stacked_nophot_H.fits") - )[0].data - stack_expected = fits.open(os.path.join(datadir, "bxy3_stack_expected.fits"))[ - 0 - ].data + stack_test = fits.open(os.path.join(reddir, "2017-07-25_Neptune_stacked_nophot_H.fits"))[0].data + stack_expected = fits.open(os.path.join(datadir, "bxy3_stack_expected.fits"))[0].data assert np.allclose(stack_test, stack_expected, rtol=1e-3) - flat_test = fits.open(os.path.join(flatdir, "2017-07-25_flat_master_1024_h.fits"))[ - 0 - ].data + 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) diff --git a/nirc2_reduce/tests/test_observation.py b/nirc2_reduce/tests/test_observation.py index a5d847b..1a2e980 100644 --- a/nirc2_reduce/tests/test_observation.py +++ b/nirc2_reduce/tests/test_observation.py @@ -1,34 +1,33 @@ +import os + +import numpy as np import pytest +from astropy.io import fits from pytest import fixture -from .. import observation -import numpy as np -import os -from astropy.io import fits +from nirc2_reduce import observation @fixture -def datadir(request, tmpdir): +def datadir(request): rootdir = request.config.rootdir path = os.path.join(rootdir, "nirc2_reduce", "tests", "data") return path @fixture -def rawdir(request, tmpdir): +def rawdir(request): rootdir = request.config.rootdir path = os.path.join(rootdir, "nirc2_reduce", "tests", "data", "raw") return path def test_datafiles_working(datadir, rawdir): - assert os.path.isfile(os.path.join(rawdir, "bxy3_1.fits")) assert os.path.isfile(os.path.join(datadir, "badpx_map_expected.fits")) def test_bxy3(datadir, rawdir): - # run the bxy3 workflow fnames = [ os.path.join(rawdir, "bxy3_1.fits"), @@ -55,9 +54,7 @@ def test_bxy3(datadir, rawdir): # test the final frame stack_test = fits.open(os.path.join(datadir, "bxy3_stack_test.fits"))[0].data - stack_expected = fits.open(os.path.join(datadir, "bxy3_stack_expected.fits"))[ - 0 - ].data + stack_expected = fits.open(os.path.join(datadir, "bxy3_stack_expected.fits"))[0].data assert np.allclose(stack_test, stack_expected, rtol=1e-3) ### cleanup: remove test fits @@ -66,13 +63,14 @@ def test_bxy3(datadir, rawdir): def test_nod(datadir, rawdir): - obs = observation.Nod( - os.path.join(rawdir, "bxy3_2.fits"), os.path.join(datadir, "sky_expected.fits"), "nirc2_pre_oct23" + os.path.join(rawdir, "bxy3_2.fits"), + os.path.join(datadir, "sky_expected.fits"), + "nirc2_pre_oct23", ) obs.apply_sky() obs.apply_flat(os.path.join(datadir, "flat_expected.fits")) - assert (not np.any(np.isnan(obs.frames))) + assert not np.any(np.isnan(obs.frames)) obs.apply_badpx_map(os.path.join(datadir, "badpx_map_expected.fits")) obs.dewarp() obs.remove_cosmic_rays() @@ -90,11 +88,12 @@ def test_nod(datadir, rawdir): def test_raises(datadir, rawdir): - obs = observation.Nod( - os.path.join(rawdir, "bxy3_2.fits"), os.path.join(datadir, "sky_expected.fits"), "nirc2" + os.path.join(rawdir, "bxy3_2.fits"), + os.path.join(datadir, "sky_expected.fits"), + "nirc2", ) obs.frames = np.empty(obs.frames.shape) obs.frames[:] = np.nan - with pytest.raises(ValueError) as e_info: - obs.remove_cosmic_rays() \ No newline at end of file + with pytest.raises(ValueError): + obs.remove_cosmic_rays() diff --git a/nirc2_reduce/tests/test_sort_rawfiles.py b/nirc2_reduce/tests/test_sort_rawfiles.py index 0ef3f82..c2948e1 100644 --- a/nirc2_reduce/tests/test_sort_rawfiles.py +++ b/nirc2_reduce/tests/test_sort_rawfiles.py @@ -1,31 +1,32 @@ -from pytest import fixture -from .. import sort_rawfiles -import numpy as np - import os -from astropy import table import warnings +import numpy as np +from astropy import table +from pytest import fixture + +from nirc2_reduce import sort_rawfiles + @fixture -def rawdir(request, tmpdir): +def rawdir(request): rootdir = request.config.rootdir path = os.path.join(rootdir, "nirc2_reduce", "tests", "data", "raw") return path def test_dfits_fitsort(rawdir): - with warnings.catch_warnings(): warnings.simplefilter("ignore") input_wildcard = os.path.join(rawdir, "o*.fits") - tab = sort_rawfiles.dfits_fitsort(input_wildcard, - ["OBJECT", "DATE-OBS", "FILTER", "AXESTAT", "FLSPECTR", "SHRNAME"]) + tab = sort_rawfiles.dfits_fitsort( + input_wildcard, ["OBJECT", "DATE-OBS", "FILTER", "AXESTAT", "FLSPECTR", "SHRNAME"] + ) flatoff, flaton = sort_rawfiles.get_flats(tab, "nirc2_pre_oct23") offs_expected = np.array([os.path.join(rawdir, f"off{i}.fits") for i in range(5)]) ons_expected = np.array([os.path.join(rawdir, f"on{i}.fits") for i in range(5)]) - assert type(tab) == table.Table + assert isinstance(tab, table.Table) assert np.all(np.sort(flatoff) == offs_expected) assert np.all(np.sort(flaton) == ons_expected) diff --git a/pyproject.toml b/pyproject.toml index 7a846f1..a9cf8be 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,6 +29,7 @@ dynamic = ["version"] test = [ "pytest", "pytest-cov", + "pre-commit" ] docs = [ "sphinx",