.. 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-28 .. 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 29-31 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 31-37 .. 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 38-39 Next lets filter by energy range in this case 12 - 25 keV .. GENERATED FROM PYTHON SOURCE LINES 39-43 .. 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 44-46 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 46-50 .. 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 51-52 Now we can create the visibility object from the filtered visibilities. .. GENERATED FROM PYTHON SOURCE LINES 52-66 .. 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 67-68 Lets have a look at the point spread function (PSF) or dirty beam .. GENERATED FROM PYTHON SOURCE LINES 68-71 .. 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 72-74 We can now make an image using the back projection algorithm essentially and inverse Fourier transform of the visibilities. .. GENERATED FROM PYTHON SOURCE LINES 74-77 .. 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 78-81 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 81-95 .. 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 96-97 MEM .. GENERATED FROM PYTHON SOURCE LINES 97-120 .. 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: 2024-07-16 20:43:20 :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/v0.2.0/lib/python3.9/site-packages/astropy/units/quantity.py:1331: ComplexWarning: Casting complex values to real discards the imaginary part self.view(np.ndarray).__setitem__(i, self._to_own_unit(value)) /home/docs/checkouts/readthedocs.org/user_builds/xrayvision/envs/v0.2.0/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/v0.2.0/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 121-122 Comparison .. GENERATED FROM PYTHON SOURCE LINES 122-137 .. 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/v0.2.0/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/v0.2.0/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/v0.2.0/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/v0.2.0/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/v0.2.0/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/v0.2.0/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 46.478 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 ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_