.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/gallery/imaging_demo.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_generated_gallery_imaging_demo.py: ============ 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 .. GENERATED FROM PYTHON SOURCE LINES 13-38 .. code-block:: Python 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__) .. GENERATED FROM PYTHON SOURCE LINES 39-40 Read science file as Product .. GENERATED FROM PYTHON SOURCE LINES 40-46 .. code-block:: Python 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 .. rst-class:: sphx-glr-script-out .. code-block:: none Files Downloaded: 0%| | 0/1 [00:00 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,_] .. GENERATED FROM PYTHON SOURCE LINES 47-48 Read background file as Product .. GENERATED FROM PYTHON SOURCE LINES 48-54 .. code-block:: Python 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 .. rst-class:: sphx-glr-script-out .. code-block:: none Files Downloaded: 0%| | 0/1 [00:00 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] .. GENERATED FROM PYTHON SOURCE LINES 55-56 Set time and energy ranges which will be considered for the science and the background file .. GENERATED FROM PYTHON SOURCE LINES 56-64 .. code-block:: Python 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 .. GENERATED FROM PYTHON SOURCE LINES 65-66 Create the meta pixel, A, B, C, D for the science and the background data .. GENERATED FROM PYTHON SOURCE LINES 66-75 .. code-block:: Python 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 ) .. GENERATED FROM PYTHON SOURCE LINES 76-77 Perform background subtraction .. GENERATED FROM PYTHON SOURCE LINES 77-86 .. code-block:: Python 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 ), } .. GENERATED FROM PYTHON SOURCE LINES 87-88 Create visibilities from the meta pixels .. GENERATED FROM PYTHON SOURCE LINES 88-91 .. code-block:: Python vis = create_visibility(meta_pixels_bkg_subtracted) .. GENERATED FROM PYTHON SOURCE LINES 92-93 Obtain the necessary ephemeris data .. GENERATED FROM PYTHON SOURCE LINES 93-103 .. code-block:: Python 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), ) .. rst-class:: sphx-glr-script-out .. code-block:: none Files Downloaded: 0%| | 0/1 [00:00 .. GENERATED FROM PYTHON SOURCE LINES 168-170 Estimate the flare location and plot on top of back projection map. Note the coordinates are automatically converted from the STIXImaging to Helioprojective .. GENERATED FROM PYTHON SOURCE LINES 170-179 .. code-block:: Python 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() .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/stixpy/checkouts/v0.3.0/examples/imaging_demo.py:177: UserWarning: The figure layout has changed to tight fig.tight_layout() .. GENERATED FROM PYTHON SOURCE LINES 180-181 Use estimated flare location to create more accurate visibilities .. GENERATED FROM PYTHON SOURCE LINES 181-197 .. code-block:: Python 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) .. GENERATED FROM PYTHON SOURCE LINES 198-200 Selected detectors 10 to 3 order by sub-collimator e.g. 10a, 10b, 10c, 9a, 9b, 9c .... .. GENERATED FROM PYTHON SOURCE LINES 200-203 .. code-block:: Python 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() .. GENERATED FROM PYTHON SOURCE LINES 204-205 Create an ``xrayvsion`` visibility object .. GENERATED FROM PYTHON SOURCE LINES 205-209 .. code-block:: Python cal_vis.meta["offset"] = max_stix vis10_3 = cal_vis[idx] .. GENERATED FROM PYTHON SOURCE LINES 210-211 Set up image parameters .. GENERATED FROM PYTHON SOURCE LINES 211-215 .. code-block:: Python imsize = [129, 129] * u.pixel # number of pixels of the map to reconstruct pixel = [2, 2] * u.arcsec / u.pixel # pixel size in arcsec .. GENERATED FROM PYTHON SOURCE LINES 216-217 Create a back projection image with natural weighting .. GENERATED FROM PYTHON SOURCE LINES 217-220 .. code-block:: Python bp_nat = vis_to_image(vis10_3, imsize, pixel_size=pixel) .. GENERATED FROM PYTHON SOURCE LINES 221-222 Create a back projection image with uniform weighting .. GENERATED FROM PYTHON SOURCE LINES 222-225 .. code-block:: Python bp_uni = vis_to_image(vis10_3, imsize, pixel_size=pixel, scheme="uniform") .. GENERATED FROM PYTHON SOURCE LINES 226-227 Create a `sunpy.map.Map` with back projection .. GENERATED FROM PYTHON SOURCE LINES 227-230 .. code-block:: Python bp_map = vis_to_map(vis10_3, imsize, pixel_size=pixel) .. GENERATED FROM PYTHON SOURCE LINES 231-232 Crete a clean image using the clean algorithm `vis_clean` .. GENERATED FROM PYTHON SOURCE LINES 232-240 .. code-block:: Python 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 ) .. rst-class:: sphx-glr-script-out .. code-block:: none 2026-06-05T12:48:17Z INFO xrayvision.clean 124: Iter: 0, strength: 1.930919100604466, location: (np.int64(59), np.int64(67)) 2026-06-05T12:48:17Z INFO xrayvision.clean 124: Iter: 25, strength: 0.47719851724349815, location: (np.int64(80), np.int64(54)) 2026-06-05T12:48:17Z INFO xrayvision.clean 145: Largest residual negative .. GENERATED FROM PYTHON SOURCE LINES 241-242 Create a sunpy map for the clean image in Helioprojective .. GENERATED FROM PYTHON SOURCE LINES 242-252 .. code-block:: Python 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, ) .. GENERATED FROM PYTHON SOURCE LINES 253-254 Crete a map using the MEM GE algorithm `mem` .. GENERATED FROM PYTHON SOURCE LINES 254-260 .. code-block:: Python 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) .. rst-class:: sphx-glr-script-out .. code-block:: none 2026-06-05T12:48:19Z INFO xrayvision.mem 159: Iter: 0, Chi2: 234.49830788341825 2026-06-05T12:48:19Z INFO xrayvision.mem 159: Iter: 25, Chi2: 24.333702665147236 2026-06-05T12:48:19Z INFO xrayvision.mem 159: Iter: 50, Chi2: 19.266467368879347 2026-06-05T12:48:19Z INFO xrayvision.mem 159: Iter: 75, Chi2: 17.124105167153616 2026-06-05T12:48:19Z INFO xrayvision.mem 159: Iter: 100, Chi2: 15.880549155749382 2026-06-05T12:48:19Z INFO xrayvision.mem 159: Iter: 125, Chi2: 15.057967256302211 2026-06-05T12:48:19Z INFO xrayvision.mem 159: Iter: 150, Chi2: 14.507035644240354 2026-06-05T12:48:19Z INFO xrayvision.mem 159: Iter: 175, Chi2: 14.118968596360151 2026-06-05T12:48:20Z INFO xrayvision.mem 159: Iter: 200, Chi2: 13.825428062747665 2026-06-05T12:48:20Z INFO xrayvision.mem 159: Iter: 225, Chi2: 13.592951927443597 2026-06-05T12:48:20Z INFO xrayvision.mem 159: Iter: 250, Chi2: 13.399299686673773 2026-06-05T12:48:20Z INFO xrayvision.mem 510: Iter: 0, Obj function: (241.44033193177933+0j) 2026-06-05T12:48:23Z INFO xrayvision.mem 510: Iter: 25, Obj function: (29.805142705550765+0j) 2026-06-05T12:48:26Z INFO xrayvision.mem 510: Iter: 50, Obj function: (28.996740805172074+0j) 2026-06-05T12:48:30Z INFO xrayvision.mem 510: Iter: 75, Obj function: (28.920500700663403+0j) 2026-06-05T12:48:33Z INFO xrayvision.mem 510: Iter: 100, Obj function: (28.911357252278577+0j) .. GENERATED FROM PYTHON SOURCE LINES 261-262 Crete a map using the EM algorithm `EM` .. GENERATED FROM PYTHON SOURCE LINES 262-280 .. code-block:: Python 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()]) .. rst-class:: sphx-glr-script-out .. code-block:: none 2026-06-05T12:48:36Z INFO stixpy.imaging.em 181: Iteration: 25, StdDeV: 0.13928589714685566, C-stat: 0.08517368133346023 2026-06-05T12:48:36Z INFO stixpy.imaging.em 181: Iteration: 50, StdDeV: 0.030663641077421023, C-stat: 0.05480387386700314 2026-06-05T12:48:36Z INFO stixpy.imaging.em 181: Iteration: 75, StdDeV: 0.0059564106733533875, C-stat: 0.04546091584825557 2026-06-05T12:48:36Z INFO stixpy.imaging.em 181: Iteration: 100, StdDeV: 0.0017573115540636135, C-stat: 0.042483320318589514 2026-06-05T12:48:36Z INFO stixpy.imaging.em 181: Iteration: 125, StdDeV: 0.0007913896525914521, C-stat: 0.04112209848866932 .. GENERATED FROM PYTHON SOURCE LINES 281-282 Finally compare the images from each algorithm .. GENERATED FROM PYTHON SOURCE LINES 282-301 .. code-block:: Python 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() .. image-sg:: /generated/gallery/images/sphx_glr_imaging_demo_002.png :alt: Back Projection, Clean, MEM GE, EM :srcset: /generated/gallery/images/sphx_glr_imaging_demo_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 31.850 seconds) .. _sphx_glr_download_generated_gallery_imaging_demo.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: imaging_demo.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: imaging_demo.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: imaging_demo.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_