Image Simulation

This notebook demonstrates how to generate synthetic ESIS images using data from the Interface Region Imaging Spectrograph (IRIS) [De Pontieu et al., 2014].

import IPython.display
import numpy as np
import matplotlib.pyplot as plt
import named_arrays as na
import astropy.units as u
import astropy.time
import astropy.visualization
import optika
import iris
import esis

astropy.visualization.quantity_support();

Download a sequence of 320-step dense rasters, despike, and crop to \(\pm 100 \text{km/s}\).

obs = iris.sg.SpectrographObservation.from_time_range(
    time_start=astropy.time.Time("2021-09-23T14:00"),
    time_stop=astropy.time.Time("2021-09-23T14:01"),
)

obs = obs.explicit

bg = iris.sg.background.estimate(
    obs=obs,
    axis_time=obs.axis_time,
    axis_wavelength=obs.axis_wavelength,
    axis_detector_x=obs.axis_detector_x,
    axis_detector_y=obs.axis_detector_y,
)

obs = obs - bg.outputs

obs = na.despike(obs, axis=(obs.axis_wavelength, obs.axis_detector_y))

velocity_thresh = 150 * u.km / u.s

velocity_centers = obs.velocity_doppler.cell_centers().mean(obs.axis_time)

index_lower = np.argmax(-velocity_thresh < velocity_centers)[obs.axis_wavelength]
index_upper = np.argmax(velocity_thresh < velocity_centers)[obs.axis_wavelength]

crop_wavelength = {obs.axis_wavelength: slice(index_lower.ndarray, index_upper.ndarray)}

obs = obs[crop_wavelength]

obs.shape

/opt/hostedtoolcache/Python/3.11.13/x64/lib/python3.11/site-packages/numpy/lib/nanfunctions.py:1217: RuntimeWarning: All-NaN slice encountered
  return function_base._ureduce(a, func=_nanmedian, keepdims=keepdims,

{'time': 3, 'detector_x': 320, 'detector_y': 548, 'wavelength': 55}

obs.to_jshtml(
    velocity_min=obs.velocity_doppler.ndarray.min(),
    velocity_max=obs.velocity_doppler.ndarray.max(),
    fps=2,
)

Once Loop Reflect

fwhm_si_iv = (obs.wavelength_center + 0.188 * u.AA).to(u.km/u.s, equivalencies=u.doppler_radio(obs.wavelength_center))
fwhm_si_iv

$40.432048 \; \mathrm{\frac{km}{s}}$

wavelength_o_v = 629.73 * u.AA
fwhm_o_v = (wavelength_o_v + 0.129 * u.AA).to(u.km/u.s, equivalencies=u.doppler_radio(wavelength_o_v))
fwhm_o_v

$61.399817 \; \mathrm{\frac{km}{s}}$

width_ratio = fwhm_o_v / fwhm_si_iv
width_ratio

$1.5185928 \; \mathrm{}$

velocity = width_ratio * obs.velocity_doppler
obs.inputs.wavelength = velocity.to(u.AA, equivalencies=u.doppler_radio(wavelength_o_v))
obs.wavelength_center = wavelength_o_v

veranazza_radiance = 690.80 * u.erg / u.cm ** 2 / u.sr / u.s

dw = np.diff(obs.inputs.wavelength, axis=obs.axis_wavelength)
avg = np.nanmean((dw * obs.outputs).sum(obs.axis_wavelength)[{obs.axis_time: 0}]).ndarray
avg

$5.0010132 \; \mathrm{\mathring{A}\,DN}$

obs.outputs.shape

{'time': 3, 'detector_x': 320, 'detector_y': 548, 'wavelength': 55}

obs.inputs.wavelength.shape

{'time': 3, 'wavelength': 56}

obs = obs / avg * veranazza_radiance
obs = obs.to(u.erg / u.cm ** 2 / u.sr / u.s / u.AA)

obs.outputs = np.nan_to_num(obs.outputs)

np.nanmean((dw * obs.outputs).sum(obs.axis_wavelength)).ndarray

$1265.1031 \; \mathrm{\frac{erg}{s\,sr\,cm^{2}}}$

obs.inputs.position.x -= obs.inputs.position.x.mean((obs.axis_detector_x, obs.axis_detector_y))
obs.inputs.position.y -= obs.inputs.position.y.mean((obs.axis_detector_x, obs.axis_detector_y))

instrument = esis.flights.f1.optics.design_single(num_distribution=0)
instrument.primary_mirror.translation = na.nominal(instrument.primary_mirror.translation)
instrument.field_stop.translation = na.nominal(instrument.field_stop.translation)
instrument.grating.sag.radius = na.nominal(instrument.grating.sag.radius)
instrument.grating.rulings.spacing.coefficients[0] = na.nominal(instrument.grating.rulings.spacing.coefficients[0])
instrument.grating.rulings.spacing.coefficients[1] = na.nominal(instrument.grating.rulings.spacing.coefficients[1])
instrument.grating.rulings.spacing.coefficients[2] = na.nominal(instrument.grating.rulings.spacing.coefficients[2])
instrument.grating.rulings.depth = na.nominal(instrument.grating.rulings.depth)
instrument.grating.rulings.ratio_duty = na.nominal(instrument.grating.rulings.ratio_duty)
instrument.grating.translation = na.nominal(instrument.grating.translation)
instrument.grating.roll = na.nominal(instrument.grating.roll)

system = instrument.system
# system.sensor.material = optika.sensors.IdealImagingSensorMaterial()

pupil = na.Cartesian2dVectorLinearSpace(
    start=-1,
    stop=1,
    axis=na.Cartesian2dVectorArray("pupil_x", "pupil_y"),
    num=2,
)

%%time
image = system.image(
    obs,
    pupil=pupil,
    axis_wavelength=obs.axis_wavelength,
    axis_field=(obs.axis_detector_x, obs.axis_detector_y),
)

WARNING: function 'sqrt' is not known to astropy's Quantity. Will run it anyway, hoping it will treat ndarray subclasses correctly. Please raise an issue at https://github.com/astropy/astropy/issues. [astropy.units.quantity]
/opt/hostedtoolcache/Python/3.11.13/x64/lib/python3.11/site-packages/astropy/units/quantity.py:653: RuntimeWarning: divide by zero encountered in divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/opt/hostedtoolcache/Python/3.11.13/x64/lib/python3.11/site-packages/astropy/units/quantity.py:653: RuntimeWarning: overflow encountered in exp
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/opt/hostedtoolcache/Python/3.11.13/x64/lib/python3.11/site-packages/astropy/units/quantity.py:653: RuntimeWarning: invalid value encountered in multiply
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/opt/hostedtoolcache/Python/3.11.13/x64/lib/python3.11/site-packages/astropy/units/quantity.py:653: RuntimeWarning: divide by zero encountered in divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/opt/hostedtoolcache/Python/3.11.13/x64/lib/python3.11/site-packages/named_arrays/_scalars/scalars.py:552: RuntimeWarning: invalid value encountered in divide
  result_ndarray = getattr(function, method)(*inputs_ndarray, **kwargs_ndarray)

CPU times: user 10min 31s, sys: 8min 3s, total: 18min 34s
Wall time: 21min 19s

r = obs.copy()
factor = 10000
r.outputs = factor * np.nan_to_num(r.outputs)
r.inputs.wavelength = (r.velocity_doppler/factor).to(u.AA, equivalencies=u.doppler_radio(wavelength_o_v))

%%time
image_undispersed = system.image(
    r,
    pupil=pupil,
    axis_wavelength=obs.axis_wavelength,
    axis_field=(obs.axis_detector_x, obs.axis_detector_y),
    noise=False,
)

WARNING: function 'sqrt' is not known to astropy's Quantity. Will run it anyway, hoping it will treat ndarray subclasses correctly. Please raise an issue at https://github.com/astropy/astropy/issues. [astropy.units.quantity]
/opt/hostedtoolcache/Python/3.11.13/x64/lib/python3.11/site-packages/astropy/units/quantity.py:653: RuntimeWarning: divide by zero encountered in divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)

CPU times: user 10min 2s, sys: 8min 10s, total: 18min 12s
Wall time: 20min 49s

crop = dict(
    detector_x=slice(1024 + 256, 1024 + 512 + 256),
    detector_y=slice(256, 512 + 256),
)
image = image[crop]
image_undispersed = image_undispersed[crop]

stacked = na.stack(
    arrays=[
        image_undispersed,
        image,
    ],
    axis="dispersion",
)

fig, ax = na.plt.subplots(
    axis_cols="dispersion",
    ncols=3,
    # sharex=True,
    # sharey=True,
    figsize=(9, 4.3),
    constrained_layout=True,
    dpi=200,
    gridspec_kw=dict(width_ratios=[.9, .9, .1]),
)
cax = ax.ndarray[~0]
cax_twin = cax.twinx()
ax = ax[dict(dispersion=slice(0, 2))]
cax.xaxis.set_ticks_position("top")
cax.xaxis.set_label_position("top")
ani, colorbar = na.plt.rgbmovie(
    obs.inputs.time,
    image.inputs.wavelength,
    image.inputs.position.x.to(u.mm),
    image.inputs.position.y.to(u.mm),
    C=stacked.outputs,
    axis_time=obs.axis_time,
    axis_wavelength=obs.axis_wavelength,
    ax=ax,
    vmin=0,
    vmax=np.nanpercentile(image.outputs, 99.9, axis=(obs.axis_time, "detector_x", "detector_y")),
)
colorbar = colorbar[{obs.axis_time: 0}]
na.plt.pcolormesh(
    colorbar.inputs.x,
    colorbar.inputs.y,
    C=colorbar.outputs.value,
    axis_rgb=obs.axis_wavelength,
    ax=cax,
)
na.plt.pcolormesh(
    colorbar.inputs.x,
    colorbar.inputs.y.to(u.km / u.s, equivalencies=u.doppler_radio(wavelength_o_v)),
    C=colorbar.outputs.value,
    axis_rgb=obs.axis_wavelength,
    ax=cax_twin,
)
na.plt.set_xlim(5.5, 9, ax=ax)
na.plt.set_ylim(-2, 2, ax=ax)
na.plt.set_aspect("equal", ax=ax)
na.plt.set_xlabel("detector $x$ (mm)", ax=ax)
na.plt.set_ylabel("detector $y$ (mm)", ax=ax.ndarray[0])
na.plt.set_ylabel("", ax=ax.ndarray[~0])
ax.ndarray[~0].tick_params(axis="y", labelleft=False)

result = ani.to_jshtml(fps=2)
result = IPython.display.HTML(result)

plt.close(ani._fig)

result

Once Loop Reflect

t0 = {obs.axis_time: ~0}
s = stacked[t0]

fig, ax = plt.subplots(
    ncols=2,
    constrained_layout=True,
    dpi=200,
    gridspec_kw=dict(width_ratios=[.9, .1]),
    figsize=(7, 5),
)
ax_twin = ax[1].twinx()
ax[1].xaxis.set_ticks_position("top")
ax[1].xaxis.set_label_position("top")
ani, colorbar = na.plt.rgbmovie(
    obs.inputs.time[t0],
    image.inputs.wavelength[t0],
    s.inputs.position.x.to(u.mm),
    s.inputs.position.y.to(u.mm),
    C=na.nominal(s.outputs),
    axis_time="dispersion",
    axis_wavelength=obs.axis_wavelength,
    ax=ax[0],
    vmin=0,
    vmax=np.nanpercentile(image.outputs, 99.9, axis=(obs.axis_time, "detector_x", "detector_y")),
)
# colorbar = colorbar[{obs.axis_time: 0}]
na.plt.pcolormesh(
    colorbar.inputs.x,
    colorbar.inputs.y,
    C=colorbar.outputs.value,
    axis_rgb=obs.axis_wavelength,
    ax=ax[1],
)
na.plt.pcolormesh(
    colorbar.inputs.x,
    colorbar.inputs.y.to(u.km / u.s, equivalencies=u.doppler_radio(wavelength_o_v)),
    C=colorbar.outputs.value,
    axis_rgb=obs.axis_wavelength,
    ax=ax_twin,
)
na.plt.set_xlim(5.5, 9, ax=ax[0])
na.plt.set_ylim(-2, 2, ax=ax[0])
na.plt.set_aspect("equal", ax=ax[0])
na.plt.set_xlabel(f"detector $x$ (mm)", ax=ax[0])
na.plt.set_ylabel("detector $y$ (mm)", ax=ax[0])

result = ani.to_jshtml(fps=.5)
result = IPython.display.HTML(result)

plt.close(ani._fig)

result

Once Loop Reflect