.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/gallery/raster.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_raster.py: ==================================== Working with IRIS spectrograph files ==================================== In this example, we will showcase how to use ``irispy-lmsal`` to open, crop and plot IRIS spectrograph data. .. GENERATED FROM PYTHON SOURCE LINES 8-18 .. code-block:: Python import astropy.units as u import matplotlib.pyplot as plt import numpy as np import pooch from astropy.coordinates import SkyCoord, SpectralCoord from sunpy.coordinates.frames import Helioprojective from irispy.io import read_files .. GENERATED FROM PYTHON SOURCE LINES 19-24 We start with getting the data. This is done by downloading the data from the IRIS archive. In this case, we will use ``pooch`` to keep this example self contained but using your browser will also work. .. GENERATED FROM PYTHON SOURCE LINES 24-30 .. code-block:: Python raster_filename = pooch.retrieve( "http://www.lmsal.com/solarsoft/irisa/data/level2_compressed/2018/01/02/20180102_153155_3610108077/iris_l2_20180102_153155_3610108077_raster.tar.gz", known_hash="0ec2b7b20757c52b02e0d92c27a5852b6e28759512c3d455f8b6505d4e1f5cd6", ) .. GENERATED FROM PYTHON SOURCE LINES 31-34 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 34-41 .. code-block:: Python raster = read_files(raster_filename, memmap=True, uncertainty=False) # Provide an overview of the data print(raster) # Will give us all the keys that corresponds to wavelengths. print(raster.keys()) .. rst-class:: sphx-glr-script-out .. code-block:: none Collection ---------- Cube keys: ('C II 1336', 'O I 1356', 'Si IV 1394', 'Si IV 1403', '2832', '2814', 'Mg II k 2796') Number of Cubes: 7 Aligned dimensions: [ ] Aligned physical types: [('meta.obs.sequence',), ('custom:pos.helioprojective.lat', 'custom:pos.helioprojective.lon', 'time'), ('custom:pos.helioprojective.lat', 'custom:pos.helioprojective.lon')] dict_keys(['C II 1336', 'O I 1356', 'Si IV 1394', 'Si IV 1403', '2832', '2814', 'Mg II k 2796']) .. GENERATED FROM PYTHON SOURCE LINES 42-43 We will now access one wavelength window and plot it. .. GENERATED FROM PYTHON SOURCE LINES 43-51 .. code-block:: Python mg_ii = raster["Mg II k 2796"][0] plt.figure() mg_ii.plot() plt.show() .. rst-class:: sphx-glr-horizontal * .. image-sg:: /generated/gallery/images/sphx_glr_raster_001.png :alt: raster :srcset: /generated/gallery/images/sphx_glr_raster_001.png :class: sphx-glr-multi-img * .. image-sg:: /generated/gallery/images/sphx_glr_raster_002.png :alt: raster :srcset: /generated/gallery/images/sphx_glr_raster_002.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 52-53 We can index to get the first index and plot it. .. GENERATED FROM PYTHON SOURCE LINES 53-59 .. code-block:: Python plt.figure() mg_ii[0].plot() plt.show() .. image-sg:: /generated/gallery/images/sphx_glr_raster_003.png :alt: raster :srcset: /generated/gallery/images/sphx_glr_raster_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 60-61 Or we can get the spectrum directly. .. GENERATED FROM PYTHON SOURCE LINES 61-67 .. code-block:: Python plt.figure() mg_ii[120, 200].plot() plt.show() .. image-sg:: /generated/gallery/images/sphx_glr_raster_004.png :alt: raster :srcset: /generated/gallery/images/sphx_glr_raster_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 68-72 The default plots take the units and labels from the FITS WCS information, and often do not come in the most useful units (e.g. wavelengths in metres). We can read the wavelengths of the Mg window by calling axis_world_coords for wl (wavelength), and redo the plot with a better scale: .. GENERATED FROM PYTHON SOURCE LINES 72-82 .. code-block:: Python (mg_wave,) = mg_ii.axis_world_coords("wl") mg_wave.to("m") fig, ax = plt.subplots() # Here we index the data directly instead of other methods to get access to the data and plot it "raw" plt.plot(mg_wave.to("nm"), mg_ii.data[120, 200]) plt.show() .. image-sg:: /generated/gallery/images/sphx_glr_raster_005.png :alt: raster :srcset: /generated/gallery/images/sphx_glr_raster_005.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 83-90 When we use the underlying data directly, we lose all the metadata and WCS information, so we would suggest not doing it normally but there will be times you will need to do this. If you are unfamiliar with WCS, the following links are quite useful: The base comes from astropy: https://docs.astropy.org/en/stable/wcs/index.html The plotting makes use of WCSAxes: https://docs.astropy.org/en/stable/visualization/wcsaxes/index.html Some of the higher-level utils are via NDCube, e.g. coordinate transformations: https://docs.sunpy.org/projects/ndcube/en/v2.0.0rc2/coordinates.html. We can also improve on the default spectrogram plot by adjusting some options: .. GENERATED FROM PYTHON SOURCE LINES 90-103 .. code-block:: Python plt.figure() mg_ii[0].plot(aspect="auto", cmap="irissjiNUV") ax = plt.gca() ax.coords[0].set_major_formatter("x.x") ax.coords[0].set_format_unit(u.nm) ax.set_xlabel("Wavelength (nm)") ax.set_ylabel("Helioprojective latitude [Solar Y] (arcsec)") # Remove longitude ticks ax.coords[2].set_ticks([] * u.degree) plt.show() .. image-sg:: /generated/gallery/images/sphx_glr_raster_006.png :alt: raster :srcset: /generated/gallery/images/sphx_glr_raster_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 104-105 What is the wavelength position that corresponds to Mg II k core (279.63 nm)? .. GENERATED FROM PYTHON SOURCE LINES 105-113 .. code-block:: Python wcs_loc = mg_ii.wcs.world_to_pixel( SpectralCoord(279.63, unit=u.nm), SkyCoord(0 * u.arcsec, 0 * u.arcsec, frame=Helioprojective), ) mg_index = int(np.round(wcs_loc[0])) print(mg_index) .. 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( /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( 111 .. GENERATED FROM PYTHON SOURCE LINES 114-118 Now we will plot spectroheliogram for Mg II k core wavelength. We can use a new feature in ndcube 2.0 called crop to get this information that will require a SkyCoord and SpectralCoord object from astropy. In this case, we only need to use a SpectralCoord. .. GENERATED FROM PYTHON SOURCE LINES 118-129 .. code-block:: Python lower_corner = [SpectralCoord(280, unit=u.nm), None] upper_corner = [SpectralCoord(280, unit=u.nm), None] mg_crop = mg_ii.crop(lower_corner, upper_corner) ax = mg_crop.plot() plt.xlabel("Solar X") plt.ylabel("Solar Y") plt.show() .. image-sg:: /generated/gallery/images/sphx_glr_raster_007.png :alt: raster :srcset: /generated/gallery/images/sphx_glr_raster_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 130-132 Imagine there's a really cool feature at (-338", 275"), how can you plot the spectrum at that location? .. GENERATED FROM PYTHON SOURCE LINES 132-139 .. code-block:: Python lower_corner = [None, SkyCoord(-338 * u.arcsec, 275 * u.arcsec, frame=Helioprojective)] upper_corner = [None, SkyCoord(-338 * u.arcsec, 275 * u.arcsec, frame=Helioprojective)] mg_ii.crop(lower_corner, upper_corner).plot() plt.show() .. image-sg:: /generated/gallery/images/sphx_glr_raster_008.png :alt: raster :srcset: /generated/gallery/images/sphx_glr_raster_008.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 140-142 Now, you may also be interested in knowing the time that was this observation taken. There is some information in ``.meta``. .. GENERATED FROM PYTHON SOURCE LINES 142-145 .. code-block:: Python print(mg_ii.meta) .. rst-class:: sphx-glr-script-out .. code-block:: none SGMeta ------ Observatory: IRIS Instrument: SPEC Detector: NUV Spectral Window: Mg II k 2796 Spectral Range: [2790.65477674 2809.95345615] Angstrom Spectral Band: NUV Dimensions: [320, 548, 380] Date: 2018-01-02T15:31:55.870 OBS ID: 3610108077 OBS Description: Very large dense 320-step raster 105.3x175 320s Deep x 8 Spatial x .. GENERATED FROM PYTHON SOURCE LINES 146-149 But this is mostly about the observation in general. Times of individual scans are saved in .extra_coords['time']. Getting access to it can be done in the following way: .. GENERATED FROM PYTHON SOURCE LINES 149-152 .. code-block:: Python # The first index is to escape the tuple that ``axis_world_coords`` returns print(mg_ii.axis_world_coords("time", wcs=mg_ii.extra_coords)[0][0].isot) .. rst-class:: sphx-glr-script-out .. code-block:: none 2018-01-02T15:31:55.870 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 14.077 seconds) .. _sphx_glr_download_generated_gallery_raster.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: raster.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: raster.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_