Calibration: Steffl method

Recreating Steffl’s Flatfields

Based on his thesis appendix

⚠️ Demo cells below are marked #| eval: false — they require live PDS index fetches, cached data products, or interactive plotting backends that are unavailable during a clean docs render. Open this notebook in JupyterLab to step through it cell-by-cell against your local data.

from functools import cached_property

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr
from tqdm.auto import tqdm, trange

import holoviews as hv
import hvplot.xarray  # noqa: F401
import panel as pn
from panel import widgets

from planetarypy.datetime_format_converters import (
    fromdoyformat,
    doyformat,  # replaces removed iso_to_nasa_date
)
from pyuvis import UVPDS, UVISObs, CatalogFilter
from pyuvis.calib.greg import filter_spica_for_date
from pyuvis.calib.steffl import (
    steffl_spica_dates,
    steffl_spica_doy_dates,
    Row2Row,
    create_detector_stack,
)

# `param` was used for class definitions in the original; the
# implementation now lives in src/pyuvis/calib/steffl.py — only
# import param here if you're prototyping new classes interactively.
# import param

# `find_mixed_type_cols` is no longer exposed by planetarypy.pds.indexes;
# any cell below that referenced it is left as #| eval: false.

hv.extension("bokeh", logo=False)

missing = [] there = [] for id in obsids: try: data = UVPDS(id, skip_download=True) except FileNotFoundError: print(id, “not there.”) missing.append(id) else: print(“Got”, id) there.append(id)

cat = CatalogFilter(steffl_spica_dates[2])

pids = list(cat.get_euv_date().query("OBSERVATION_TYPE=='CALIB'").index)

pids

cat.set_next_day()

pids.extend(list(cat.get_euv_date().query("OBSERVATION_TYPE=='CALIB'").index))

pids

kwargs = {"x": "nx", "y": "ny", "cmap": "viridis", "clim": (0, 50)}

Row2Row correction

pids

r2r = Row2Row(pids[0])

r2r.plot_set

r2r.plot_averaged

r2r.plot_ff

r2r.plot_integrated

r2r.plot_column_std

r2r.ff

Col2Col sensitivity variation.

Universal detector stack xarray creator

This function allows to create detector-shaped stacks of data for different purposes, like: * along_slit scans * across_slit scans * sets of corrections for the Steffl flatfield

c2c = Col2Col(pids)

c2c.calculate_simple_correction()

c2c.plot_all_Fm()

• It is neat to see how the grid has a weaker signal in the last array, where only 2 columns could be combined

c2c.plot_simple_correction()

Row2Row flats

Every across_slit scan comes with an inner along_slit scan that can be used to create a row2row flatfield correction.

So we have 14 of them in this case, that can be averaged.

c2c.plot_r2r_flats()

c2c.r2r_flat_average_plot()

• Nice to see how the SNR below 155 nm has improved drastically from having the average of the 14 FF corrections.

c2c.simple_both_plot()

c2c.i=100
c2c.m=0

c2c.plot_triplet()

c2c.plot_triplet(corrected=True)

def compare_before_after(i=0):
    p1 = (
        c2c.arr.where(c2c.arr > 0, np.nan)
        .isel(across_slit=i, drop=True)
        .hvplot(
            x="spectral",
            y="spatial",
            cmap="viridis",
            logz=True,
            label="Original",
            # widget_type="scrubber",
            widget_location="bottom",
        )
    )
    p2 = (
        c2c.corrected_arr.where(c2c.corrected_arr > 0, np.nan)
        .isel(across_slit=i, drop=True)
        .hvplot.image(
            x="spectral",
            y="spatial",
            cmap="viridis",
            logz=True,
            label="Corrected",
            # widget_type="scrubber",
            widget_location="bottom",
        )
    )
    return pn.Column(p1) + pn.Column(p2)

slider = pn.widgets.IntSlider(end=14)

pn.interact(compare_before_after, i=(0, 14))

p2 = c2c.corrected_arr.where(c2c.corrected_arr > 0, np.nan).isel(across_slit=i, drop=True).hvplot.image(
    x="spectral",
    y="spatial",
    cmap="viridis",
    logz=True,
    label="Corrected",
    # widget_type="scrubber",
    widget_location="bottom",
)

import panel as pn

pn.Column(p1) + pn.Column(p2)

(p1+p2)

Mean of inner area is 1:

c2c.simple_both.isel(spectral=slice(15, 998), spatial=slice(3,61)).mean().data

iarr = c2c.arr.interactive()

### import panel as pn
from panel import widgets

ival = widgets.IntSlider(start=0, end=1024)
mval = widgets.IntSlider()

iarr.isel(spectral=ival, drop=True).hvplot().layout()

c2c.calculate_simple_correction()

c2c.simple_correction

c2c.i = 200
c2c.m = 0

c2c.column_set_mean() / c2c.column_set()[0]

c2c.column_set()[0]

c2c.i = 200

c2c.m = 3

c2c.plot_averaged_triplet()

c2c.plot_triplet(corrected=False, i=300)

c2c.plot_simple_correction()

c2c.simple_correction.shape

c2c.plot_Fm(0)

c2c.plot_triplet()

c2c.column_set_mean()

c2c.plot()

archive_df.loc["EUV2001_093_08_35_28"]

p = obsdir / "index_repaired.tab"

df = pd.read_csv(p, quotechar='"', skipinitialspace=True)

from planetarypy.pds.indexes import find_mixed_type_cols

find_mixed_type_cols(df, fix=False)

df.columns

index.columns

index[index.filename.str.startswith("EUV")].iloc[0]

obs.head()

cols = ["index start_time stop_time detector target obsid_time unknown type comment "]