Spectro flat and mask to improve spectral extraction

config/auxtel.ini

@@ -19,6 +19,10 @@ SPECTRACTOR_DECONVOLUTION_SIGMA_CLIP = 100
 SPECTRACTOR_FIT_TIMEOUT_PER_ITER = 1200
 # maximum time per gradient descent before TimeoutError in seconds
 SPECTRACTOR_FIT_TIMEOUT = 7200
+# library to compute atmospheric transmission: none, libradtran, getobsatmo
+SPECTRACTOR_ATMOSPHERE_SIM = getobsatmo
+# simulate star field with Gaia catalog: False, True
+SPECTRACTOR_SIMULATE_STARFIELD = False
 
 [instrument]
 # instrument name
@@ -53,7 +57,7 @@ CCD_PIXEL2ARCSEC = 0.0952
 # approximate maximum ADU output of the CCD
 CCD_MAXADU = 170000
 # electronic gain : elec/ADU
-CCD_GAIN = 1.1
+CCD_GAIN = 1.3
 # rebinning of the image in pixel
 CCD_REBIN = 2
 

python/lsst/atmospec/processStar.py

@@ -304,6 +304,8 @@
             -1: "The first order spectrum in the negative y direction",
             2: "The second order spectrum in the positive y direction",
             -2: "The second order spectrum in the negative y direction",
+            3: "The third order spectrum in the positive y direction",
+            -3: "The third order spectrum in the negative y direction",
         }
     )
     signalWidth = pexConfig.Field(  # TODO: change this to be set wrt the focus/seeing, i.e. FWHM from imChar
@@ -509,11 +511,17 @@
     def validate(self):
         super().validate()
         uvspecPath = shutil.which('uvspec')
-        if uvspecPath is None and self.doFitAtmosphere is True:
-            raise FieldValidationError(self.__class__.doFitAtmosphere, self, "uvspec is not in the path,"
+        try:
+            import getObsAtmo
+        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,"
                                        " but doFitAtmosphere is True.")
-        if uvspecPath is None and self.doFitAtmosphereOnSpectrogram is True:
-            raise FieldValidationError(self.__class__.doFitAtmosphereOnSpectrogram, self, "uvspec is not in"
+        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.")
 
 

python/lsst/atmospec/spectraction.py

@@ -44,6 +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
 
 
 class SpectractorShim:
@@ -213,11 +214,18 @@
         if parameters.DEBUG:
             self.debugPrintTargetCentroidValue(image)
 
-        self._setReadNoiseFromExp(image, exp, 1)
+        self._setReadNoiseFromExp(image, exp, 8.5)
         # xxx remove hard coding of 1 below!
-        image.gain = self._setGainFromExp(image, exp, .85)  # gain required for calculating stat err
+        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'])
+        PTCGainDict = getPTCGainDict(butler)
+        certifiedFlat = getCertifiedFlat(butler, dataId=exp.visitInfo.id, filter="empty")
+        image.gain = self._transformArrayFromExpToImage(makeGainFlat(certifiedFlat, PTCGainDict).image.array)
         self._setStatErrorInImage(image, exp, useExpVariance=False)
-        # image.coord as an astropy SkyCoord - currently unused
+        self._setMask(image, exp)
+        sensorFlat = makeSensorFlat(certifiedFlat, PTCGainDict, kernel="mean", invertGains=True)
+        image.flat = self._transformArrayFromExpToImage(sensorFlat.image.array)
 
         self._setImageAndHeaderInfo(image, exp)  # sets image attributes
 
@@ -271,7 +279,11 @@
             data = exp.image.array[0:4000, 0:4000]
         else:
             data = exp.image.array
-        return data.T[::, ::]
+        return self._transformArrayFromExpToImage(data)
+
+    def _transformArrayFromExpToImage(self, array):
+        # apply transformation on an exposure array to have correct orientation in Spectractor Image
+        return array.T[::, ::]
 
     def _setReadNoiseFromExp(self, spectractorImage, exp, constValue=None):
         # xxx need to implement this properly
@@ -281,16 +293,32 @@
             # 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_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()
+                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,')')
+
 
     def _setReadNoiseToNone(self, spectractorImage):
         spectractorImage.read_out_noise = None
 
     def _setStatErrorInImage(self, image, exp, useExpVariance=False):
         if useExpVariance:
-            image.stat_errors = exp.maskedImage.variance.array  # xxx need to deal with TRANSPOSE here
+            image.err = self._transformArrayFromExpToImage(np.sqrt(exp.getVariance().array))
         else:
             image.compute_statistical_error()
 
+    def _setMask(self, image, exp):
+        BAD = exp.getMask().getPlaneBitMask("BAD")
+        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
@@ -448,9 +476,9 @@
                                                               signal_width=parameters.PIXWIDTH_SIGNAL,
                                                               ws=(parameters.PIXDIST_BACKGROUND,
                                                                   parameters.PIXDIST_BACKGROUND
-                                                                  + parameters.PIXWIDTH_BACKGROUND),
-                                                              right_edge=image.data.shape[1])
+                                                                  + parameters.PIXWIDTH_BACKGROUND))
         spectrum.atmospheric_lines = atmospheric_lines
+        spectrum.plot_spectrum()
 
         # PSF2D deconvolution
         if parameters.SPECTRACTOR_DECONVOLUTION_PSF2D:
@@ -467,8 +495,8 @@
 
         # not necessarily set during fit but required to be present for astropy
         # fits writing to work (required to be in keeping with upstream)
-        spectrum.data_order2 = np.zeros_like(spectrum.lambdas)
-        spectrum.err_order2 = np.zeros_like(spectrum.lambdas)
+        spectrum.data_next_order = np.zeros_like(spectrum.lambdas)
+        spectrum.err_next_order = np.zeros_like(spectrum.lambdas)
 
         # Full forward model extraction:
         # adds transverse ADR and order 2 subtraction

python/lsst/atmospec/spectroFlat.py

@@ -0,0 +1,201 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.colors import LogNorm
+import logging
+
+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}'"
+    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=}')
+
+    flat_dates = np.sort(np.unique([flat_.day_obs for flat_ in flat_records]))
+    ids_ = [flat_.id for flat_ in flat_records]
+
+    flat_dict_ids = {}
+    for date_ in flat_dates:
+        flat_dict_ids[date_] = []
+        for id_ in ids_:
+            if str(date_) in str(id_):
+                flat_dict_ids[date_].append(id_)
+        flat_dict_ids[date_] = np.sort(flat_dict_ids[date_])
+
+    return flat_dates, flat_dict_ids
+
+
+def getCertifiedFlat(butler, dataId, filter='empty', disperser='empty'):
+    flat_dates, flat_ids = find_flats(butler, filter=filter, disperser=disperser)
+    obs_day = dataId // 100_000
+    date_diff = int(obs_day) - flat_dates
+    closest_date = flat_dates[np.argmax(date_diff[date_diff <= 0])]
+    # load flat
+    flat_id = flat_ids[closest_date][-1]
+    certifiedFlat = butler.get('flat', instrument='LATISS', exposure=flat_id, detector=0,
+                               collections=['LATISS/calib', 'LATISS/raw/all'])
+    return certifiedFlat
+
+
+def getGainDict(exposure):
+    det = exposure.getDetector()
+    gainDict = {}
+    for amp in det:
+        gainDict[amp.getName()] = amp.getGain()
+    return gainDict
+
+
+def getPTCGainDict(butler):
+    ptc = butler.get('ptc', instrument="LATISS", detector=0, collections='u/cslage/sdf/latiss/ptc_20220927J')
+    PTCGainDict = ptc.gain
+    return PTCGainDict
+
+
+def makeGainFlat(certifiedFlat, gainDict, invertGains=False):
+    """Given an exposure, make a flat from the gains.
+
+    Construct an exposure where the image array data contains the gain of the
+    pixel, such that dividing by (or mutiplying by) the image will convert
+    an image from ADU to electrons.
+
+    Parameters
+    ----------
+    certifiedFlat : `lsst.afw.image.exposure`
+        The certifiedFlat from which the flat is to be made.
+
+    gainDict : `dict` of `float`
+        A dict of the amplifiers' gains, keyed by the amp names.
+
+    invertGains : `bool`
+        Gains are specified in inverted units and should be applied as such.
+
+    Returns
+    -------
+    gainFlat : `lsst.afw.image.exposure`
+        The gain flat
+    """
+    flat = certifiedFlat.clone()
+    detector = flat.getDetector()
+    ampNames = set(list(a.getName() for a in detector))
+    assert set(gainDict.keys()) == ampNames
+
+    gainDict_norm = {}
+    mean = np.mean(list(gainDict.values()))
+    for amp in gainDict.keys():
+        gainDict_norm[amp] = gainDict[amp] / mean
+
+    for amp in detector:
+        bbox = amp.getBBox()
+        if invertGains:
+            flat[bbox].maskedImage.image.array[:, :] = 1. / gainDict_norm[amp.getName()]
+        else:
+            flat[bbox].maskedImage.image.array[:, :] = gainDict_norm[amp.getName()]
+    flat.maskedImage.mask[:] = 0x0
+    flat.maskedImage.variance[:] = 0.0
+
+    return flat
+
+
+def makeOpticFlat(certifiedFlat, kernel='mean', window_size=40, mode='mirror', percentile=1.):
+    logger = logging.getLogger()
+    opticFlat = certifiedFlat.clone()
+
+    detector = opticFlat.getDetector()
+    for amp in detector:
+        bbox = amp.getBBox()
+        data = np.copy(opticFlat[bbox].maskedImage.image.array[:, :])
+        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)
+        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))
+            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)
+        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),
+                                                                          mode=mode)
+        else:
+            raise ValueError(f'Kernel must be within ["mean", "median", "gauss"]. Got {kernel=}.')
+
+    opticFlat.maskedImage.mask[:] = 0x0
+    opticFlat.maskedImage.variance[:] = 0.0
+    return opticFlat
+
+
+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)
+
+    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[:, :]
+        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
+
+
+
+
+    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):
+    logger = logging.getLogger()
+    logger.info(f'Window size for {kernel} smoothing = {window_size}')
+    pixelFlat = makePixelFlat(certifiedFlat, kernel=kernel, window_size=window_size, mode=mode,
+                              percentile=percentile, remove_background=remove_background)
+    gainFlat = makeGainFlat(certifiedFlat, gainDict, invertGains=invertGains)
+    sensorFlat = pixelFlat.clone()
+
+    detector = sensorFlat.getDetector()
+    ampNames = set(list(a.getName() for a in detector))
+    assert set(gainDict.keys()) == ampNames
+
+    for amp in detector:
+        bbox = amp.getBBox()
+        sensorFlat[bbox].maskedImage.image.array[:, :] *= gainFlat[bbox].maskedImage.image.array[:, :]
+
+    sensorFlat.maskedImage.mask[:] = 0x0
+    sensorFlat.maskedImage.variance[:] = 0.0
+    return sensorFlat
+
+
+def plotFlat(flat, title=None, figsize=(10, 10), cmap='gray', vmin=0.9, vmax=1.1, lognorm=False):
+    fig = plt.figure(figsize=figsize)
+    ax = fig.add_subplot(111)
+    if lognorm:
+        im = ax.imshow(flat.image.array, cmap=cmap, origin='lower', norm=LogNorm())
+    else:
+        im = ax.imshow(flat.image.array, cmap=cmap, origin='lower', vmin=vmin, vmax=vmax)
+    if title is not None:
+        ax.set_title(title)
+    plt.colorbar(im, ax=ax)
+    plt.tight_layout()
+    return

python/lsst/atmospec/utils.py

@@ -59,6 +59,13 @@
 from astropy.coordinates import SkyCoord, Distance
 
 
+def getGainDict(exposure):
+    det = exposure.getDetector()
+    gainDict = {}
+    for amp in det:
+        gainDict[amp.getName()] = amp.getGain()
+    return gainDict
+
 def makeGainFlat(exposure, gainDict, invertGains=False):
     """Given an exposure, make a flat from the gains.