Note
Go to the end to download the full example code.
Imaging Demo#
How to create visibility from pixel data and make images.
The example uses stixpy to obtain STIX pixel data and convert these into visibilities and xrayvisim
to make the images.
Imports
import logging
import matplotlib.pyplot as plt
import numpy as np
from xrayvision.clean import vis_clean
from xrayvision.imaging import vis_to_image, vis_to_map
from xrayvision.mem import mem, resistant_mean
import astropy.units as u
from astropy.coordinates import SkyCoord
from sunpy.coordinates import HeliographicStonyhurst, Helioprojective
from sunpy.map import Map, make_fitswcs_header
from sunpy.time import TimeRange
from stixpy.calibration.visibility import calibrate_visibility, create_meta_pixels, create_visibility
from stixpy.coordinates.frames import STIXImaging
from stixpy.coordinates.transforms import get_hpc_info
from stixpy.imaging.em import em
from stixpy.map.stix import STIXMap # noqa
from stixpy.product import Product
logger = logging.getLogger(__name__)
Read science file as Product
cpd_sci = Product(
"http://pub099.cs.technik.fhnw.ch/fits/L1/2021/09/23/SCI/solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits"
)
cpd_sci
Files Downloaded: 0%| | 0/1 [00:00<?, ?file/s]
solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits: 0%| | 0.00/48.4M [00:00<?, ?B/s]
solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits: 0%| | 1.02k/48.4M [00:00<5:53:04, 2.28kB/s]
solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits: 0%| | 8.94k/48.4M [00:00<39:29, 20.4kB/s]
solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits: 0%| | 64.6k/48.4M [00:00<05:17, 152kB/s]
solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits: 0%| | 225k/48.4M [00:00<01:32, 518kB/s]
solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits: 1%| | 569k/48.4M [00:00<00:37, 1.27MB/s]
solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits: 3%|▎ | 1.23M/48.4M [00:00<00:17, 2.66MB/s]
solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits: 5%|▌ | 2.55M/48.4M [00:01<00:08, 5.44MB/s]
solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits: 11%|█ | 5.17M/48.4M [00:01<00:03, 10.9MB/s]
solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits: 22%|██▏ | 10.4M/48.4M [00:01<00:01, 22.0MB/s]
solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits: 40%|████ | 19.5M/48.4M [00:01<00:00, 40.7MB/s]
solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits: 58%|█████▊ | 27.9M/48.4M [00:01<00:00, 52.9MB/s]
solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits: 77%|███████▋ | 37.5M/48.4M [00:01<00:00, 65.0MB/s]
solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits: 96%|█████████▌| 46.4M/48.4M [00:01<00:00, 71.9MB/s]
Files Downloaded: 100%|██████████| 1/1 [00:02<00:00, 2.55s/file]
Files Downloaded: 100%|██████████| 1/1 [00:02<00:00, 2.55s/file]
CompressedPixelData <sunpy.time.timerange.TimeRange object at 0x7e2f981321e0>
Start: 2021-09-23 15:20:15
End: 2021-09-23 15:26:39
Center:2021-09-23 15:23:27
Duration:0.004439814814814813 days or
0.10655555555555551 hours or
6.393333333333331 minutes or
383.59999999999985 seconds
DetectorMasks
[0...697]: [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31]
PixelMasks
[0...697]: [['1' '1' '1' '1' '1' '1' '1' '1' '1' '1' '1' '1']]
EnergyEdgeMasks
[0]: [_,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,_]
Read background file as Product
cpd_bkg = Product(
"http://pub099.cs.technik.fhnw.ch/fits/L1/2021/09/23/SCI/solo_L1_stix-sci-xray-cpd_20210923T095923-20210923T113523_V02_2109230083-57078.fits"
)
cpd_bkg
Files Downloaded: 0%| | 0/1 [00:00<?, ?file/s]
solo_L1_stix-sci-xray-cpd_20210923T095923-20210923T113523_V02_2109230083-57078.fits: 0%| | 0.00/121k [00:00<?, ?B/s]
solo_L1_stix-sci-xray-cpd_20210923T095923-20210923T113523_V02_2109230083-57078.fits: 1%| | 1.02k/121k [00:00<00:53, 2.26kB/s]
solo_L1_stix-sci-xray-cpd_20210923T095923-20210923T113523_V02_2109230083-57078.fits: 7%|▋ | 8.95k/121k [00:00<00:05, 20.4kB/s]
solo_L1_stix-sci-xray-cpd_20210923T095923-20210923T113523_V02_2109230083-57078.fits: 47%|████▋ | 56.9k/121k [00:00<00:00, 131kB/s]
Files Downloaded: 100%|██████████| 1/1 [00:01<00:00, 1.49s/file]
Files Downloaded: 100%|██████████| 1/1 [00:01<00:00, 1.49s/file]
CompressedPixelData <sunpy.time.timerange.TimeRange object at 0x7e2f898fb770>
Start: 2021-09-23 09:59:23
End: 2021-09-23 11:35:23
Center:2021-09-23 10:47:23
Duration:0.06666666666666665 days or
1.5999999999999996 hours or
95.99999999999997 minutes or
5759.999999999999 seconds
DetectorMasks
[0]: [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31]
PixelMasks
[0]: [['1' '1' '1' '1' '1' '1' '1' '1' '1' '1' '1' '1']]
EnergyEdgeMasks
[0]: [_,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32]
Set time and energy ranges which will be considered for the science and the background file
time_range_sci = ["2021-09-23T15:20:00", "2021-09-23T15:23:00"]
time_range_bkg = [
"2021-09-23T09:00:00",
"2021-09-23T12:00:00",
] # Set this range larger than the actual observation time
energy_range = [25, 28] * u.keV
Create the meta pixel, A, B, C, D for the science and the background data
meta_pixels_sci = create_meta_pixels(
cpd_sci, time_range=time_range_sci, energy_range=energy_range, flare_location=[0, 0] * u.arcsec, no_shadowing=True
)
meta_pixels_bkg = create_meta_pixels(
cpd_bkg, time_range=time_range_bkg, energy_range=energy_range, flare_location=[0, 0] * u.arcsec, no_shadowing=True
)
Perform background subtraction
meta_pixels_bkg_subtracted = {
**meta_pixels_sci,
"abcd_rate_kev_cm": meta_pixels_sci["abcd_rate_kev_cm"] - meta_pixels_bkg["abcd_rate_kev_cm"],
"abcd_rate_error_kev_cm": np.sqrt(
meta_pixels_sci["abcd_rate_error_kev_cm"] ** 2 + meta_pixels_bkg["abcd_rate_error_kev_cm"] ** 2
),
}
Create visibilities from the meta pixels
vis = create_visibility(meta_pixels_bkg_subtracted)
Obtain the necessary ephemeris data
vis_tr = TimeRange(vis.meta["time_range"])
roll, solo_xyz, pointing = get_hpc_info(vis_tr.start, vis_tr.end)
solo = HeliographicStonyhurst(*solo_xyz, obstime=vis_tr.center, representation_type="cartesian")
center_stix = SkyCoord(
0 * u.deg,
0 * u.deg,
frame=STIXImaging(obstime=vis_tr.start, obstime_end=vis_tr.end, observer=solo),
)
Files Downloaded: 0%| | 0/1 [00:00<?, ?file/s]
solo_ANC_stix-asp-ephemeris_20210923_V02.fits: 0%| | 0.00/397k [00:00<?, ?B/s]
solo_ANC_stix-asp-ephemeris_20210923_V02.fits: 0%| | 1.02k/397k [00:00<00:45, 8.68kB/s]
solo_ANC_stix-asp-ephemeris_20210923_V02.fits: 2%|▏ | 8.95k/397k [00:00<00:08, 43.4kB/s]
solo_ANC_stix-asp-ephemeris_20210923_V02.fits: 12%|█▏ | 48.9k/397k [00:00<00:01, 190kB/s]
solo_ANC_stix-asp-ephemeris_20210923_V02.fits: 28%|██▊ | 112k/397k [00:00<00:00, 349kB/s]
solo_ANC_stix-asp-ephemeris_20210923_V02.fits: 52%|█████▏ | 208k/397k [00:00<00:00, 541kB/s]
Files Downloaded: 100%|██████████| 1/1 [00:00<00:00, 1.07file/s]
Files Downloaded: 100%|██████████| 1/1 [00:00<00:00, 1.07file/s]
Calibrate the visibilities
If not given will default to center of STIX field-of-view
cal_vis = calibrate_visibility(vis, flare_location=center_stix)
Selected detectors 10 to 7
# order by sub-collimator e.g. 10a, 10b, 10c, 9a, 9b, 9c ....
isc_10_7 = [3, 20, 22, 16, 14, 32, 21, 26, 4, 24, 8, 28]
idx = np.argwhere(np.isin(cal_vis.meta["isc"], isc_10_7)).ravel()
Slice the visibilities to detectors 10 - 7
vis10_7 = cal_vis[idx]
Set up image parameters
imsize = [512, 512] * u.pixel # number of pixels of the map to reconstruct
pixel = [10, 10] * u.arcsec / u.pixel # pixel size in arcsec
Make a full disk back projection (inverse transform) map
bp_image = vis_to_image(vis10_7, imsize, pixel_size=pixel)
Obtain the necessary ephemeris data
header = make_fitswcs_header(
bp_image, center_stix, telescope="STIX", observatory="Solar Orbiter", scale=[10, 10] * u.arcsec / u.pix
)
fd_bp_map = Map((bp_image, header))
Convert the coordinates and make a map in Helioprojective and rotate so “North” is “up” Center of STIX pointing in HPC
center_hpc = center_stix.transform_to(Helioprojective(obstime=vis_tr.center))
header_hp = make_fitswcs_header(
bp_image, center_hpc, scale=[10, 10] * u.arcsec / u.pix, rotation_angle=90 * u.deg + roll
)
hp_map = Map((bp_image, header_hp))
hp_map_rotated = hp_map.rotate()
Plot the both maps
fig = plt.figure(layout="constrained", figsize=(12, 6))
ax = fig.subplot_mosaic(
[["stix", "hpc"]], per_subplot_kw={"stix": {"projection": fd_bp_map}, "hpc": {"projection": hp_map_rotated}}
)
fd_bp_map.plot(axes=ax["stix"])
fd_bp_map.draw_limb()
fd_bp_map.draw_grid()
hp_map_rotated.plot(axes=ax["hpc"])
hp_map_rotated.draw_limb()
hp_map_rotated.draw_grid()
<CoordinatesMap with 2 world coordinates:
index aliases type unit wrap format_unit visible
----- ------- --------- ---- --------- ----------- -------
0 lon longitude deg 180.0 deg deg yes
1 lat latitude deg None deg yes
>
Estimate the flare location and plot on top of back projection map. Note the coordinates are automatically converted from the STIXImaging to Helioprojective
max_pixel = np.argwhere(fd_bp_map.data == fd_bp_map.data.max()).ravel() * u.pixel
# because WCS axes and array are reversed
max_stix = fd_bp_map.pixel_to_world(max_pixel[1], max_pixel[0])
ax["stix"].plot_coord(max_stix, marker=".", markersize=50, fillstyle="none", color="r", markeredgewidth=2)
ax["hpc"].plot_coord(max_stix, marker=".", markersize=50, fillstyle="none", color="r", markeredgewidth=2)
fig.tight_layout()
/home/docs/checkouts/readthedocs.org/user_builds/stixpy/checkouts/latest/examples/imaging_demo.py:177: UserWarning: The figure layout has changed to tight
fig.tight_layout()
Use estimated flare location to create more accurate visibilities
meta_pixels_sci = create_meta_pixels(
cpd_sci, time_range=time_range_sci, energy_range=energy_range, flare_location=max_stix, no_shadowing=True
)
meta_pixels_bkg_subtracted = {
**meta_pixels_sci,
"abcd_rate_kev_cm": meta_pixels_sci["abcd_rate_kev_cm"] - meta_pixels_bkg["abcd_rate_kev_cm"],
"abcd_rate_error_kev_cm": np.sqrt(
meta_pixels_sci["abcd_rate_error_kev_cm"] ** 2 + meta_pixels_bkg["abcd_rate_error_kev_cm"] ** 2
),
}
vis = create_visibility(meta_pixels_bkg_subtracted)
cal_vis = calibrate_visibility(vis, flare_location=max_stix)
Selected detectors 10 to 3 order by sub-collimator e.g. 10a, 10b, 10c, 9a, 9b, 9c ….
isc_10_3 = [3, 20, 22, 16, 14, 32, 21, 26, 4, 24, 8, 28, 15, 27, 31, 6, 30, 2, 25, 5, 23, 7, 29, 1]
idx = np.argwhere(np.isin(cal_vis.meta["isc"], isc_10_3)).ravel()
Create an xrayvsion visibility object
cal_vis.meta["offset"] = max_stix
vis10_3 = cal_vis[idx]
Set up image parameters
imsize = [129, 129] * u.pixel # number of pixels of the map to reconstruct
pixel = [2, 2] * u.arcsec / u.pixel # pixel size in arcsec
Create a back projection image with natural weighting
bp_nat = vis_to_image(vis10_3, imsize, pixel_size=pixel)
Create a back projection image with uniform weighting
bp_uni = vis_to_image(vis10_3, imsize, pixel_size=pixel, scheme="uniform")
Create a sunpy.map.Map with back projection
bp_map = vis_to_map(vis10_3, imsize, pixel_size=pixel)
Crete a clean image using the clean algorithm vis_clean
niter = 200 # number of iterations
gain = 0.1 # gain used in each clean iteration
beam_width = 20.0 * u.arcsec
clean_map, model_map, resid_map = vis_clean(
vis10_3, imsize, pixel_size=pixel, gain=gain, niter=niter, clean_beam_width=20 * u.arcsec
)
2026-06-05T12:39:18Z INFO xrayvision.clean 124: Iter: 0, strength: 1.930919100604466, location: (np.int64(59), np.int64(67))
2026-06-05T12:39:18Z INFO xrayvision.clean 124: Iter: 25, strength: 0.47719851724349815, location: (np.int64(80), np.int64(54))
2026-06-05T12:39:18Z INFO xrayvision.clean 145: Largest residual negative
Create a sunpy map for the clean image in Helioprojective
header = make_fitswcs_header(
clean_map.data,
max_stix.transform_to(Helioprojective(obstime=vis_tr.center, observer=solo)),
telescope="STIX",
observatory="Solar Orbiter",
scale=pixel,
rotation_angle=90 * u.deg + roll,
)
Crete a map using the MEM GE algorithm mem
snr_value, _ = resistant_mean((np.abs(vis10_3.visibilities) / vis10_3.amplitude_uncertainty).flatten(), 3)
percent_lambda = 2 / (snr_value**2 + 90)
mem_map = mem(vis10_3, shape=imsize, pixel_size=pixel, percent_lambda=percent_lambda)
2026-06-05T12:39:20Z INFO xrayvision.mem 159: Iter: 0, Chi2: 234.49830788341825
2026-06-05T12:39:20Z INFO xrayvision.mem 159: Iter: 25, Chi2: 24.333702665147236
2026-06-05T12:39:20Z INFO xrayvision.mem 159: Iter: 50, Chi2: 19.266467368879347
2026-06-05T12:39:20Z INFO xrayvision.mem 159: Iter: 75, Chi2: 17.124105167153616
2026-06-05T12:39:20Z INFO xrayvision.mem 159: Iter: 100, Chi2: 15.880549155749382
2026-06-05T12:39:20Z INFO xrayvision.mem 159: Iter: 125, Chi2: 15.057967256302211
2026-06-05T12:39:20Z INFO xrayvision.mem 159: Iter: 150, Chi2: 14.507035644240354
2026-06-05T12:39:20Z INFO xrayvision.mem 159: Iter: 175, Chi2: 14.118968596360151
2026-06-05T12:39:20Z INFO xrayvision.mem 159: Iter: 200, Chi2: 13.825428062747665
2026-06-05T12:39:20Z INFO xrayvision.mem 159: Iter: 225, Chi2: 13.592951927443597
2026-06-05T12:39:20Z INFO xrayvision.mem 159: Iter: 250, Chi2: 13.399299686673773
2026-06-05T12:39:20Z INFO xrayvision.mem 510: Iter: 0, Obj function: (241.44033193177933+0j)
2026-06-05T12:39:23Z INFO xrayvision.mem 510: Iter: 25, Obj function: (29.805142705550765+0j)
2026-06-05T12:39:26Z INFO xrayvision.mem 510: Iter: 50, Obj function: (28.996740805172074+0j)
2026-06-05T12:39:30Z INFO xrayvision.mem 510: Iter: 75, Obj function: (28.920500700663403+0j)
2026-06-05T12:39:33Z INFO xrayvision.mem 510: Iter: 100, Obj function: (28.911357252278577+0j)
Crete a map using the EM algorithm EM
em_map = em(
meta_pixels_bkg_subtracted["abcd_rate_kev_cm"],
cal_vis,
shape=imsize,
pixel_size=pixel,
flare_location=max_stix,
idx=idx,
)
clean_map = Map((clean_map.data, header)).rotate()
bp_map = Map((bp_nat, header)).rotate()
mem_map = Map((mem_map.data, header)).rotate()
em_map = Map((em_map, header)).rotate()
vmax = max([clean_map.data.max(), mem_map.data.max(), em_map.data.max()])
2026-06-05T12:39:37Z INFO stixpy.imaging.em 181: Iteration: 25, StdDeV: 0.13928589714685566, C-stat: 0.08517368133346023
2026-06-05T12:39:37Z INFO stixpy.imaging.em 181: Iteration: 50, StdDeV: 0.030663641077421023, C-stat: 0.05480387386700314
2026-06-05T12:39:37Z INFO stixpy.imaging.em 181: Iteration: 75, StdDeV: 0.0059564106733533875, C-stat: 0.04546091584825557
2026-06-05T12:39:37Z INFO stixpy.imaging.em 181: Iteration: 100, StdDeV: 0.0017573115540636135, C-stat: 0.042483320318589514
2026-06-05T12:39:37Z INFO stixpy.imaging.em 181: Iteration: 125, StdDeV: 0.0007913896525914521, C-stat: 0.04112209848866932
Finally compare the images from each algorithm
fig = plt.figure(figsize=(12, 12))
ax = fig.subplot_mosaic(
[
["bp", "clean"],
["mem", "em"],
],
subplot_kw={"projection": clean_map},
)
a = bp_map.plot(axes=ax["bp"])
ax["bp"].set_title("Back Projection")
b = clean_map.plot(axes=ax["clean"])
ax["clean"].set_title("Clean")
c = mem_map.plot(axes=ax["mem"])
ax["mem"].set_title("MEM GE")
d = em_map.plot(axes=ax["em"])
ax["em"].set_title("EM")
fig.colorbar(a, ax=ax.values())
plt.show()
Total running time of the script: (0 minutes 31.608 seconds)