.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/gallery/rhessi.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_rhessi.py: ====================================== RHESSI Visibility Imaging ====================================== Create images from RHESSI visibility data .. GENERATED FROM PYTHON SOURCE LINES 8-19 .. code-block:: Python import astropy.units as apu import matplotlib.pyplot as plt import numpy as np from astropy.io import fits from xrayvision.clean import vis_clean from xrayvision.imaging import vis_psf_map, vis_to_map from xrayvision.mem import mem, resistant_mean from xrayvision.visibility import Visibilities, VisMeta .. GENERATED FROM PYTHON SOURCE LINES 20-22 We will use `astropy.io.fits` to download and open the RHESSI visibility fits file .. GENERATED FROM PYTHON SOURCE LINES 22-27 .. code-block:: Python hdul = fits.open( "https://hesperia.gsfc.nasa.gov/rhessi_extras/visibility_fits_v2/2002/02/21/hsi_vis_20020221_2357_0054_46tx3e.fits" ) .. GENERATED FROM PYTHON SOURCE LINES 28-30 No lets extract the visibility data the first thing we will do is filter the visibilities by time. We will extract the 3rd integration time range. .. GENERATED FROM PYTHON SOURCE LINES 30-36 .. code-block:: Python times = np.unique(hdul[-1].data["TRANGE"], axis=0) time_index, _ = np.where(hdul[-1].data["TRANGE"] == times[6]) vis_data = hdul[3].data[time_index] .. GENERATED FROM PYTHON SOURCE LINES 37-38 Next lets filter by energy range in this case 12 - 25 keV .. GENERATED FROM PYTHON SOURCE LINES 38-42 .. code-block:: Python energy_index, _ = np.where(vis_data["ERANGE"] == [12.0, 25.0]) vis_data = vis_data[energy_index] .. GENERATED FROM PYTHON SOURCE LINES 43-45 Now lets filter by ISC or detector to remove possibly bad data in this case need to remove ISC 0 and 1. .. GENERATED FROM PYTHON SOURCE LINES 45-49 .. code-block:: Python vis_data = vis_data[vis_data["isc"] > 1] vis_data = vis_data[vis_data["obsvis"] != 0 + 0j] .. GENERATED FROM PYTHON SOURCE LINES 50-51 Now we can create the visibility object from the filtered visibilities. .. GENERATED FROM PYTHON SOURCE LINES 51-65 .. code-block:: Python meta = VisMeta({"vis_labels": vis_data["isc"]}) vunit = apu.Unit("photon/(cm**2 s)") vis = Visibilities( visibilities=vis_data["obsvis"] * vunit, u=vis_data["u"] / apu.arcsec, v=vis_data["v"] / apu.arcsec, phase_center=vis_data["xyoffset"][0] * apu.arcsec, meta=meta, amplitude_uncertainty=vis_data["sigamp"] * vunit, ) .. GENERATED FROM PYTHON SOURCE LINES 66-67 Lets have a look at the point spread function (PSF) or dirty beam .. GENERATED FROM PYTHON SOURCE LINES 67-70 .. code-block:: Python psf_map = vis_psf_map(vis, shape=(101, 101) * apu.pixel, pixel_size=1.5 * apu.arcsec / apu.pixel, scheme="uniform") .. GENERATED FROM PYTHON SOURCE LINES 71-73 We can now make an image using the back projection algorithm essentially and inverse Fourier transform of the visibilities. .. GENERATED FROM PYTHON SOURCE LINES 73-76 .. code-block:: Python backproj_map = vis_to_map(vis, shape=[101, 101] * apu.pixel, pixel_size=1.5 * apu.arcsec / apu.pix) .. GENERATED FROM PYTHON SOURCE LINES 77-80 Back projection contain many artifact due to the incomplete sampling of the u-v plane as a result various algorithms have been developed to remove or deconvolve this effect. CLEAN is one of the oldest and simplest, a CLEAN image can be made. .. GENERATED FROM PYTHON SOURCE LINES 80-94 .. code-block:: Python # vis_data_59 = vis_data[vis_data['isc'] > 3] # # vis_59 = Visibility(vis=vis_data_59['obsvis']*apu.Unit('ph/cm*s'), u=vis_data_59['u']/apu.arcsec, # v=vis_data_59['v']/apu.arcsec, offset=vis_data_59['xyoffset'][0]*apu.arcsec) clean_map, model_map, residual_map = vis_clean( vis, shape=[101, 101] * apu.pixel, pixel_size=[1.5, 1.5] * apu.arcsec / apu.pix, clean_beam_width=10 * apu.arcsec, niter=100, ) .. rst-class:: sphx-glr-script-out .. code-block:: none Max iterations reached .. GENERATED FROM PYTHON SOURCE LINES 95-96 MEM .. GENERATED FROM PYTHON SOURCE LINES 96-119 .. code-block:: Python # Compute percent_lambda # Loop through ISCs starting with 6-9, but if we don't have at least 2 vis, lower isc_min to include next one down, etc. isc_min = 6 nbig = 0 while isc_min >= 0 and nbig < 2: ibig = np.argwhere(vis.meta.vis_labels >= isc_min) nbig = len(ibig) isc_min = isc_min - 1 # If still don't have at least 2 vis, return -1, otherwise calculate mean (but reject points > sigma away from mean) if nbig < 2: snr_value = -1 else: snr_value, _ = resistant_mean((np.abs(vis.visibilities[ibig]) / vis.amplitude_uncertainty[ibig]).flatten(), 3) percent_lambda = 11.0 / (snr_value**2 + 383.0) mem_map = mem(vis, shape=[101, 101] * apu.pixel, pixel_size=[1.5, 1.5] * apu.arcsec / apu.pix) mem_map.plot() .. image-sg:: /generated/gallery/images/sphx_glr_rhessi_001.png :alt: 2025-04-04 15:14:26 :srcset: /generated/gallery/images/sphx_glr_rhessi_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/xrayvision/envs/stable/lib/python3.9/site-packages/sunpy/map/mapbase.py:892: SunpyMetadataWarning: Missing metadata for observation time, setting observation time to current time. Set the 'DATE-AVG' FITS keyword to prevent this warning. warn_metadata("Missing metadata for observation time, " INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase] /home/docs/checkouts/readthedocs.org/user_builds/xrayvision/envs/stable/lib/python3.9/site-packages/sunpy/map/mapbase.py:633: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer. For frame 'heliographic_stonyhurst' the following metadata is missing: dsun_obs,hgln_obs,hglt_obs For frame 'heliographic_carrington' the following metadata is missing: dsun_obs,crln_obs,crlt_obs obs_coord = self.observer_coordinate .. GENERATED FROM PYTHON SOURCE LINES 120-121 Comparison .. GENERATED FROM PYTHON SOURCE LINES 121-136 .. code-block:: Python fig = plt.figure(figsize=(10, 10)) fig.add_subplot(221, projection=psf_map) fig.add_subplot(222, projection=backproj_map) fig.add_subplot(223, projection=clean_map) fig.add_subplot(224, projection=mem_map) axs = fig.get_axes() psf_map.plot(axes=axs[0]) axs[0].set_title("PSF") backproj_map.plot(axes=axs[1]) axs[1].set_title("Back Projection") clean_map.plot(axes=axs[2]) axs[2].set_title("Clean") mem_map.plot(axes=axs[3]) axs[3].set_title("MEM") plt.show() .. image-sg:: /generated/gallery/images/sphx_glr_rhessi_002.png :alt: PSF, Back Projection, Clean, MEM :srcset: /generated/gallery/images/sphx_glr_rhessi_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/xrayvision/envs/stable/lib/python3.9/site-packages/sunpy/map/mapbase.py:892: SunpyMetadataWarning: Missing metadata for observation time, setting observation time to current time. Set the 'DATE-AVG' FITS keyword to prevent this warning. warn_metadata("Missing metadata for observation time, " INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase] /home/docs/checkouts/readthedocs.org/user_builds/xrayvision/envs/stable/lib/python3.9/site-packages/sunpy/map/mapbase.py:633: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer. For frame 'heliographic_stonyhurst' the following metadata is missing: dsun_obs,hgln_obs,hglt_obs For frame 'heliographic_carrington' the following metadata is missing: dsun_obs,crln_obs,crlt_obs obs_coord = self.observer_coordinate /home/docs/checkouts/readthedocs.org/user_builds/xrayvision/envs/stable/lib/python3.9/site-packages/sunpy/map/mapbase.py:892: SunpyMetadataWarning: Missing metadata for observation time, setting observation time to current time. Set the 'DATE-AVG' FITS keyword to prevent this warning. warn_metadata("Missing metadata for observation time, " INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase] /home/docs/checkouts/readthedocs.org/user_builds/xrayvision/envs/stable/lib/python3.9/site-packages/sunpy/map/mapbase.py:633: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer. For frame 'heliographic_stonyhurst' the following metadata is missing: dsun_obs,hgln_obs,hglt_obs For frame 'heliographic_carrington' the following metadata is missing: dsun_obs,crln_obs,crlt_obs obs_coord = self.observer_coordinate /home/docs/checkouts/readthedocs.org/user_builds/xrayvision/envs/stable/lib/python3.9/site-packages/sunpy/map/mapbase.py:892: SunpyMetadataWarning: Missing metadata for observation time, setting observation time to current time. Set the 'DATE-AVG' FITS keyword to prevent this warning. warn_metadata("Missing metadata for observation time, " INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase] /home/docs/checkouts/readthedocs.org/user_builds/xrayvision/envs/stable/lib/python3.9/site-packages/sunpy/map/mapbase.py:633: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer. For frame 'heliographic_stonyhurst' the following metadata is missing: dsun_obs,hgln_obs,hglt_obs For frame 'heliographic_carrington' the following metadata is missing: dsun_obs,crln_obs,crlt_obs obs_coord = self.observer_coordinate .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 32.837 seconds) .. _sphx_glr_download_generated_gallery_rhessi.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: rhessi.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: rhessi.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: rhessi.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_