.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/gallery/mg_ii_dopplergrams.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_mg_ii_dopplergrams.py: ================== Mg II Dopplergrams ================== In this example we are going to produce a Dopplergram for the Mg II k line from a IRIS 400-step raster. The Dopplergram is obtained by subtracting the intensities at symmetrical velocity shifts from the line core (e.g. ±50 km/s). This very large dense raster took more than three hours to complete the 400 scans (30 s exposures), which means that the spacecraft's orbital velocity changes during the observations. This means that any precise wavelength calibration will need to correct for those shifts. .. GENERATED FROM PYTHON SOURCE LINES 12-27 .. code-block:: Python import tarfile import astropy.units as u import matplotlib.pyplot as plt import numpy as np import pooch from astropy.coordinates import SpectralCoord from astropy.io import fits from scipy.constants import c from scipy.interpolate import interp1d from irispy.io import read_files from irispy.utils import image_clipping .. GENERATED FROM PYTHON SOURCE LINES 28-35 We start with getting the data. This is done by downloading the data from the IRIS archive. In this case, we will use ``pooch`` as to keep this example self contained but using your browser will also work. Using the url: http://www.lmsal.com/solarsoft/irisa/data/level2_compressed/2014/07/08/ we are after `iris_l2_20140708_114109_3824262996_raster.tar.gz` which is ~730 MB in size. .. GENERATED FROM PYTHON SOURCE LINES 35-41 .. code-block:: Python downloaded_tar_iris_file = pooch.retrieve( "http://www.lmsal.com/solarsoft/irisa/data/level2_compressed/2014/07/08/20140708_114109_3824262996/iris_l2_20140708_114109_3824262996_raster.tar.gz", known_hash=None, ) .. GENERATED FROM PYTHON SOURCE LINES 42-46 Now to open the file using ``irispy-lmsal``. Note that when ``memmap=True``, the data values are read from the FITS file directly without the scaling to Float32, the data values are no longer in DN, but in scaled integer units that start at -2$^{16}$/2. .. GENERATED FROM PYTHON SOURCE LINES 46-49 .. code-block:: Python raster = read_files(downloaded_tar_iris_file, memmap=True, uncertainty=False) .. GENERATED FROM PYTHON SOURCE LINES 50-51 We are after the Mg II k window, which we can select using a key. .. GENERATED FROM PYTHON SOURCE LINES 51-62 .. code-block:: Python mg_ii = raster["Mg II k 2796"][0] (mg_wave,) = mg_ii.axis_world_coords("wl") # We will plot the spatially averaged spectrum plt.plot(mg_wave.to("nm"), mg_ii.data.mean((0, 1))) plt.ylabel("DN (Memory Mapped Value)") plt.xlabel("Wavelength (nm)") plt.show() .. image-sg:: /generated/gallery/images/sphx_glr_mg_ii_dopplergrams_001.png :alt: mg ii dopplergrams :srcset: /generated/gallery/images/sphx_glr_mg_ii_dopplergrams_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/irispy-lmsal/conda/latest/lib/python3.10/site-packages/astropy/wcs/wcsapi/fitswcs.py:533: AstropyUserWarning: target cannot be converted to ICRS, so will not be set on SpectralCoord warnings.warn( .. GENERATED FROM PYTHON SOURCE LINES 63-68 To better understand the orbital velocity problem let us look at how the line intensity varies for a strong Mn I line at around 280.2 nm, in between the Mg II k and h lines. For this dataset, the line core of this line falls around index 350. To plot a spectrogram in the correct orientation we will transpose the data. .. GENERATED FROM PYTHON SOURCE LINES 68-87 .. code-block:: Python lower_corner = [SpectralCoord(280.2, unit=u.nm), None] upper_corner = [SpectralCoord(280.2, unit=u.nm), None] mg_crop = mg_ii.crop(lower_corner, upper_corner) vmin, vmax = image_clipping(mg_ii.data[..., 350]) plt.figure(figsize=(6, 10)) plt.imshow( mg_ii.data[..., 350].T, origin="lower", aspect=0.5, vmin=vmin, vmax=vmax, ) plt.xlabel("Solar X (arcsec)") plt.ylabel("Solar Y (arcsec)") plt.show() .. image-sg:: /generated/gallery/images/sphx_glr_mg_ii_dopplergrams_002.png :alt: mg ii dopplergrams :srcset: /generated/gallery/images/sphx_glr_mg_ii_dopplergrams_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/irispy-lmsal/conda/latest/lib/python3.10/site-packages/astropy/wcs/wcsapi/fitswcs.py:671: AstropyUserWarning: No observer defined on WCS, SpectralCoord will be converted without any velocity frame change warnings.warn( .. GENERATED FROM PYTHON SOURCE LINES 88-94 You can see a regular bright-dark pattern along the Y axis, an indication that its intensities are not taken at the same position in the line because of wavelength shifts. The shifts are caused by the orbital velocity changes, and we can find these in the auxiliary metadata which are to be found in the extension past the "last" window in the FITS file. .. GENERATED FROM PYTHON SOURCE LINES 94-115 .. code-block:: Python # astropy.io.fits does not support opening tar files, so we need to extract. tar_iris_file = tarfile.open(downloaded_tar_iris_file, "r:gz") tar_iris_file.extractall("./") tar_iris_file.close() raster_filename = "iris_l2_20140708_114109_3824262996_raster_t000_r00000.fits" # In this case, it is the 9th HDU, which we access directly aux_data = fits.getdata(raster_filename, 9) aux_header = fits.getheader(raster_filename, 9) v_obs = aux_data[:, aux_header["OBS_VRIX"]] # Convert to km/s v_obs /= 1000.0 plt.figure() plt.plot(v_obs) plt.ylabel("Orbital velocity (km/s)") plt.xlabel("Scan number") plt.show() .. image-sg:: /generated/gallery/images/sphx_glr_mg_ii_dopplergrams_003.png :alt: mg ii dopplergrams :srcset: /generated/gallery/images/sphx_glr_mg_ii_dopplergrams_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 116-121 To look at intensities at any given scan we only need to subtract this velocity shift from the wavelength scale, but to look at the whole image at a given wavelength we must interpolate the original data to take this shift into account. Here is a way to do it (note that array dimensions apply to this specific example). .. GENERATED FROM PYTHON SOURCE LINES 121-129 .. code-block:: Python c_kms = c / 1000.0 wave_shift = -v_obs * mg_wave[350] / (c_kms) # Linear interpolation in wavelength, for each scan for i in range(mg_ii.data.shape[0]): method = interp1d(mg_wave - wave_shift[i], mg_ii.data[:, i, :], bounds_error=False) mg_ii.data[:, i, :] = method(mg_wave) .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/irispy-lmsal/checkouts/latest/examples/mg_ii_dopplergrams.py:127: RuntimeWarning: invalid value encountered in cast mg_ii.data[:, i, :] = method(mg_wave) .. GENERATED FROM PYTHON SOURCE LINES 130-132 Now we can plot the shifted data to see that the large scale shifts have disappeared .. GENERATED FROM PYTHON SOURCE LINES 132-147 .. code-block:: Python vmin, vmax = image_clipping(mg_ii.data[..., 350]) plt.figure(figsize=(6, 10)) plt.imshow( mg_ii.data[..., 350].T, origin="lower", aspect=0.5, vmin=vmin, vmax=vmax, ) plt.xlabel("Solar X (arcsec)") plt.ylabel("Solar Y (arcsec)") plt.show() .. image-sg:: /generated/gallery/images/sphx_glr_mg_ii_dopplergrams_004.png :alt: mg ii dopplergrams :srcset: /generated/gallery/images/sphx_glr_mg_ii_dopplergrams_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 148-162 Some residual shift remains, but we will not correct for it here. A more elaborate correction can be obtained by the IDL routine ``iris_prep_wavecorr_l2``, but this has not yet been ported to Python see the `IDL version of this tutorial `__ for more details. We can use the calibrated data for example to calculate Dopplergrams. A Dopplergram is here defined as the difference between the intensities at two wavelength positions at the same (and opposite) distance from the line core. For example, at +/- 50 km/s from the Mg II k3 core. To do this, let us first calculate a velocity scale for the k line and find the indices of the -50 and +50 km/s velocity positions (here using the convention of negative velocities for up flows): .. GENERATED FROM PYTHON SOURCE LINES 162-170 .. code-block:: Python mg_k_centre = 279.6351 * u.nm pos = 50 # in km/s around line centre velocity = (mg_wave - mg_k_centre) * c_kms / mg_k_centre index_p = np.argmin(np.abs(velocity - pos)) index_m = np.argmin(np.abs(velocity + pos)) doppl = mg_ii.data[..., index_m] - mg_ii.data[..., index_p] .. GENERATED FROM PYTHON SOURCE LINES 171-173 And now we can plot this as before (intensity units are again arbitrary because of the unscaled DNs): .. GENERATED FROM PYTHON SOURCE LINES 173-188 .. code-block:: Python vmin, vmax = image_clipping(doppl) plt.figure(figsize=(6, 10)) plt.imshow( doppl.T, cmap="gist_gray", origin="lower", aspect=0.5, vmin=vmin, vmax=vmax, ) plt.xlabel("Solar X (arcsec)") plt.ylabel("Solar Y (arcsec)") plt.show() .. image-sg:: /generated/gallery/images/sphx_glr_mg_ii_dopplergrams_005.png :alt: mg ii dopplergrams :srcset: /generated/gallery/images/sphx_glr_mg_ii_dopplergrams_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 43.326 seconds) .. _sphx_glr_download_generated_gallery_mg_ii_dopplergrams.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: mg_ii_dopplergrams.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: mg_ii_dopplergrams.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_