Skip to content

Commit

Permalink
flake8 fixups and f-string bugfix
Browse files Browse the repository at this point in the history
  • Loading branch information
mfisherlevine committed Aug 16, 2024
1 parent 0eff5c2 commit 99b720e
Show file tree
Hide file tree
Showing 4 changed files with 69 additions and 39 deletions.
8 changes: 4 additions & 4 deletions python/lsst/atmospec/processStar.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.")

Expand Down
26 changes: 15 additions & 11 deletions python/lsst/atmospec/spectraction.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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):
Expand All @@ -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,

Check failure on line 296 in python/lsst/atmospec/spectraction.py

View workflow job for this annotation

GitHub Actions / call-workflow / lint

F821

undefined name 'butler'

Check failure on line 296 in python/lsst/atmospec/spectraction.py

View workflow job for this annotation

GitHub Actions / call-workflow / lint

F821

undefined name 'dataId'
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
Expand All @@ -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
Expand Down
73 changes: 49 additions & 24 deletions python/lsst/atmospec/spectroFlat.py
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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=}.')
Expand All @@ -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,
Expand Down
1 change: 1 addition & 0 deletions python/lsst/atmospec/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down

0 comments on commit 99b720e

Please sign in to comment.