This notebook shows fitted EIS line intensities in a polar coronal hole and the quiet Sun. Link to Figure 5.

(The internal hyperlink only works on GitHub Pages or nbviewer. Do not click when viewing the notebook on GitHub.)

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from matplotlib import rcParams
from matplotlib import ticker
from matplotlib import patches
rcParams["figure.figsize"] = (10,8)
plt.style.use("science")
import eispac
from glob import glob
import sunpy
import sunpy.map
from sunpy.io.special import read_genx
from sunpy.image.coalignment import calculate_shift
import cmcrameri.cm as cmcm
import matplotlib.colors as colors
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from astropy.visualization import (ImageNormalize, SqrtStretch, LogStretch,
                                   ZScaleInterval, AsinhStretch)
from scipy import interpolate
from scipy import ndimage
import juanfit
import importlib
importlib.reload(juanfit)
from juanfit import SpectrumFitSingle, SpectrumFitRow, SpectrumFit2D
rcParams['axes.labelsize'] = 18
rcParams['xtick.labelsize'] = 16
rcParams['ytick.labelsize'] = 16
rcParams['font.size'] = 18
rcParams['figure.figsize'] = [10,10]
rcParams['axes.linewidth'] = 2
rcParams['axes.titlesize'] = 18
from scipy.interpolate import interp1d
import astropy.constants as const
import copy
import astropy.units as u
from astropy.coordinates import SkyCoord
In [ ]:
def plot_colorbar(im, ax, width="3%", height="100%",loc="lower left",fontsize=14,
                bbox_to_anchor=(1.02, 0., 1, 1)):
    clb_ax = inset_axes(ax,width=width,height=height,loc=loc,
                bbox_to_anchor=bbox_to_anchor,
                 bbox_transform=ax.transAxes,
                 borderpad=0)
    clb = plt.colorbar(im,pad = 0.05,orientation='vertical',ax=ax,cax=clb_ax)
    clb_ax.yaxis.set_minor_locator(ticker.AutoMinorLocator(5))
    clb_ax.yaxis.get_offset_text().set_fontsize(fontsize)
    clb_ax.tick_params(labelsize=fontsize)
    return clb, clb_ax
In [ ]:
ch_fex_184_file = "../../sav/EIS/NPCHDB/eispac_fit/eis_20170821_110818.fe_10_184_536.1c-0.fit.h5"
qs_fex_184_file = "../../sav/EIS/EQSPY/eispac_fit/eis_20170821_205401.fe_10_184_536.1c-0.fit.h5"
ch_fexii_195_file = "../../sav/EIS/NPCHDB/eispac_fit/eis_20170821_110818.fe_12_195_119.1c-0.fit.h5"
qs_fexii_195_file = "../../sav/EIS/EQSPY/eispac_fit/eis_20170821_205401.fe_12_195_119.1c-0.fit.h5"
ch_feviii_185_file = "../../sav/EIS/NPCHDB/eispac_fit/eis_20170821_110818.fe_08_185_213.1c-0.fit.h5"
qs_fexiv_264_file = "../../sav/EIS/EQSPY/eispac_fit/eis_20170821_205401.fe_14_264_787.1c-0.fit.h5"
In [ ]:
ch_fex_184_fitres = eispac.read_fit(ch_fex_184_file) 
ch_fex_184_intmap = ch_fex_184_fitres.get_map(component=0, measurement='intensity')
qs_fex_184_fitres = eispac.read_fit(qs_fex_184_file)
qs_fex_184_intmap = qs_fex_184_fitres.get_map(component=0, measurement='intensity')
ch_fexii_195_fitres = eispac.read_fit(ch_fexii_195_file)
ch_fexii_195_intmap = ch_fexii_195_fitres.get_map(component=0, measurement='intensity')
qs_fexii_195_fitres = eispac.read_fit(qs_fexii_195_file)
qs_fexii_195_intmap = qs_fexii_195_fitres.get_map(component=0, measurement='intensity')
ch_feviii_185_fitres = eispac.read_fit(ch_feviii_185_file)
ch_feviii_185_intmap = ch_feviii_185_fitres.get_map(component=0, measurement='intensity')
qs_fexiv_264_fitres = eispac.read_fit(qs_fexiv_264_file)
qs_fexiv_264_intmap = qs_fexiv_264_fitres.get_map(component=0, measurement='intensity')
Reading fit result from, 
   ../../sav/EIS/NPCHDB/eispac_fit/eis_20170821_110818.fe_10_184_536.1c-0.fit.h5
INFO: uncertainty should have attribute uncertainty_type. [astropy.nddata.nddata]
Reading fit result from, 
   ../../sav/EIS/EQSPY/eispac_fit/eis_20170821_205401.fe_10_184_536.1c-0.fit.h5
INFO: uncertainty should have attribute uncertainty_type. [astropy.nddata.nddata]
Reading fit result from, 
   ../../sav/EIS/NPCHDB/eispac_fit/eis_20170821_110818.fe_12_195_119.1c-0.fit.h5
INFO: uncertainty should have attribute uncertainty_type. [astropy.nddata.nddata]
Reading fit result from, 
   ../../sav/EIS/EQSPY/eispac_fit/eis_20170821_205401.fe_12_195_119.1c-0.fit.h5
INFO: uncertainty should have attribute uncertainty_type. [astropy.nddata.nddata]
Reading fit result from, 
   ../../sav/EIS/NPCHDB/eispac_fit/eis_20170821_110818.fe_08_185_213.1c-0.fit.h5
INFO: uncertainty should have attribute uncertainty_type. [astropy.nddata.nddata]
Reading fit result from, 
   ../../sav/EIS/EQSPY/eispac_fit/eis_20170821_205401.fe_14_264_787.1c-0.fit.h5
INFO: uncertainty should have attribute uncertainty_type. [astropy.nddata.nddata]
In [ ]:
qs_shiftx, qs_shifty = 4.5*u.arcsec, -4*u.arcsec
ch_shiftx, ch_shifty = -4.09618*u.arcsec, -2.50511*u.arcsec

ch_fex_184_intmap_shifted = ch_fex_184_intmap.shift(ch_shiftx, ch_shifty)
qs_fex_184_intmap_shifted = qs_fex_184_intmap.shift(qs_shiftx, qs_shifty)
ch_fexii_195_intmap_shifted = ch_fexii_195_intmap.shift(ch_shiftx, ch_shifty)
qs_fexii_195_intmap_shifted = qs_fexii_195_intmap.shift(qs_shiftx, qs_shifty)
ch_feviii_185_intmap_shifted = ch_feviii_185_intmap.shift(ch_shiftx, ch_shifty)
qs_fexiv_264_intmap_shifted = qs_fexiv_264_intmap.shift(qs_shiftx, qs_shifty)
In [ ]:
def plot_eismap(eismap,ax,cmap,xticklabel=True,yticklabel=False,norm=None,
                linename=None):
    if norm is None:
        norm = ImageNormalize(eismap.data, interval=ZScaleInterval(), stretch=AsinhStretch())
    ax.imshow(eismap.data,aspect=eismap.meta["cdelt2"]/eismap.meta["cdelt1"],
              cmap=cmap,norm=norm)
    
    ax.set_xlabel(" ")
    ax.set_ylabel(" ")

    ax.tick_params(which="both",direction="in")
    ax.tick_params(which="major",length=6)
    ax.tick_params(which="minor",length=4)
    ax.tick_params(axis="x",labelbottom=xticklabel)
    ax.tick_params(axis="y",labelleft=yticklabel)
    ax.coords[0].set_ticks(number=4)

    if linename is not None:
        title = r"\textbf{{EIS}} {:}".format(linename)+"\n"+\
                r"\textbf{{{:}}}".format(eismap.date_average.strftime("%Y-%m-%d %H:%M:%S"))
        title = ax.set_title(title,fontsize=17)
        title.set_multialignment("center")
In [ ]:
fig = plt.figure(figsize=(7,8),constrained_layout=True)

ax1 = fig.add_subplot(2,3,1,projection=ch_feviii_185_intmap_shifted)
plot_eismap(ch_feviii_185_intmap, ax1, plt.get_cmap("sohoeit171"),yticklabel=True,
            linename=r"\textbf{{Fe \textsc{{viii}} 18.52 nm}}")

ax2 = fig.add_subplot(2,3,2,projection=ch_fex_184_intmap_shifted)
plot_eismap(ch_fex_184_intmap, ax2, plt.get_cmap("sdoaia171"),
            linename=r"\textbf{{Fe \textsc{{x}} 18.45 nm}}")

ax3 = fig.add_subplot(2,3,3,projection=ch_fexii_195_intmap_shifted)
plot_eismap(ch_fexii_195_intmap, ax3, plt.get_cmap("sdoaia193"),
            linename=r"\textbf{{Fe \textsc{{xii}} 19.51 nm}}")

ax4 = fig.add_subplot(2,3,4,projection=qs_fex_184_intmap_shifted)
plot_eismap(qs_fex_184_intmap, ax4, plt.get_cmap("sdoaia171"),yticklabel=True,
            linename=r"\textbf{{Fe \textsc{{x}} 18.45 nm}}")

ax5 = fig.add_subplot(2,3,5,projection=qs_fexii_195_intmap_shifted)
plot_eismap(qs_fexii_195_intmap, ax5, plt.get_cmap("sdoaia193"),
            norm=ImageNormalize(qs_fexii_195_intmap.data, interval=ZScaleInterval(),
                                stretch=AsinhStretch(0.5)),
            linename=r"\textbf{{Fe \textsc{{xii}} 19.51 nm}}")

ax6 = fig.add_subplot(2,3,6,projection=qs_fexiv_264_intmap_shifted)
plot_eismap(qs_fexiv_264_intmap, ax6, plt.get_cmap("sdoaia211"),
            linename=r"\textbf{{Fe \textsc{{xiv}} 26.48 nm}}")

for ax_ in (ax1,ax2,ax3):
    ax_.axvline(77,color="red",alpha=1,zorder=6,linewidth=2)

for ax_ in (ax4,ax5,ax6):
    ax_.add_patch(patches.Rectangle((qs_fexiv_264_intmap_shifted.bottom_left_coord.Tx.to("deg").value,
                                     qs_fexiv_264_intmap_shifted.bottom_left_coord.Ty.to("deg").value),32/3600,
                                     qs_fex_184_intmap_shifted.top_right_coord.Ty.to("deg").value - \
                                     qs_fexiv_264_intmap_shifted.bottom_left_coord.Ty.to("deg").value,linewidth=2,edgecolor="red",
                            facecolor="none",alpha=1,zorder=6,transform=ax_.get_transform("world")))
    



ax1.set_ylabel(r"\textbf{SOLAR-Y [arcsec]}",fontsize=17)
ax4.set_ylabel(r"\textbf{SOLAR-Y [arcsec]}",fontsize=17)
ax4.set_xlabel(r"\textbf{SOLAR-X [arcsec]}",fontsize=17)
ax5.set_xlabel(r"\textbf{SOLAR-X [arcsec]}",fontsize=17)
ax6.set_xlabel(r"\textbf{SOLAR-X [arcsec]}",fontsize=17)

# plt.savefig(fname="../../figs/ms/eis_quicklook.pdf",format="pdf",dpi=300,
#             bbox_inches="tight")
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
In [ ]:
ch_select_xcoord = np.linspace(0,154,78) + ch_feviii_185_intmap_shifted.bottom_left_coord.Tx.value
ch_select_ycoord_098 = np.sqrt((950.*0.98)**2 - ch_select_xcoord**2)
ch_select_ycoord_103 = np.sqrt((950.*1.03)**2 - ch_select_xcoord**2)
ch_select_ycoord_108 = np.sqrt((950.*1.08)**2 - ch_select_xcoord**2)
ch_select_ycoord_113 = np.sqrt((950.*1.13)**2 - ch_select_xcoord**2)
ch_select_ycoord_118 = np.sqrt((950.*1.18)**2 - ch_select_xcoord**2)
In [ ]:
qs_select_xcoord = np.linspace(0,30,16) + qs_fex_184_intmap_shifted.bottom_left_coord.Tx.value
qs_select_xcoord_bottom = np.linspace(5,30,26) + qs_fex_184_intmap_shifted.bottom_left_coord.Tx.value
qs_select_ycoord_1035 = np.sqrt((950.*1.035)**2 - qs_select_xcoord_bottom**2)
qs_select_ycoord_106 = np.sqrt((950.*1.06)**2 - qs_select_xcoord**2)
qs_select_ycoord_110 = np.sqrt((950.*1.10)**2 - qs_select_xcoord**2)

Figure 5¶

back to top

In [ ]:
fig = plt.figure(figsize=(7,8),constrained_layout=True)

ax1 = fig.add_subplot(2,3,1,projection=ch_feviii_185_intmap_shifted)
plot_eismap(ch_feviii_185_intmap, ax1, plt.get_cmap("sohoeit171"),yticklabel=True,
            linename=r"\textbf{{Fe \textsc{{viii}} 18.52 nm}}")

ax2 = fig.add_subplot(2,3,2,projection=ch_fex_184_intmap_shifted)
plot_eismap(ch_fex_184_intmap, ax2, plt.get_cmap("sdoaia171"),
            linename=r"\textbf{{Fe \textsc{{x}} 18.45 nm}}")

ax3 = fig.add_subplot(2,3,3,projection=ch_fexii_195_intmap_shifted)
plot_eismap(ch_fexii_195_intmap, ax3, plt.get_cmap("sdoaia193"),
            linename=r"\textbf{{Fe \textsc{{xii}} 19.51 nm}}")

ax4 = fig.add_subplot(2,3,4,projection=qs_fex_184_intmap_shifted)
plot_eismap(qs_fex_184_intmap, ax4, plt.get_cmap("sdoaia171"),yticklabel=True,
            linename=r"\textbf{{Fe \textsc{{x}} 18.45 nm}}")

ax5 = fig.add_subplot(2,3,5,projection=qs_fexii_195_intmap_shifted)
plot_eismap(qs_fexii_195_intmap, ax5, plt.get_cmap("sdoaia193"),
            norm=ImageNormalize(qs_fexii_195_intmap.data, interval=ZScaleInterval(),
                                stretch=AsinhStretch(0.5)),
            linename=r"\textbf{{Fe \textsc{{xii}} 19.51 nm}}")

ax6 = fig.add_subplot(2,3,6,projection=qs_fexiv_264_intmap_shifted)
plot_eismap(qs_fexiv_264_intmap, ax6, plt.get_cmap("sdoaia211"),
            linename=r"\textbf{{Fe \textsc{{xiv}} 26.48 nm}}")

for ax_ in (ax1,ax2,ax3):
    ax_.plot([ch_select_xcoord[-1]/3600.,ch_select_xcoord[-1]/3600.],
             [ch_select_ycoord_103[-1]/3600.,ch_select_ycoord_118[-1]/3600.],
             color="red",alpha=1,zorder=6,linewidth=2,transform=ax_.get_transform("world"))
    ax_.plot([ch_select_xcoord[0]/3600.,ch_select_xcoord[0]/3600.],
             [ch_select_ycoord_103[0]/3600.,ch_select_ycoord_118[0]/3600.],
             color="red",alpha=1,zorder=6,linewidth=2,transform=ax_.get_transform("world"))
    
    ax_.plot(ch_select_xcoord/3600.,ch_select_ycoord_103/3600.,transform=ax_.get_transform("world"),
             linewidth=2,color="red")
    ax_.plot(ch_select_xcoord/3600.,ch_select_ycoord_108/3600.,transform=ax_.get_transform("world"),
             linewidth=2,color="red")
    ax_.plot(ch_select_xcoord/3600.,ch_select_ycoord_113/3600.,transform=ax_.get_transform("world"),
             linewidth=2,color="red")
    ax_.plot(ch_select_xcoord/3600.,ch_select_ycoord_118/3600.,transform=ax_.get_transform("world"),
             linewidth=2,color="red")
    ax_.plot(ch_select_xcoord/3600.,ch_select_ycoord_098/3600.,transform=ax_.get_transform("world"),
             linewidth=2,color="#86C166")
    
    ax_.plot([ch_select_xcoord[-1]/3600.,ch_select_xcoord[-1]/3600.],
             [ch_fexii_195_intmap_shifted.bottom_left_coord.Ty.value/3600.,ch_select_ycoord_098[-1]/3600.],
             color="#86C166",alpha=1,zorder=6,linewidth=2,transform=ax_.get_transform("world"))
    ax_.plot([ch_select_xcoord[0]/3600.,ch_select_xcoord[-1]/3600.],
             [ch_fexii_195_intmap_shifted.bottom_left_coord.Ty.value/3600.,ch_fexii_195_intmap_shifted.bottom_left_coord.Ty.value/3600.],
             color="#86C166",alpha=1,zorder=6,linewidth=2,transform=ax_.get_transform("world"))
    ax_.plot([ch_select_xcoord[0]/3600.,ch_select_xcoord[0]/3600.],
             [ch_fexii_195_intmap_shifted.bottom_left_coord.Ty.value/3600.,ch_select_ycoord_098[0]/3600.],
             color="#86C166",alpha=1,zorder=6,linewidth=2,transform=ax_.get_transform("world"))
    
for ax_ in (ax4,ax5,ax6):
    ax_.plot(qs_select_xcoord_bottom/3600.,qs_select_ycoord_1035/3600.,transform=ax_.get_transform("world"),
             linewidth=2,color="red")
    ax_.plot(qs_select_xcoord/3600.,qs_select_ycoord_106/3600.,transform=ax_.get_transform("world"),
             linewidth=2,color="red")
    
    ax_.plot([qs_select_xcoord[-1]/3600.,qs_select_xcoord[-1]/3600.],
            [qs_select_ycoord_1035[-1]/3600.,qs_fex_184_intmap_shifted.top_right_coord.Ty.to("deg").value],
            color="red",alpha=1,zorder=6,linewidth=2,transform=ax_.get_transform("world"))
    ax_.plot([qs_select_xcoord[0]/3600.,qs_select_xcoord[0]/3600.],
            [qs_select_ycoord_1035[0]/3600.,qs_fex_184_intmap_shifted.top_right_coord.Ty.to("deg").value],
            color="red",alpha=1,zorder=6,linewidth=2,transform=ax_.get_transform("world"))
    ax_.plot([qs_select_xcoord[0]/3600.,qs_select_xcoord[-1]/3600.],
            [qs_fex_184_intmap_shifted.top_right_coord.Ty.to("deg").value,qs_fex_184_intmap_shifted.top_right_coord.Ty.to("deg").value],
            color="red",alpha=1,zorder=6,linewidth=2,transform=ax_.get_transform("world"))
    ax_.plot([qs_select_xcoord[0]/3600.,qs_select_xcoord_bottom[0]/3600.],
            [qs_fexiv_264_intmap_shifted.bottom_left_coord.Ty.to("deg").value,qs_fexiv_264_intmap_shifted.bottom_left_coord.Ty.to("deg").value],
            color="red",alpha=1,zorder=6,linewidth=2,transform=ax_.get_transform("world"))
    # ax_.add_patch(patches.Rectangle((qs_fexiv_264_intmap_shifted.bottom_left_coord.Tx.to("deg").value,
    #                                  qs_fexiv_264_intmap_shifted.bottom_left_coord.Ty.to("deg").value),30/3600,
    #                                  qs_fex_184_intmap_shifted.top_right_coord.Ty.to("deg").value - \
    #                                  qs_fexiv_264_intmap_shifted.bottom_left_coord.Ty.to("deg").value,linewidth=2,edgecolor="red",
    #                         facecolor="none",alpha=1,zorder=6,transform=ax_.get_transform("world")))
    


ax1.set_ylabel(r"\textbf{SOLAR-Y [arcsec]}",fontsize=17)
ax4.set_ylabel(r"\textbf{SOLAR-Y [arcsec]}",fontsize=17)
ax4.set_xlabel(r"\textbf{SOLAR-X [arcsec]}",fontsize=17)
ax5.set_xlabel(r"\textbf{SOLAR-X [arcsec]}",fontsize=17)
ax6.set_xlabel(r"\textbf{SOLAR-X [arcsec]}",fontsize=17)

# plt.savefig(fname="../../figs/ms/eis_quicklook.pdf",format="pdf",dpi=300,
#             bbox_inches="tight")
In [ ]:
ax4.wcs
Out[ ]:
WCS Keywords

Number of WCS axes: 2
CTYPE : 'HPLN-TAN'  'HPLT-TAN'  
CRVAL : -0.25882405784395  0.086270489162869  
CRPIX : 1.0  1.0  
PC1_1 PC1_2  : 1.0  0.0  
PC2_1 PC2_2  : 0.0  1.0  
CDELT : 0.00055466665161981  0.00027777777777778  
NAXIS : 60  160
In [ ]:
qs_fexiv_264_intmap_shifted.bottom_left_coord.Ty.to("deg").value
Out[ ]:
0.09117479960123376
In [ ]:
qs_fexiv_264_intmap_shifted.top_right_coord.Ty - qs_fex_184_intmap_shifted.top_right_coord.Ty
Out[ ]:
$17.6555\mathrm{{}^{\prime\prime}}$
In [ ]:
qs_fexiv_264_intmap_shifted.top_right_coord.Ty - qs_fexii_195_intmap_shifted.top_right_coord.Ty
Out[ ]:
$15.8477\mathrm{{}^{\prime\prime}}$