diff --git a/python/lsst/atmospec/processStar.py b/python/lsst/atmospec/processStar.py index 132fba3..ffbe428 100644 --- a/python/lsst/atmospec/processStar.py +++ b/python/lsst/atmospec/processStar.py @@ -518,11 +518,11 @@ def validate(self): except ModuleNotFoundError: getObsAtmo = None if uvspecPath is None and getObsAtmo is None and self.doFitAtmosphere is True: - raise FieldValidationError(self.__class__.doFitAtmosphere, self, - "uvspec is not in the path nor getObsAtmo is installed," + raise FieldValidationError(self.__class__.doFitAtmosphere, self, + "uvspec is not in the path nor getObsAtmo is instaled," " but doFitAtmosphere is True.") - if uvspecPath is None and getObsAtmo is None and self.doFitAtmosphereOnSpectrogram is True: - raise FieldValidationError(self.__class__.doFitAtmosphereOnSpectrogram, self, + if uvspecPath is None and getObsAtmo is None and self.doFitAtmosphereOnSpectrogram is True: + raise FieldValidationError(self.__class__.doFitAtmosphereOnSpectrogram, self, "uvspec is not in the path nor getObsAtmo is installed," " the path, but doFitAtmosphere is True.") diff --git a/python/lsst/atmospec/spectraction.py b/python/lsst/atmospec/spectraction.py index fd5d678..c99443d 100644 --- a/python/lsst/atmospec/spectraction.py +++ b/python/lsst/atmospec/spectraction.py @@ -44,7 +44,7 @@ from lsst.daf.base import DateTime # noqa: E402 from .utils import getFilterAndDisperserFromExp # noqa: E402 -from .spectroFlat import getGainDict, getPTCGainDict, makeSensorFlat, getCertifiedFlat, makeGainFlat # noqa: E402 +from .spectroFlat import getPTCGainDict, makeSensorFlat, getCertifiedFlat, makeGainFlat # noqa: E402 class SpectractorShim: @@ -217,8 +217,7 @@ def _translateCentroid(dmXpos, dmYpos): self._setReadNoiseFromExp(image, exp, 8.5) # xxx remove hard coding of 1 below! import lsst.daf.butler as dafButler - # butler = dafButler.Butler("/sdf/group/rubin/repo/oga/", collections=['LATISS/calib','LATISS/raw/all']) - butler = dafButler.Butler("/repo/embargo", collections=['LATISS/calib','LATISS/raw/all']) + butler = dafButler.Butler("/repo/embargo", collections=['LATISS/calib', 'LATISS/raw/all']) PTCGainDict = getPTCGainDict(butler) certifiedFlat = getCertifiedFlat(butler, dataId=exp.visitInfo.id, filter="empty") image.gain = self._transformArrayFromExpToImage(makeGainFlat(certifiedFlat, PTCGainDict).image.array) @@ -282,7 +281,8 @@ def _getImageData(self, exp, trimToSquare=False): return self._transformArrayFromExpToImage(data) def _transformArrayFromExpToImage(self, array): - # apply transformation on an exposure array to have correct orientation in Spectractor Image + # apply transformation on an exposure array to have correct orientation + # in Spectractor Image return array.T[::, ::] def _setReadNoiseFromExp(self, spectractorImage, exp, constValue=None): @@ -293,16 +293,20 @@ def _setReadNoiseFromExp(self, spectractorImage, exp, constValue=None): # TODO: Check with Jeremy if we want the raw read noise # or the per-pixel variance. Either is doable, just need to know. raise NotImplementedError("Setting noise image from exp variance not implemented") - raw = butler.get('raw',instrument='LATISS', exposure=dataId, detector=0, collections=['LATISS/calib','LATISS/raw/all']) - det=raw.getDetector() + raw = butler.get('raw', instrument='LATISS', exposure=dataId, detector=0, + collections=['LATISS/calib', 'LATISS/raw/all']) + det = raw.getDetector() raw_image = raw.getImage() spectractorImage.read_out_noise = np.zeros_like(spectractorImage.data) - for amp_cur in det : - noise_over_serial_fliped=(raw_image[amp_cur.getRawSerialOverscanBBox()].array)[:,:-2].std(axis=0).mean() + for amp in det: + ampData = raw_image[amp.getRawSerialOverscanBBox()].array + noise_over_serial_fliped = ampData[:, :-2].std(axis=0).mean() bbox = amp.getBBox() spectractorImage.read_out_noise[bbox] = noise_over_serial_fliped - print('amp name=',amp_cur.getName(),'header noise ? =',amp_cur.getReadNoise(),'header gain ? =',amp_cur.getGain(),'(flip =',noise_over_serial_fliped,')') - + print('amp name=', amp.getName(), + ' header noise ? =', amp.getReadNoise(), + 'header gain ? =', amp.getGain(), + '(flip =', noise_over_serial_fliped, ')') def _setReadNoiseToNone(self, spectractorImage): spectractorImage.read_out_noise = None @@ -318,7 +322,7 @@ def _setMask(self, image, exp): CR = exp.getMask().getPlaneBitMask("CR") badPixels = np.logical_or((exp.getMask().array & BAD) > 0, (exp.getMask().array & CR) > 0) image.mask = self._transformArrayFromExpToImage(badPixels) - + def _setGainFromExp(self, spectractorImage, exp, constValue=None): # xxx need to implement this properly # Note that this is array-like and per-amplifier diff --git a/python/lsst/atmospec/spectroFlat.py b/python/lsst/atmospec/spectroFlat.py index 0158586..9a59932 100644 --- a/python/lsst/atmospec/spectroFlat.py +++ b/python/lsst/atmospec/spectroFlat.py @@ -6,20 +6,19 @@ from scipy.ndimage import median_filter, uniform_filter, gaussian_filter -# Butler -# import lsst.daf.butler as dafButler -# butler = dafButler.Butler(repo, collections=calibCollections) - - def find_flats(butler, cameraName='LATISS', filter='empty', disperser='empty', obs_type='flat'): registry = butler.registry physical_filter = f'{filter}~{disperser}' - where = f"instrument='{cameraName}' AND physical_filter='{physical_filter}' AND exposure.observation_type='{obs_type}'" + where = ( + f"instrument='{cameraName}' AND physical_filter='{physical_filter}' " + f"AND exposure.observation_type='{obs_type}'" + ) flat_records = list(registry.queryDimensionRecords('exposure', where=where)) if len(flat_records) == 0: raise ValueError( - 'WARNING: No flats found with the selected settings: \n{cameraName=} {physical_filter=} {obs_type=}') + 'WARNING: No flats found with the selected settings: \n' + f'{cameraName=} {physical_filter=} {obs_type=}') flat_dates = np.sort(np.unique([flat_.day_obs for flat_ in flat_records])) ids_ = [flat_.id for flat_ in flat_records] @@ -117,19 +116,27 @@ def makeOpticFlat(certifiedFlat, kernel='mean', window_size=40, mode='mirror', p data /= np.median(data) if kernel == 'gauss': logger.info(f'Pixel flat smoothing with Gaussian filter with {window_size=}.') - # 'ATTENTION: window size should be equal to Gaussian standard deviation (sigma)' - opticFlat[bbox].maskedImage.image.array[:, :] = gaussian_filter(data, sigma=window_size, mode=mode) + # 'ATTENTION: window size should be equal to Gaussian standard + # deviation (sigma)' + opticFlat[bbox].maskedImage.image.array[:, :] = gaussian_filter(data, sigma=window_size, + mode=mode) elif kernel == 'mean': logger.info( - f'Pixel flat smoothing with mean filter.\nMasking outliers <{percentile:.2f} and >{100. - percentile:.2f} percentiles') - mask = (data >= np.percentile(data.flatten(), percentile)) * ( - data <= np.percentile(data.flatten(), 100. - percentile)) + f'Pixel flat smoothing with mean filter.\n' + f'Masking outliers <{percentile:.2f} and >{100. - percentile:.2f} percentiles' + ) + mask = ( + data >= np.percentile(data.flatten(), percentile)) * ( + data <= np.percentile(data.flatten(), 100. - percentile) + ) interp_array = np.copy(data) interp_array[~mask] = np.median(data) - opticFlat[bbox].maskedImage.image.array[:, :] = uniform_filter(interp_array, size=window_size, mode=mode) + opticFlat[bbox].maskedImage.image.array[:, :] = uniform_filter(interp_array, size=window_size, + mode=mode) elif kernel == 'median': logger.info(f'Pixel flat smoothing with median filter with {window_size=}.') - opticFlat[bbox].maskedImage.image.array[:, :] = median_filter(data, size=(window_size, window_size), + opticFlat[bbox].maskedImage.image.array[:, :] = median_filter(data, + size=(window_size, window_size), mode=mode) else: raise ValueError(f'Kernel must be within ["mean", "median", "gauss"]. Got {kernel=}.') @@ -139,34 +146,52 @@ def makeOpticFlat(certifiedFlat, kernel='mean', window_size=40, mode='mirror', p return opticFlat -def makePixelFlat(certifiedFlat, kernel='mean', window_size=40, mode='mirror', percentile=1., remove_background=True): +def makePixelFlat( + certifiedFlat, + kernel='mean', + window_size=40, + mode='mirror', + percentile=1., + remove_background=True +): logger = logging.getLogger() logger.info(f'Window size for {kernel} smoothing = {window_size}') pixelFlat = certifiedFlat.clone() - opticFlat = makeOpticFlat(certifiedFlat, kernel=kernel, window_size=window_size, mode=mode, percentile=percentile) + opticFlat = makeOpticFlat(certifiedFlat, kernel=kernel, window_size=window_size, mode=mode, + percentile=percentile) detector = opticFlat.getDetector() for amp in detector: bbox = amp.getBBox() - pixelFlat[bbox].maskedImage.image.array[:, :] /= np.median(pixelFlat[bbox].maskedImage.image.array[:, :]) - pixelFlat[bbox].maskedImage.image.array[:, :] /= opticFlat[bbox].maskedImage.image.array[:, :] + ampData = pixelFlat[bbox].maskedImage.image.array[:, :] + ampData /= np.median(ampData) + ampData /= opticFlat[bbox].maskedImage.image.array[:, :] if remove_background: from astropy.stats import SigmaClip from photutils.background import Background2D, MedianBackground sigma_clip = SigmaClip(sigma=3.0) bkg_estimator = MedianBackground() - bkg = Background2D(pixelFlat[bbox].maskedImage.image.array[:, :], (20, 20), filter_size=(3, 3), sigma_clip=sigma_clip, bkg_estimator=bkg_estimator) - pixelFlat[bbox].maskedImage.image.array[:, :] /= bkg.background - - - + bkg = Background2D(ampData, (20, 20), + filter_size=(3, 3), + sigma_clip=sigma_clip, + bkg_estimator=bkg_estimator) + ampData /= bkg.background pixelFlat.maskedImage.mask[:] = 0x0 pixelFlat.maskedImage.variance[:] = 0.0 return pixelFlat -def makeSensorFlat(certifiedFlat, gainDict, kernel='mean', window_size=40, mode='mirror', percentile=1., invertGains=False, remove_background=True): +def makeSensorFlat( + certifiedFlat, + gainDict, + kernel='mean', + window_size=40, + mode='mirror', + percentile=1., + invertGains=False, + remove_background=True +): logger = logging.getLogger() logger.info(f'Window size for {kernel} smoothing = {window_size}') pixelFlat = makePixelFlat(certifiedFlat, kernel=kernel, window_size=window_size, mode=mode, diff --git a/python/lsst/atmospec/utils.py b/python/lsst/atmospec/utils.py index 7a2ebd1..46cf2b3 100644 --- a/python/lsst/atmospec/utils.py +++ b/python/lsst/atmospec/utils.py @@ -66,6 +66,7 @@ def getGainDict(exposure): gainDict[amp.getName()] = amp.getGain() return gainDict + def makeGainFlat(exposure, gainDict, invertGains=False): """Given an exposure, make a flat from the gains.