From 18a63972fff73db85dba6bba3afd3b45f192d141 Mon Sep 17 00:00:00 2001 From: Russell Ryan Date: Thu, 23 Jan 2025 12:10:20 -0500 Subject: [PATCH] quick test with ROsteen & BKuhn on WR96 (#83) * quick test with ROsteen & BKuhn on WR96 * Bump the actions group across 1 directory with 2 updates Bumps the actions group with 2 updates in the /.github/workflows directory: [codecov/codecov-action](https://github.com/codecov/codecov-action) and [OpenAstronomy/github-actions-workflows](https://github.com/openastronomy/github-actions-workflows). Updates `codecov/codecov-action` from 4.6.0 to 5.1.2 - [Release notes](https://github.com/codecov/codecov-action/releases) - [Changelog](https://github.com/codecov/codecov-action/blob/main/CHANGELOG.md) - [Commits](https://github.com/codecov/codecov-action/compare/b9fd7d16f6d7d1b5d2bec1a2887e65ceed900238...1e68e06f1dbfde0e4cefc87efeba9e4643565303) Updates `OpenAstronomy/github-actions-workflows` from 1.13.0 to 1.15.0 - [Release notes](https://github.com/openastronomy/github-actions-workflows/releases) - [Commits](https://github.com/openastronomy/github-actions-workflows/compare/924441154cf3053034c6513d5e06c69d262fb9a6...9f1f43251dde69da8613ea8e11144f05cdea41d5) --- updated-dependencies: - dependency-name: codecov/codecov-action dependency-type: direct:production update-type: version-update:semver-major dependency-group: actions - dependency-name: OpenAstronomy/github-actions-workflows dependency-type: direct:production update-type: version-update:semver-minor dependency-group: actions ... Signed-off-by: dependabot[bot] * fix codestyle updates * quick test with ROsteen & BKuhn on WR96 * fix codestyle updates * Remove 3.10 CI test --------- Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Co-authored-by: Ricky O'Steen --- .github/workflows/ci_tests.yml | 6 - pyproject.toml | 2 +- .../core/modules/extract/single/single.py | 175 ++++++++++-------- slitlessutils/core/modules/module.py | 2 +- .../core/modules/tabulate/tabulate.py | 5 +- .../core/preprocess/background/background.py | 107 ++++++----- slitlessutils/core/wfss/config/flatfield.py | 5 +- .../core/wfss/config/instrumentconfig.py | 6 +- slitlessutils/core/wfss/config/sensitivity.py | 15 +- slitlessutils/core/wfss/config/wfssconfig.py | 5 +- slitlessutils/core/wfss/data/wfss.py | 14 +- slitlessutils/examples/wr96.py | 53 ++++-- tox.ini | 2 +- 13 files changed, 236 insertions(+), 161 deletions(-) diff --git a/.github/workflows/ci_tests.yml b/.github/workflows/ci_tests.yml index b004e31b..6cdea9c2 100644 --- a/.github/workflows/ci_tests.yml +++ b/.github/workflows/ci_tests.yml @@ -30,12 +30,6 @@ jobs: strategy: matrix: include: - - os: ubuntu-latest - python: '3.10' - tox_env: 'py310-test' - allow_failure: false - prefix: '' - - os: ubuntu-latest python: '3.11' tox_env: 'py311-test-cov' diff --git a/pyproject.toml b/pyproject.toml index 888dfb50..7b6ac1f5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,7 +31,7 @@ classifiers = [ 'Topic :: Scientific/Engineering :: Astronomy', ] dynamic = ['version'] -requires-python = '>=3.10' +requires-python = '>=3.11' dependencies = [ 'astropy>=5.0.4', 'astroquery', diff --git a/slitlessutils/core/modules/extract/single/single.py b/slitlessutils/core/modules/extract/single/single.py index deeec531..5ab804e5 100644 --- a/slitlessutils/core/modules/extract/single/single.py +++ b/slitlessutils/core/modules/extract/single/single.py @@ -48,10 +48,10 @@ class Single(Module): to set this. Default is None. outpath : str, optional - A path where output products will be written. If set to None or invalid, - then CWD will be used. If set to valid str, then path will be created - (if possible), then used. If still not valid, then will use CWD. - Default is None. + A path where output products will be written. If set to None or + invalid, then CWD will be used. If set to valid str, then path + will be created (if possible), then used. If still not valid, then + will use CWD. Default is None. writecsv : bool, optional Flag to write light-weight boolean files. Default is True. @@ -275,7 +275,7 @@ def _combine(self, results, data, sources, **kwargs): func = np.array(res['func']) wave = np.array(res['wave']) cont = np.array(res['cont']) - file = np.array(res['file']) + nimg = np.array(res['file']) # find the spectral bins that these elements belong to lamb = pars.indices(wave) @@ -290,7 +290,7 @@ def _combine(self, results, data, sources, **kwargs): cont = cont[g] # contamination model wave = wave[g] # wavelengths in A lamb = lamb[g] # wavelength indices - file = file[g] + nimg = nimg[g] # get reverse elements ri = indices.reverse(lamb) @@ -390,7 +390,7 @@ def apply_bitmask(dqa, *args, bitmask=None): out = tuple(a[g] for a in args) else: LOGGER.warning(f'Bitmask ({bitmask}) removes all pixels') - out = tuple([] for a in args) + out = tuple(np.empty(0, dtype=a.dtype) for a in args) if len(args) == 1: out = out[0] return out @@ -446,7 +446,11 @@ def extract(self, data, sources, **kwargs): # sort out optional inputs cartesian = kwargs.get('cartesian', True) - profile = kwargs.get('profile', 'uniform').lower() + profile = kwargs.get('profile', 'uniform') + if isinstance(profile, str): + profile = profile.lower() + else: + profile = 'none' # padding for the contamination models padx = (5, 5) @@ -484,7 +488,6 @@ def extract(self, data, sources, **kwargs): time = phdr['EXPTIME'] sci /= time unc /= time - # initialize the contamination modeling self.contamination.make_model(sources, h5) @@ -537,7 +540,7 @@ def extract(self, data, sources, **kwargs): if cartesian: # for making the residuals - # mod = np.zeros_like(sci) + mod = np.zeros_like(sci) for x in range(x0, x1 + 1, 1): g = np.where(x == xg)[0] xx = xg[g] @@ -549,75 +552,93 @@ def extract(self, data, sources, **kwargs): # apply the bitmask xx, yy, vv, ww, dw = self.apply_bitmask( dqa[yy, xx], xx, yy, vv, ww, dw, bitmask=bitmask) - - ww, _y = indices.decimate(ww * vv, yy) - vv, yy = indices.decimate(vv, yy) - ww /= vv - xx = np.full_like(yy, x, dtype=int) - - # the average wavelength in this column - wave = np.average(ww, weights=vv) - - # compute the calibrations - disp = detdata.config[ordname].dispersion( - *source.xyc, wavelength=ww) - sens = detdata.config[ordname].sensitivity(ww) - flat = flatfield(xx, yy, ww) - area = detdata.relative_pixelarea(xx, yy) - den = flat * area * sens * fluxscale * disp - - # apply calibrations - ss = sci[yy, xx] / den - uu = unc[yy, xx] / den - - # set the cross dispersion profile - if profile == 'uniform': - # simple summing over pixels - flam = np.sum(ss) - func = np.sqrt(np.sum(uu**2)) - else: - # weighting by profile (a la Horne) - if profile == 'forward': - prof = vv.copy() - elif profile == 'data': - prof = np.maximum(sci[yy, x], 0.) + if len(xx) > 0: + ww, _y = indices.decimate(ww * vv, yy) + vv, yy = indices.decimate(vv, yy) + ww /= vv + xx = np.full_like(yy, x, dtype=int) + + # the average wavelength in this column + wave = np.average(ww, weights=vv) + + # compute the calibrations + disp = detdata.config[ordname].dispersion( + *source.xyc, wavelength=ww) + sens = detdata.config[ordname].sensitivity(ww) + flat = flatfield(xx, yy, ww) + area = detdata.relative_pixelarea(xx, yy) + den = flat * area * sens * fluxscale * disp + + # apply calibrations, but guard against + # divide by zero errors + bad = np.where(den == 0)[0] + den[bad] = 1. + ss = sci[yy, xx] / den + uu = unc[yy, xx] / den + ss[bad] = np.nan + uu[bad] = np.nan + + # set the cross dispersion profile + if profile is None or profile == 'none': + # simple summing over pixels + flam = np.nansum(ss) + func = np.sqrt(np.nansum(uu**2)) + + # set the profile to a dummy value + prof = 1.0 else: - msg = f'Profile setting ({profile}) is invalid' - LOGGER.error(msg) - raise RuntimeError(msg) - - # normalize the profile - prof /= np.sum(prof) - - wht = prof / uu**2 - norm = np.sum(prof * wht) - flam = np.sum(ss * wht) / norm - func = np.sqrt(np.sum(prof) / norm) - - # this can happen if sens==0 - if np.isnan(flam): - flam = 0.0 - - # update the model - # mod[yy, x] = flam*den*prof - - # compute contamination - if self.contamination: - cc = chdu.data[yy - yoff, xx - xoff] / den - if profile == 'uniform': - cont = np.sum(cc) + # if using a profile + if profile == 'forward': + prof = vv.copy() + elif profile == 'data': + prof = np.maximum(sci[yy, x], 0.) + elif profile == 'uniform': + prof = np.ones_like(ss, dtype=float) + else: + msg = f'Profile setting ({profile}) is invalid' + LOGGER.error(msg) + raise RuntimeError(msg) + + # normalize the profile + profnorm = np.nansum(prof) + prof /= profnorm + + wht = prof / uu**2 + norm = np.nansum(prof * wht) + + # take out some protection from + # divide by zero + if norm == 0: + flam = np.nan + func = np.nan + else: + flam = np.nansum(ss * wht) / norm + func = np.sqrt(profnorm / norm) + + # update the model + mod[yy, x] = flam * den * prof + + # compute contamination + if self.contamination: + cc = chdu.data[yy - yoff, xx - xoff] / den + if profile == 'uniform': + cont = np.sum(cc) + else: + cont = np.sum(cc * wht) / norm else: - cont = np.sum(cc * wht) / norm - else: - cont = 0.0 - - # save the results - results[segid]['flam'].append(flam) - results[segid]['func'].append(func) - results[segid]['wave'].append(wave) - results[segid]['dwav'].append(0.0) - results[segid]['cont'].append(cont) - + cont = 0.0 + + # save the results + results[segid]['flam'].append(flam) + results[segid]['func'].append(func) + results[segid]['wave'].append(wave) + results[segid]['dwav'].append(0.0) + results[segid]['cont'].append(cont) + with open(f'{data.dataset}_{segid}.sed', 'w') as fp: + for args in zip(results[segid]['wave'], results[segid]['flam']): + print(*args, file=fp) + + # fits.writeto(f'{data.dataset}_res.fits', (sci-mod)/unc, overwrite=True) else: di = 0.5 diff --git a/slitlessutils/core/modules/module.py b/slitlessutils/core/modules/module.py index 1b350161..64b5a114 100644 --- a/slitlessutils/core/modules/module.py +++ b/slitlessutils/core/modules/module.py @@ -50,7 +50,7 @@ class Module: """ - def __init__(self, func, path='tables', ncpu=None, postfunc=None, + def __init__(self, func, path='su_tables', ncpu=None, postfunc=None, multiprocess=True, **kwargs): self.multiprocess = multiprocess diff --git a/slitlessutils/core/modules/tabulate/tabulate.py b/slitlessutils/core/modules/tabulate/tabulate.py index 5685c7a1..18bf1763 100644 --- a/slitlessutils/core/modules/tabulate/tabulate.py +++ b/slitlessutils/core/modules/tabulate/tabulate.py @@ -37,9 +37,10 @@ class Tabulate(Module): """ + # should be +0.5 for uvis # define the pixel footprint - DX = np.array([0, 0, 1, 1], dtype=float) - 0.5 - DY = np.array([0, 1, 1, 0], dtype=float) - 0.5 + DX = np.array([0, 0, 1, 1], dtype=float) + 0.5 + DY = np.array([0, 1, 1, 0], dtype=float) + 0.5 DESCRIPTION = 'Tabulating' diff --git a/slitlessutils/core/preprocess/background/background.py b/slitlessutils/core/preprocess/background/background.py index 76e35e35..35d70221 100644 --- a/slitlessutils/core/preprocess/background/background.py +++ b/slitlessutils/core/preprocess/background/background.py @@ -48,16 +48,20 @@ def wrapper(self, filename, newfile=None, inplace=False, **kwargs): LOGGER.warning(msg) return - # check that this is *NOT* a subarray. Eventually, we can excise - # the master sky image(s) and this check can be removed, but that is - # going to take a little more care. Specifically, the concerns will be - # which CCD/detector is the subarray on? This affects the primary - # for-loop over the HDUL. This shouldn't be too hard. + # check that this is *NOT* a subarray. Eventually, we can + # excise the master sky image(s) and this check can be + # removed, but that is going to take a little more care. + # Specifically, the concerns will be which CCD/detector is + # the subarray on? This affects the primary for-loop + # over the HDUL. This shouldn't be too hard. if mastersky and phdr.get('SUBARRAY', False): msg = f"Master-sky subtraction is not supported for subarrays: {filename}" LOGGER.knownissue(msg) return + # dummy variable for writing a new file + wrote = False + # look at each HDU for hdu in hdul: @@ -70,55 +74,64 @@ def wrapper(self, filename, newfile=None, inplace=False, **kwargs): hdr = hdul[('SCI', vers)].header unc = hdul[('ERR', vers)].data dqa = hdul[('DQ', vers)].data - - # get an estimate of the sky from statistics - gpx = (dqa == 0) # these are good pixels - ave, med, sig = sigma_clipped_stats(sci[gpx], sigma=self.skysigma) - - # if doing master sky, then need to prep the sky - if mastersky: - if isinstance(backfile, str) and os.path.exists(backfile): - # get the model and check normalization - mod = fits.getdata(backfile, ('SKY', vers)) - - a, m, s = sigma_clipped_stats(mod) - if np.abs(a - 1) > 1e-2: - msg = (f'The master sky image ({backfile}) is ' - f'unnormalized {a}. Results will be fine, ' - f'but the units may not make sense.') + skycorr = hdr.get('SKYCORR', False) + + if not skycorr: + # get an estimate of the sky from statistics + gpx = (dqa == 0) # these are good pixels + ave, med, sig = sigma_clipped_stats(sci[gpx], sigma=self.skysigma) + + # if doing master sky, then need to prep the sky + if mastersky: + if isinstance(backfile, str) and os.path.exists(backfile): + # get the model and check normalization + mod = fits.getdata(backfile, ('SKY', vers)) + + a, m, s = sigma_clipped_stats(mod) + if np.abs(a - 1) > 1e-2: + msg = (f'The master sky image ({backfile}) is ' + f'unnormalized {a}. Results will be fine, ' + f'but the values may be suspect.') + LOGGER.warning(msg) + + # excise if need-be + if phdr['INSTRUME'] == 'WFC3': + x0 = -int(hdr.get('LTV1', 0)) + y0 = -int(hdr.get('LTV2', 0)) + nx = hdr['NAXIS1'] + ny = hdr['NAXIS2'] + mod = mod[y0:ny + y0, x0:nx + x0] + + ret, src = func(self, sci, hdr, unc, med, gpx, mod, **kwargs) + else: + msg = f'The master-sky image is invalid: {backfile}. ' + \ + 'No subtraction performed.' LOGGER.warning(msg) + ret = 0. + else: + ret, src = func(self, sci, hdr, unc, med, gpx, **kwargs) - # excise if need-be - if phdr['INSTRUME'] == 'WFC3': - x0 = -int(hdr.get('LTV1', 0)) - y0 = -int(hdr.get('LTV2', 0)) - nx = hdr['NAXIS1'] - ny = hdr['NAXIS2'] - mod = mod[y0:ny + y0, x0:nx + x0] + # have the model now, so subtract it + hdul[('SCI', vers)].data = sci - ret + hdul[('SCI', vers)].header = hdr - ret, src = func(self, sci, hdr, unc, med, gpx, mod, **kwargs) - else: - msg = f'The master-sky image is invalid: {backfile}. ' + \ - 'No subtraction performed.' - LOGGER.warning(msg) - ret = 0. + # update the dummy + wrote = True else: - ret, src = func(self, sci, hdr, unc, med, gpx, **kwargs) - - # have the model now, so subtract it - hdul[('SCI', vers)].data = sci - ret - hdul[('SCI', vers)].header = hdr + LOGGER.warn('Images are already background subtracted. Ignoring.') # do something for each type of inplace writing, but always # return a filename - if inplace: - outfile = filename + if wrote: + if inplace: + outfile = filename + else: + if newfile is None: + newfile = f'{data.dataset}_sub.fits' + hdul.writeto(newfile, overwrite=True) + outfile = newfile else: - if newfile is None: - newfile = f'{data.dataset}_sub.fits' - hdul.writeto(newfile, overwrite=True) - outfile = newfile - + outfile = filename return outfile return wrapper return background @@ -548,9 +561,9 @@ def update_header(self, hdr): The fits header to update """ + hdr.set('SKYCORR', value=True, comment='Master sky subtracted?') hdr.set('SKYNSIG', value=self.skysigma, comment='Nsigma for constant sky model') - hdr.set('SRCNSIG', value=self.srcsigma, comment='Nsigma for thresholding sources') hdr.set('OUTNSIG', value=self.outliersigma, diff --git a/slitlessutils/core/wfss/config/flatfield.py b/slitlessutils/core/wfss/config/flatfield.py index 70121163..8bb492b8 100644 --- a/slitlessutils/core/wfss/config/flatfield.py +++ b/slitlessutils/core/wfss/config/flatfield.py @@ -291,7 +291,7 @@ def __call__(self, x, y, l): # super().update_header(hdr) -def load_flatfield(*args, unity=False): +def load_flatfield(*args, **kwargs): """ A factory function to load different flat field types @@ -316,6 +316,9 @@ def load_flatfield(*args, unity=False): """ + # sort out the inputs + unity = kwargs.get('unity', False) + n = len(args) if n == 0 or unity: obj = UnityFlatField() diff --git a/slitlessutils/core/wfss/config/instrumentconfig.py b/slitlessutils/core/wfss/config/instrumentconfig.py index 33d7cd57..27665523 100644 --- a/slitlessutils/core/wfss/config/instrumentconfig.py +++ b/slitlessutils/core/wfss/config/instrumentconfig.py @@ -363,8 +363,8 @@ class DetectorConfig: Class to configure a WFSS detector """ - def __init__(self, name, data, disperser, path, exptime=1000., background=0., - **kwargs): + def __init__(self, name, data, disperser, path, exptime=1000., + background=0., **kwargs): """ Initializer method @@ -700,7 +700,7 @@ def from_fitsfile(cls, filename, **kwargs): disperser = (h0['FILTER'], h0['PUPIL']) kwargs['fwcpos'] = h0['FWCPOS'] else: - raise NotImplementedError(f"Unsupported: {tel = } {ins = }") + raise NotImplementedError(f"Unsupported: {tel=} {ins=}") # instantiate the object insconf = cls(tel, ins, disperser, subarray=subarray, **kwargs) diff --git a/slitlessutils/core/wfss/config/sensitivity.py b/slitlessutils/core/wfss/config/sensitivity.py index 95c0dbf3..6a210a45 100644 --- a/slitlessutils/core/wfss/config/sensitivity.py +++ b/slitlessutils/core/wfss/config/sensitivity.py @@ -1,4 +1,6 @@ +import os import numpy as np + from astropy.io import fits from ...photometry import Band @@ -30,11 +32,14 @@ def __init__(self, sensfile, senslimit=1e10): self.senslimit = senslimit # read the file - data, header = fits.getdata(self.sensfile, 1, header=True) - - g = np.where(data['SENSITIVITY'] > self.senslimit) - Band.__init__(self, data['WAVELENGTH'], data['SENSITIVITY'], - where=g, unit=header.get('TUNIT1', '')) + if os.path.exists(self.sensfile): + data, header = fits.getdata(self.sensfile, 1, header=True) + + g = np.where(data['SENSITIVITY'] > self.senslimit) + Band.__init__(self, data['WAVELENGTH'], data['SENSITIVITY'], + where=g, unit=header.get('TUNIT1', '')) + else: + Band.__init__(self, np.array([0, np.inf]), np.array([1., 1.])) def __str__(self): return f'Sensitivity curve: {self.sensfile}' diff --git a/slitlessutils/core/wfss/config/wfssconfig.py b/slitlessutils/core/wfss/config/wfssconfig.py index c35f4425..d21977de 100644 --- a/slitlessutils/core/wfss/config/wfssconfig.py +++ b/slitlessutils/core/wfss/config/wfssconfig.py @@ -161,7 +161,7 @@ def add_shift(self, dx, dy): self.xshift += dx self.yshift += dy - def load_flatfield(self, unity=False): + def load_flatfield(self, **kwargs): """ Method to read a flat field from the configuration file. @@ -176,6 +176,9 @@ def load_flatfield(self, unity=False): The two-dimensional flat field """ + # sort out the defaults + unity = kwargs.get('unity', False) + if unity or (self.ffname is None): ff = load_flatfield() else: diff --git a/slitlessutils/core/wfss/data/wfss.py b/slitlessutils/core/wfss/data/wfss.py index f3f8782f..6fa7a60e 100644 --- a/slitlessutils/core/wfss/data/wfss.py +++ b/slitlessutils/core/wfss/data/wfss.py @@ -1,4 +1,5 @@ import os +import re import warnings import numpy as np @@ -667,9 +668,18 @@ def get_parameters(self): @property def dataset(self): + if isinstance(self.config.suffix, (list, tuple, np.ndarray)): + regex = '|'.join([f'_{s}' for s in self.config.suffix]) + else: + regex = '_' + self.config.suffix + filename = os.path.basename(self.filename) - tokens = filename.split(f'_{self.config.suffix}') - return tokens[0] + tokens = re.split(regex, filename) + dataset = tokens[0] + if len(tokens) == 1: + LOGGER.warning('Filename does not contain suffix.') + + return dataset @property def telescope(self): diff --git a/slitlessutils/examples/wr96.py b/slitlessutils/examples/wr96.py index f2d91899..1f44a017 100644 --- a/slitlessutils/examples/wr96.py +++ b/slitlessutils/examples/wr96.py @@ -4,6 +4,7 @@ import matplotlib.pyplot as plt import numpy as np from astropy.io import fits +from astropy.modeling import models, fitting from astropy.wcs import WCS from astroquery.mast import Observations from drizzlepac import astrodrizzle @@ -42,7 +43,7 @@ # position of source RA = 264.1019010 # in deg DEC = -32.9087315 # in deg -RAD = 0.25 # in arcsec +RAD = 0.5 # in arcsec SUFFIX, DRZSUF = 'flc', 'drc' SCALE = 0.05 # driz image pix scale @@ -91,7 +92,6 @@ def preprocess_direct(): imgfile = f'{imgdset}_{SUFFIX}.fits' # su.core.preprocess.crrej.laplace(imgfile,inplace=True) files.append(imgfile) - # mosaic data via astrodrizzle astrodrizzle.AstroDrizzle(files, output=ROOT, build=False, static=False, skysub=True, driz_separate=False, @@ -109,11 +109,27 @@ def preprocess_direct(): wcs = WCS(hdr) x, y = wcs.all_world2pix(RA, DEC, 0) + xim, yim = np.meshgrid(np.arange(hdr['NAXIS1']), + np.arange(hdr['NAXIS2'])) + + sz = 10 + xmin = int(x) - sz + xmax = int(x) + sz + ymin = int(y) - sz + ymax = int(y) + sz + pix = (slice(ymin, ymax), slice(xmin, xmax)) + + mod = models.Gaussian2D(x_mean=sz, y_mean=sz, amplitude=np.amax(img[pix])) + + fitter = fitting.LevMarLSQFitter() + yy, xx = np.indices(img[pix].shape, dtype=float) + res = fitter(mod, xx, yy, img[pix]) + x = res.x_mean.value + xmin + y = res.y_mean.value + ymin - xx, yy = np.meshgrid(np.arange(hdr['NAXIS1']), - np.arange(hdr['NAXIS2'])) - rr = np.hypot(xx - x, yy - y) + rr = np.hypot(xim - x, yim - y) seg = rr < (RAD / SCALE) + seg = seg.astype(int) # add some things for SU hdr['TELESCOP'] = TELESCOPE @@ -122,7 +138,7 @@ def preprocess_direct(): # write the files to disk fits.writeto(f'{ROOT}_{DRZSUF}_sci.fits', img, hdr, overwrite=True) - fits.writeto(f'{ROOT}_{DRZSUF}_seg.fits', seg.astype(int), hdr, overwrite=True) + fits.writeto(f'{ROOT}_{DRZSUF}_seg.fits', seg, hdr, overwrite=True) def extract_single(): @@ -137,12 +153,12 @@ def extract_single(): zeropoint=ZEROPOINT) # project the sources onto the grism images - tab = su.modules.Tabulate(ncpu=1) + tab = su.modules.Tabulate(ncpu=1, orders=('+1',), remake=True) pdtfiles = tab(data, sources) # noqa: F841 # run the single-orient extraction - ext = su.modules.Single('+1', mskorders=None, root=ROOT) - res = ext(data, sources) # noqa: F841 + ext = su.modules.Single('+1', mskorders=None, root=ROOT, ncpu=1) + res = ext(data, sources, profile='forward') # noqa: F841 def plot_spectra(): @@ -156,7 +172,12 @@ def plot_spectra(): # load and smooth the Larsen spectrum l, f = np.loadtxt(reffile, unpack=True, usecols=(0, 1)) f /= 1e-13 - ff = gaussian_filter1d(f, 40) + + # smooth the Larsen spectrum. the Scale factor is from comparing the + # notional ACS dispersion (40A/pix) to the spectrum quoted from + # Pasquali+ 2002 which has 1.26 A/pix. Therefore smoothing factor is + # 40./1.26 = 31.7 + ff = gaussian_filter1d(f, 31.7) # load the data and change the units dat, hdr = fits.getdata(f'{ROOT}_x1d.fits', header=True) @@ -164,18 +185,22 @@ def plot_spectra(): dat['func'] *= cfg.fluxscale / 1e-13 # get a good range of points to compute a (variance-weighted) scale factor - g = np.where((dat['lamb'] > 5800) & (dat['lamb'] < 9900))[0] + lmin = 7350. + lmax = 8120. + g = np.where((lmin <= dat['lamb']) & (dat['lamb'] <= lmax))[0] ff2 = np.interp(dat['lamb'], l, ff) den = np.sum((dat['flam'][g] / dat['func'][g])**2) num = np.sum((ff2[g] / dat['func'][g]) * (dat['flam'][g] / dat['func'][g])) scl = num / den + # scl = np.median(ff2[g] / dat['flam'][g]) + + print(f'Scaling factor: {scl}') # plot the SU spectrum + plt.axvspan(lmin, lmax, color='lightgrey') plt.plot(dat['lamb'], dat['flam'] * scl, label=GRATING) plt.plot(l, ff, label='Larsen et al. (smoothed)') - print('scaling factor', scl, 1. / scl) - # uncomment this to see the hi-res file. # plt.plot(l, f, label='Larsen et al. (high-res)') @@ -185,7 +210,7 @@ def plot_spectra(): # put on legend and change limits plt.legend(loc='upper left') - plt.ylim(0.05, 0.55) + plt.ylim(0.0, 0.85) plt.xlim(5500, 10000) # write the file to disk diff --git a/tox.ini b/tox.ini index c5fc5e73..20b7970d 100644 --- a/tox.ini +++ b/tox.ini @@ -1,6 +1,6 @@ [tox] envlist = - py{310,311,312}-test{,-devdeps}{,-cov} + py{311,312}-test{,-devdeps}{,-cov} build_docs linkcheck codestyle