This notebook plots the magnetic field lines extrapolated by the Potential Field Source Surface (PFSS) model in sunkit_magex. The field lines are traced nearby the east and west upflow regions and overplotted on a SDO/AIA image and a SDO/HMI daily squashing factor Q map. Link to Figure 6.
Note: The internal hyperlink only works on GitHub Pages or nbviewer. Do not click when viewing the notebook on GitHub.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc_context
import matplotlib.patheffects as path_effects
import sunpy
import sunpy.map
from sunpy.coordinates import propagate_with_solar_surface
import astropy
from astropy.coordinates import SkyCoord
import astropy.units as u
import astropy.constants as const
from astropy.io import fits
from astropy.time import Time
from astropy.convolution import convolve, Gaussian2DKernel
from sunkit_magex import pfss as pfsspy
import pfsspy.utils
import pfsspy.tracing as tracing
import cmcrameri.cm as cmcm
from matplotlib.ticker import (AutoLocator, AutoMinorLocator,
FixedLocator, FixedFormatter, LogLocator, StrMethodFormatter)
from astropy.visualization import (AsinhStretch, LinearStretch,
LogStretch, ImageNormalize)
from regions import (CircleSkyRegion, Regions, EllipseSkyRegion)
from fancy_colorbar import plot_colorbar
from PIL import Image
In [2]:
ms_style_dict = {'text.usetex': True, 'font.family': 'serif', 'axes.linewidth': 1.2,
'xtick.major.width': 1.2, 'xtick.major.size': 4,
'ytick.major.width': 1.2, 'ytick.major.size': 4,
'xtick.minor.width': 1.2, 'xtick.minor.size': 2,
'ytick.minor.width': 1.2, 'ytick.minor.size': 2,
'xtick.direction': 'in', 'ytick.direction': 'in',
'text.latex.preamble': r'\usepackage[T1]{fontenc}'
r'\usepackage{amsmath}' r'\usepackage{siunitx}'
r'\sisetup{detect-all=True}'}
In [3]:
eis_195_velmap_shift_1025 = sunpy.map.Map("../../src/EIS/DHB_007_v2/20221025T2011/sunpymaps/eis_195_velmap_shift.fits")
eis_195_velmap_shift_1020 = sunpy.map.Map("../../src/EIS/DHB_007_v2/20221020T2343/sunpymaps/eis_195_velmap_shift.fits")
In [4]:
eis_east_region = Regions.read("../../sav/regions/eis_1025_east_pixel.reg")[0].to_sky(eis_195_velmap_shift_1025.wcs)
eis_west_region = Regions.read("../../sav/regions/eis_1020_west_pixel.reg")[0].to_sky(eis_195_velmap_shift_1020.wcs)
WARNING: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer. For frame 'heliographic_stonyhurst' the following metadata is missing: dsun_obs,hglt_obs,hgln_obs For frame 'heliographic_carrington' the following metadata is missing: dsun_obs,crln_obs,crlt_obs [sunpy.map.mapbase] WARNING: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer. For frame 'heliographic_stonyhurst' the following metadata is missing: dsun_obs,hglt_obs,hgln_obs For frame 'heliographic_carrington' the following metadata is missing: dsun_obs,crln_obs,crlt_obs [sunpy.map.mapbase]
In [5]:
aia_193_map_1025 = sunpy.map.Map("../../src/AIA/20221025/193/lvl15/aia.lev1_euv_12s.2022-10-25T192003Z.193.image.fits")
aia_193_map_1025_crop = aia_193_map_1025.submap(SkyCoord(-400*u.arcsec, 100*u.arcsec,frame=aia_193_map_1025.coordinate_frame),
top_right=SkyCoord(0*u.arcsec, 400*u.arcsec,frame=aia_193_map_1025.coordinate_frame))
aia_193_map_1025_crop.plot()
Out[5]:
<matplotlib.image.AxesImage at 0x7f52ff2e4f50>
In [6]:
fig = plt.figure(constrained_layout=True)
ax = fig.add_subplot(111, projection=aia_193_map_1025_crop)
aia_193_map_1025_crop.plot(axes=ax, title=False)
bounds = ax.axis()
with propagate_with_solar_surface(rotation_model='rigid'):
eis_east_region.to_pixel(aia_193_map_1025_crop.wcs).plot(ax=ax, edgecolor='#1B813E', lw=1.5,alpha=0.8)
eis_west_region.to_pixel(aia_193_map_1025_crop.wcs).plot(ax=ax, edgecolor='#1B813E', lw=1.5,alpha=0.8)
ax.axis(bounds)
Out[6]:
(-0.5, 666.5, -0.5, 500.5)
In [7]:
# based on mag2pfsspy https://github.com/STBadman/magnetograms2pfsspy/blob/master/mag2pfsspy.py#L12
# written by Samuel Badman, but updated to fit pfsspy 1.2.0
# pfsspy is archived at this moment, but I still use it when the succeeding sunkit-magex is under development
# def extract_br(m):
# br = m.data
# br = br - np.nanmean(br)
# # GONG maps have their LH edge at -180deg, so roll to get it at 0deg
# br = np.roll(br, int((m.meta['CRVAL1'] + 180)/np.abs(m.meta['cdelt1'])), axis=1)
# br = np.nan_to_num(br)
# return br*1e5 # Gauss to nT
def hmi2pfsspy(filepath,rss=2.5,nr=60,ret_magnetogram=False,resample=None, map_type='daily',
ext=0):
if map_type == 'daily':
hmi_fits_data = fits.getdata(filepath,ext=ext)
hmi_fits_header = fits.getheader(filepath,ext=ext)
hmi_fits_header['CUNIT2'] = 'Sine Latitude'
# for card in ['HGLN_OBS','CRDER1','CRDER2','CSYSER1','CSYSER2'] :
# hmi_fits_header[card] = 0
hmi_fits_header['T_OBS'] = hmi_fits_header['T_OBS'][:-8]+"_TAI"
hmi_fits_data = np.nan_to_num(hmi_fits_data, nan=0.0)
# since the cdelt1 is flipped, one should also update the crval1 either through sunpy.L0 (like sam) or manually from FITS header
# hmi_fits_header['CRVAL1'] = 120 + sun.L0(time=hmi_fits_header['T_OBS']).to_value(u.deg)
hmi_fits_header['CRVAL1'] = (2262*360*u.deg - \
(hmi_fits_header["CRVAL1"] + hmi_fits_header["NAXIS1"]*hmi_fits_header["CDELT1"])*u.deg).to_value(u.deg)
# print(hmi_fits_header['CRVAL1'])
hmi_map = sunpy.map.sources.HMISynopticMap(hmi_fits_data,hmi_fits_header)
elif map_type == 'cr':
hmi_fits_data = fits.getdata(filepath,ext=ext)
hmi_fits_header = fits.getheader(filepath,ext=ext)
hmi_map = sunpy.map.Map(hmi_fits_data, hmi_fits_header)
hmi_map = sunpy.map.Map(np.nan_to_num(hmi_map.data, nan=0.0), hmi_map.meta)
if resample is not None:
hmi_map_size = hmi_map.data.shape
x_superpixel = int(hmi_map_size[1]/resample[0].value)
y_superpixel = int(hmi_map_size[0]/resample[1].value)
x_std = x_superpixel
y_std = y_superpixel
kernel = Gaussian2DKernel(x_stddev=x_std, y_stddev=y_std)
hmi_map_conv = convolve(hmi_map.data, kernel, boundary='extend')
hmi_map_conv = sunpy.map.sources.HMISynopticMap(hmi_map_conv, hmi_fits_header)
hmi_map = hmi_map_conv.superpixel([x_superpixel, y_superpixel]*u.pix, func=np.nanmean)
# br_hmi = extract_br(hmi_map)
if ret_magnetogram:
return hmi_map
pfss_input = pfsspy.Input(hmi_map, nr, rss,)
pfss_output = pfsspy.pfss(pfss_input)
return pfss_output
In [8]:
hmi_daily_synoptic_map = hmi2pfsspy("../../src/HMI/CR2263/hmi.Mrdailysynframe_720s.20221025_120000_TAI.data.fits",
resample=[1440,720]*u.pix,ret_magnetogram=True, map_type='daily')
In [9]:
hmi_daily_synoptic_map.plot()
Out[9]:
<matplotlib.image.AxesImage at 0x7f52fe42a450>
In [10]:
hmi_cr_synoptic_map = hmi2pfsspy("../../src/HMI/CR2263/hmi.synoptic_mr_polfil_720s.2263.Mr_polfil.fits",
resample=[1440,720]*u.pix,ret_magnetogram=True, map_type='cr', ext=1)
In [11]:
hmi_cr_synoptic_map.plot()
INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]
WARNING: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer. For frame 'heliographic_stonyhurst' the following metadata is missing: dsun_obs,hglt_obs,hgln_obs For frame 'heliographic_carrington' the following metadata is missing: dsun_obs,crln_obs,crlt_obs [sunpy.map.mapbase]
Out[11]:
<matplotlib.image.AxesImage at 0x7f52fe495d00>
In [12]:
pfss_output = hmi2pfsspy("../../src/HMI/CR2263/hmi.Mrdailysynframe_720s.20221025_120000_TAI.data.fits",
resample=[720,360]*u.pix,nr=100, rss=2.5)
In [13]:
pfss_output_qsl = hmi2pfsspy("../../src/HMI/CR2263/hmi.Mrdailysynframe_720s.20221025_120000_TAI.data.fits",
resample=[360,180]*u.pix,nr=50, rss=2.5)
In [14]:
pfss_bg_qsl = pfss_output_qsl.bg.value
radius = np.exp(pfss_output_qsl.grid.rg[1:])
theta = np.arccos(pfss_output_qsl.grid.sg)
phi = pfss_output_qsl.grid.pg
theta = np.flip(theta)
br = pfss_bg_qsl[:,:,1:,2].transpose(2,1,0)
br = np.flip(br, axis=1)
bt = pfss_bg_qsl[:,:,1:,1].transpose(2,1,0)
bt = np.flip(bt, axis=1)
bp = pfss_bg_qsl[:,:,1:,0].transpose(2,1,0)
bp = np.flip(bp, axis=1)
with open("../../sav/PFSS/hmi_pfss_20221025.binary", 'wb') as file:
file.write(radius.tobytes('F'))
file.write(theta.tobytes('F'))
file.write(phi.tobytes('F'))
file.write(br.tobytes('F'))
file.write(bt.tobytes('F'))
file.write(bp.tobytes('F'))
/home/yjzhu/anaconda3/envs/sunpy/lib/python3.12/site-packages/pfsspy/output.py:95: UserWarning: Could not parse unit string "Mx/cm^2" as a valid FITS unit.
See https://fits.gsfc.nasa.gov/fits_standard.html for the FITS unit standards.
warnings.warn(f'Could not parse unit string "{unit_str}" as a valid FITS unit.\n'
In [15]:
with fits.open("../../src/HMI/20221025/qmap/hmi.pfss_synframe.20221025_120000_TAI.Br.fits") as hdul:
br_jsoc = hdul[1].data.copy()
radius_jsoc = np.array([1.000000, 1.008178, 1.016532, 1.025195, 1.034242,
1.043711, 1.053626, 1.064008, 1.074879, 1.086263,
1.098183, 1.110665, 1.123736, 1.137422, 1.151753,
1.166760, 1.182474, 1.198928, 1.216158, 1.234200,
1.253092, 1.272875, 1.293590, 1.315281, 1.337995,
1.361778, 1.386683, 1.412762, 1.440070, 1.468664,
1.498606, 1.529960, 1.562791, 1.597169, 1.633168,
1.670863, 1.710334, 1.751666, 1.794945, 1.840265,
1.887720, 1.937412, 1.989445, 2.043931, 2.100985,
2.160728, 2.223282, 2.288744, 2.357061, 2.427810,
2.500000])
In [16]:
fig, (ax1,ax2) = plt.subplots(2,1,layout='constrained')
ax1.imshow(br_jsoc[0,:,:], origin='lower', cmap='hmimag', vmin=-1500, vmax=1500)
ax2.imshow(br[0,:,:], cmap='hmimag', vmin=-1500, vmax=1500)
Out[16]:
<matplotlib.image.AxesImage at 0x7f52fdefe5a0>
In [17]:
def trace_flines(pfss_output, seed, max_steps='auto'):
tracer = tracing.FortranTracer(max_steps=max_steps)
with propagate_with_solar_surface():
flines = tracer.trace(seed, pfss_output)
return flines
In [18]:
fig = plt.figure(constrained_layout=True)
ax = fig.add_subplot(111, projection=aia_193_map_1025_crop)
aia_193_map_1025_crop.plot(axes=ax, title=False)
bounds = ax.axis()
with propagate_with_solar_surface(rotation_model='rigid'):
eis_east_region.to_pixel(aia_193_map_1025_crop.wcs).plot(ax=ax, edgecolor='#1B813E', lw=1.5,alpha=0.8)
eis_west_region.to_pixel(aia_193_map_1025_crop.wcs).plot(ax=ax, edgecolor='#1B813E', lw=1.5,alpha=0.8)
ax.axis(bounds)
seed_east_lon, seed_east_lat = np.meshgrid(np.linspace(-350,-270,9), np.linspace(150,240,11))
seed_east = SkyCoord(seed_east_lon.ravel()*u.arcsec,seed_east_lat.ravel()*u.arcsec,frame=aia_193_map_1025_crop.coordinate_frame)
with propagate_with_solar_surface(rotation_model='rigid'):
seed_east = seed_east[np.where(eis_east_region.contains(seed_east, aia_193_map_1025_crop.wcs))]
ax.plot_coord(seed_east,marker='.',color='red',markersize=5,alpha=0.5, lw=0)
seed_west_lon, seed_west_lat = np.meshgrid(np.linspace(-54.5,24.5,9), np.linspace(120.5,199.5,9))
seed_west = SkyCoord(seed_west_lon.ravel()*u.arcsec,seed_west_lat.ravel()*u.arcsec,frame=aia_193_map_1025_crop.coordinate_frame)
eis_west_region_circle = CircleSkyRegion(center=SkyCoord(-15*u.arcsec, 160*u.arcsec, frame=aia_193_map_1025_crop.coordinate_frame), radius=40*u.arcsec)
with propagate_with_solar_surface(rotation_model='rigid'):
eis_west_region_circle.to_pixel(aia_193_map_1025_crop.wcs).plot(ax=ax, edgecolor='#1B813E', lw=1.5,alpha=0.8)
seed_west = seed_west[np.where(eis_west_region_circle.contains(seed_west, aia_193_map_1025_crop.wcs))]
ax.plot_coord(seed_west,marker='.',color='red',markersize=5,alpha=0.5, lw=0)
Out[18]:
[<matplotlib.lines.Line2D at 0x7f52f9137530>]
In [19]:
# with propagate_with_solar_surface(rotation_model='rigid'):
east_flines = trace_flines(pfss_output, seed_east,)
west_flines = trace_flines(pfss_output, seed_west,)
/home/yjzhu/anaconda3/envs/sunpy/lib/python3.12/site-packages/pfsspy/output.py:95: UserWarning: Could not parse unit string "Mx/cm^2" as a valid FITS unit.
See https://fits.gsfc.nasa.gov/fits_standard.html for the FITS unit standards.
warnings.warn(f'Could not parse unit string "{unit_str}" as a valid FITS unit.\n'
/home/yjzhu/anaconda3/envs/sunpy/lib/python3.12/site-packages/pfsspy/tracing.py:180: UserWarning: At least one field line ran out of steps during tracing.
You should probably increase max_steps (currently set to auto) and try again.
warnings.warn(
In [20]:
def get_loop_length(line):
c = line.coords.cartesian.xyz
s = np.append(0., np.linalg.norm(np.diff(c.value, axis=1), axis=0).cumsum()) * c.unit
return np.diff(s).sum()
In [21]:
fig = plt.figure(constrained_layout=True)
ax = fig.add_subplot(111, projection=aia_193_map_1025_crop)
aia_193_map_1025_crop.plot(axes=ax, title=False)
bounds = ax.axis()
with propagate_with_solar_surface(rotation_model='rigid'):
eis_east_region.to_pixel(aia_193_map_1025_crop.wcs).plot(ax=ax, edgecolor='#1B813E', lw=1.5,alpha=0.8)
eis_west_region.to_pixel(aia_193_map_1025_crop.wcs).plot(ax=ax, edgecolor='#1B813E', lw=1.5,alpha=0.8)
ax.axis(bounds)
with propagate_with_solar_surface():
for fline, seed in zip(east_flines, seed_east):
if get_loop_length(fline) > 100 * u.Mm:
if fline.is_open:
ax.plot_coord(fline.coords, color='#89c3eb', linewidth=1,alpha=0.5)
ax.plot_coord(seed, marker='.', color='#89c3eb', markersize=5,alpha=0.5, lw=0)
else:
ax.plot_coord(fline.coords, color='red', linewidth=1,alpha=0.5)
ax.plot_coord(seed, marker='.', color='red', markersize=5,alpha=0.5, lw=0)
for fline, seed in zip(west_flines, seed_west):
if get_loop_length(fline) > 100 * u.Mm:
if fline.is_open:
ax.plot_coord(fline.coords, color='#89c3eb', linewidth=1,alpha=0.5)
ax.plot_coord(seed, marker='.', color='#89c3eb', markersize=5,alpha=0.5, lw=0)
else:
ax.plot_coord(fline.coords, color='red', linewidth=1,alpha=0.5)
ax.plot_coord(seed, marker='.', color='red', markersize=5,alpha=0.5, lw=0)
In [22]:
fig = plt.figure(constrained_layout=True)
ax = fig.add_subplot(111, projection=aia_193_map_1025)
aia_193_map_1025.plot(axes=ax, title=False)
bounds = ax.axis()
with propagate_with_solar_surface(rotation_model='rigid'):
eis_east_region.to_pixel(aia_193_map_1025.wcs).plot(ax=ax, edgecolor='#1B813E', lw=1.5,alpha=0.8)
eis_west_region.to_pixel(aia_193_map_1025.wcs).plot(ax=ax, edgecolor='#1B813E', lw=1.5,alpha=0.8)
ax.axis(bounds)
with propagate_with_solar_surface():
for fline, seed in zip(east_flines, seed_east):
if get_loop_length(fline) > 100 * u.Mm:
if fline.is_open:
ax.plot_coord(fline.coords, color='#89c3eb', linewidth=1,alpha=0.5)
ax.plot_coord(seed, marker='.', color='#89c3eb', markersize=5,alpha=0.5, lw=0)
else:
ax.plot_coord(fline.coords, color='red', linewidth=1,alpha=0.5)
ax.plot_coord(seed, marker='.', color='red', markersize=5,alpha=0.5, lw=0)
for fline, seed in zip(west_flines, seed_west):
if get_loop_length(fline) > 100 * u.Mm:
if fline.is_open:
ax.plot_coord(fline.coords, color='#89c3eb', linewidth=1,alpha=0.5)
ax.plot_coord(seed, marker='.', color='#89c3eb', markersize=5,alpha=0.5, lw=0)
else:
ax.plot_coord(fline.coords, color='red', linewidth=1,alpha=0.5)
ax.plot_coord(seed, marker='.', color='red', markersize=5,alpha=0.5, lw=0)
In [23]:
aia_193_map_1025_crop_large = aia_193_map_1025.submap(SkyCoord(-600*u.arcsec, -200*u.arcsec,frame=aia_193_map_1025.coordinate_frame),
top_right=SkyCoord(100*u.arcsec, 400*u.arcsec,frame=aia_193_map_1025.coordinate_frame))
fig = plt.figure(constrained_layout=True)
ax = fig.add_subplot(111, projection=aia_193_map_1025_crop_large)
aia_193_map_1025_crop_large.plot(axes=ax, title=False)
bounds = ax.axis()
with propagate_with_solar_surface(rotation_model='rigid'):
eis_east_region.to_pixel(aia_193_map_1025_crop_large.wcs).plot(ax=ax, edgecolor='#C1328E', lw=1.5,alpha=0.8)
eis_west_region.to_pixel(aia_193_map_1025_crop_large.wcs).plot(ax=ax, edgecolor='#1B813E', lw=1.5,alpha=0.8)
ax.axis(bounds)
with propagate_with_solar_surface():
for fline, seed in zip(east_flines, seed_east):
if get_loop_length(fline) > 100 * u.Mm:
if fline.is_open:
ax.plot_coord(fline.coords, color='#58B2DC', linewidth=1.5,alpha=0.5)
# ax.plot_coord(seed, marker='.', color='#89c3eb', markersize=5,alpha=0.5, lw=0)
else:
ax.plot_coord(fline.coords, color='#E83015', linewidth=1.5,alpha=0.5)
# ax.plot_coord(seed, marker='.', color='red', markersize=5,alpha=0.5, lw=0)
for fline, seed in zip(west_flines, seed_west):
if get_loop_length(fline) > 100 * u.Mm:
if fline.is_open:
ax.plot_coord(fline.coords, color='#58B2DC', linewidth=1.5,alpha=0.5)
# ax.plot_coord(seed, marker='.', color='#89c3eb', markersize=5,alpha=0.5, lw=0)
else:
ax.plot_coord(fline.coords, color='#E83015', linewidth=1.5,alpha=0.5)
# ax.plot_coord(seed, marker='.', color='red', markersize=5,alpha=0.5, lw=0)
region_ch_boundary = EllipseSkyRegion(center=SkyCoord(-440*u.arcsec, 70*u.arcsec, frame=aia_193_map_1025_crop_large.coordinate_frame),
width=200*u.arcsec, height=70*u.arcsec, angle=-50*u.deg)
region_ch_boundary.to_pixel(aia_193_map_1025_crop_large.wcs).plot(ax=ax, edgecolor='#BDC0BA', lw=1.5,alpha=0.8)
seeds_ch_boundary_lon, seeds_ch_boundary_lat = np.meshgrid(np.linspace(-570,-350,11), np.linspace(-10,160,13))
seeds_ch_boundary = SkyCoord(seeds_ch_boundary_lon.ravel()*u.arcsec,seeds_ch_boundary_lat.ravel()*u.arcsec,frame=aia_193_map_1025_crop_large.coordinate_frame)
with propagate_with_solar_surface(rotation_model='rigid'):
seeds_ch_boundary = seeds_ch_boundary[np.where(region_ch_boundary.contains(seeds_ch_boundary, aia_193_map_1025_crop_large.wcs))]
ax.plot_coord(seeds_ch_boundary,marker='.',color='#BDC0BA',markersize=5,alpha=0.5, lw=0)
Out[23]:
[<matplotlib.lines.Line2D at 0x7f52f833acc0>]
In [24]:
flines_ch_boundary = trace_flines(pfss_output, seeds_ch_boundary,max_steps=5000)
/home/yjzhu/anaconda3/envs/sunpy/lib/python3.12/site-packages/pfsspy/tracing.py:180: UserWarning: At least one field line ran out of steps during tracing. You should probably increase max_steps (currently set to 5000) and try again. warnings.warn(
In [25]:
aia_193_map_1025_crop_large = aia_193_map_1025.submap(SkyCoord(-600*u.arcsec, -200*u.arcsec,frame=aia_193_map_1025.coordinate_frame),
top_right=SkyCoord(100*u.arcsec, 400*u.arcsec,frame=aia_193_map_1025.coordinate_frame))
fig = plt.figure(constrained_layout=True)
ax = fig.add_subplot(111, projection=aia_193_map_1025_crop_large)
aia_193_map_1025_crop_large.plot(axes=ax, title=False)
bounds = ax.axis()
with propagate_with_solar_surface(rotation_model='rigid'):
eis_east_region.to_pixel(aia_193_map_1025_crop_large.wcs).plot(ax=ax, edgecolor='#C1328E', lw=1.5,alpha=0.8)
eis_west_region.to_pixel(aia_193_map_1025_crop_large.wcs).plot(ax=ax, edgecolor='#1B813E', lw=1.5,alpha=0.8)
ax.axis(bounds)
with propagate_with_solar_surface():
for ii, (fline, seed) in enumerate(zip(east_flines, seed_east)):
if get_loop_length(fline) > 100 * u.Mm:
if fline.is_open:
ax.plot_coord(fline.coords, color='#58B2DC', linewidth=1.5,alpha=0.5)
# ax.plot_coord(seed, marker='.', color='#89c3eb', markersize=5,alpha=0.5, lw=0)
else:
try:
vel_at_another_footpoint = sunpy.map.sample_at_coords(eis_195_velmap_shift_1025,
SkyCoord([fline.coords[0],]))
except:
vel_at_another_footpoint = np.nan
fline_coords_pixels = aia_193_map_1025_crop_large.wcs.world_to_pixel(fline.coords)
ax.text(fline_coords_pixels[0][0],fline_coords_pixels[1][0],str(ii),color='yellow',fontsize=10, transform=ax.transData,
va='center',ha='center')
if vel_at_another_footpoint > 0:
ax.plot_coord(fline.coords, color='#7BA23F', linewidth=1.5,alpha=0.5)
else:
ax.plot_coord(fline.coords, color='#E83015', linewidth=1.5,alpha=0.5)
# ax.plot_coord(seed, marker='.', color='red', markersize=5,alpha=0.5, lw=0)
for ii, (fline, seed) in enumerate(zip(west_flines, seed_west)):
if get_loop_length(fline) > 100 * u.Mm:
if fline.is_open:
ax.plot_coord(fline.coords, color='#58B2DC', linewidth=1.5,alpha=0.5)
# ax.plot_coord(seed, marker='.', color='#89c3eb', markersize=5,alpha=0.5, lw=0)
else:
ax.plot_coord(fline.coords, color='#E83015', linewidth=1.5,alpha=0.5)
# ax.plot_coord(seed, marker='.', color='red', markersize=5,alpha=0.5, lw=0)
for fline, seed in zip(flines_ch_boundary, seeds_ch_boundary):
if get_loop_length(fline) > 100 * u.Mm:
if fline.is_open:
ax.plot_coord(fline.coords, color='#58B2DC', linewidth=1.5,alpha=0.5)
# ax.plot_coord(seed, marker='.', color='#BDC0BA', markersize=5,alpha=0.5, lw=0)
else:
pass
# ax.plot_coord(fline.coords, color='#E83015', linewidth=1.5,alpha=0.5)
# ax.plot_coord(seed, marker='.', color='#BDC0BA', markersize=5,alpha=0.5, lw=0)
In [26]:
fig = plt.figure(constrained_layout=True)
ax = fig.add_subplot(111, projection=eis_195_velmap_shift_1025)
eis_195_velmap_shift_1025.plot(axes=ax, title=False,
norm=ImageNormalize(vmin=-20, vmax=20),
cmap='coolwarm',
aspect=eis_195_velmap_shift_1025.scale.axis2/eis_195_velmap_shift_1025.scale.axis1)
bounds = ax.axis()
with propagate_with_solar_surface(rotation_model='rigid'):
eis_east_region.to_pixel(eis_195_velmap_shift_1025.wcs).plot(ax=ax, edgecolor='#1B813E', lw=1.5,alpha=0.8)
eis_west_region.to_pixel(eis_195_velmap_shift_1025.wcs).plot(ax=ax, edgecolor='#1B813E', lw=1.5,alpha=0.8)
ax.axis(bounds)
with propagate_with_solar_surface(rotation_model='rigid'):
for ii, (fline, seed) in enumerate(zip(east_flines, seed_east)):
if get_loop_length(fline) > 100 * u.Mm:
if fline.is_open:
ax.plot_coord(fline.coords, color='#89c3eb', linewidth=1,alpha=0.5)
ax.plot_coord(seed, marker='.', color='#89c3eb', markersize=5,alpha=0.5, lw=0)
else:
ax.plot_coord(fline.coords, color='red', linewidth=1,alpha=0.5)
fline_coords_pixels = eis_195_velmap_shift_1025.wcs.world_to_pixel(fline.coords)
ax.plot_coord(seed, marker='.', color='red', markersize=5,alpha=0.5, lw=0)
footpoint_pixel = eis_195_velmap_shift_1025.wcs.world_to_pixel(fline.source_surface_footpoint)
# ax.plot_coord(fline.coords[:10], marker='.', color='red', markersize=5,alpha=0.5, lw=0)
# print(fline.source_surface_footpoint)
# print(fline.coords[0])
# print(fline.coords[-1])
# ax.text(fline_coords_pixels[0][0],fline_coords_pixels[1][0],str(ii),color='yellow',fontsize=10, transform=ax.transData)
# print(footpoint_pixel)
# for fline, seed in zip(west_flines, seed_west):
# if get_loop_length(fline) > 100 * u.Mm:
# if fline.is_open:
# ax.plot_coord(fline.coords, color='#89c3eb', linewidth=1,alpha=0.5)
# ax.plot_coord(seed, marker='.', color='#89c3eb', markersize=5,alpha=0.5, lw=0)
# else:
# ax.plot_coord(fline.coords, color='red', linewidth=1,alpha=0.5)
# ax.plot_coord(seed, marker='.', color='red', markersize=5,alpha=0.5, lw=0)
In [27]:
def read_qmap(filepath,layer=0):
with fits.open(filepath) as hdul:
data = hdul[1].data[layer,:,:-1]
header = hdul[1].header
header['T_OBS'] = header['T_OBS'][:-8]+"_TAI"
header["NAXIS1"] = header["NAXIS1"] - 1
old_crpix1 = header["CRPIX1"]
# for some unknown reason, sunpy cannot accept a heliogaographic WCS with CRPIX1 not at the center of the map
# so I have to manually shift the CRPIX1 and CRVAL1 to the center of the map
header["CRPIX1"] = (1 + header["NAXIS1"])/2
header["CRVAL1"] = header["CRVAL1"] + (old_crpix1 - header["CRPIX1"])*header["CDELT1"]
hmi_qmap_header_new = sunpy.map.make_fitswcs_header(data=data,
coordinate=SkyCoord(2262*360*u.deg - (header["CRVAL1"] + (header["NAXIS1"])*header["CDELT1"])*u.deg,
header["CRVAL2"]*u.deg,
rsun=696000000.*u.m,
frame="heliographic_carrington",
obstime=header["T_OBS"]),
reference_pixel=u.Quantity([header["CRPIX1"] - 1,header["CRPIX2"] - 1]*u.pix),
scale=u.Quantity([-header["CDELT1"],header["CDELT2"]]*u.deg/u.pix),
projection_code="CAR",
)
map_size = data.shape[::-1]
# cut the -90 and 90 degree pixel to fit into sunpy map (avoiding the overflow of the pixel boundary at the poles)
return sunpy.map.Map(data, hmi_qmap_header_new).submap([0, 1]*u.pix, top_right=[map_size[0], map_size[1]-2]*u.pix)
hmi_qmap_1025 = read_qmap("../../src/HMI/20221025/qmap/hmi.q_synframe.20221025_120000_TAI.slogQ.fits")
In [28]:
aia_193_map_1025_crop_large = aia_193_map_1025.submap(SkyCoord(-600*u.arcsec, -200*u.arcsec,frame=aia_193_map_1025.coordinate_frame),
top_right=SkyCoord(100*u.arcsec, 400*u.arcsec,frame=aia_193_map_1025.coordinate_frame))
with propagate_with_solar_surface(rotation_model='rigid'):
hmi_qmap_1025_repro_aia_crop_large = hmi_qmap_1025.reproject_to(aia_193_map_1025_crop_large.wcs, algorithm='adaptive')
with Image.open("../../figs/test_figs/20221025_flines.png") as img:
pyvista_img = np.asarray(img)
with rc_context(ms_style_dict):
fig = plt.figure(figsize=(4,15.5*4/7),
constrained_layout=True)
gs = fig.add_gridspec(3,1, height_ratios = [3.5,6,6], hspace=0.)
ax1 = fig.add_subplot(gs[0])
ax1.imshow(pyvista_img)
ax1.annotate(r'\textbf{East}', xy=(330, 510), xytext=(150, 700),
arrowprops={"arrowstyle":"-|>",
"connectionstyle":"arc3,rad=0.4","color":"#C1328E"},
fontsize=10, ha='center', va='center', color='#C1328E',
bbox=dict(boxstyle="round", fc="w", ec="#91989F", alpha=0.8),
)
ax1.annotate(r'\textbf{West}', xy=(650, 700), xytext=(880, 570),
arrowprops={"arrowstyle":"-|>",
"connectionstyle":"arc3,rad=0.4","color":"#1B813E"},
fontsize=10, ha='center', va='center', color='#1B813E',
bbox=dict(boxstyle="round", fc="w", ec="#91989F", alpha=0.8),
)
ax2 = fig.add_subplot(gs[1], projection=aia_193_map_1025_crop_large)
im2 = aia_193_map_1025_crop_large.plot(axes=ax2, title=False, cmap='sdoaia193',)
ax3 = fig.add_subplot(gs[2], projection=aia_193_map_1025_crop_large)
im3 = hmi_qmap_1025_repro_aia_crop_large.plot(axes=ax3, title=False, cmap='PuOr',
norm=ImageNormalize(vmin=-4,vmax=4))
ax1.tick_params(left=False, bottom=False, labelleft=False, labelbottom=False)
for ax_ in (ax2, ax3):
bounds = ax_.axis()
with propagate_with_solar_surface(rotation_model='rigid'):
eis_east_region.to_pixel(aia_193_map_1025_crop_large.wcs).plot(ax=ax_, edgecolor='#C1328E', lw=1.5,alpha=0.8)
eis_west_region.to_pixel(aia_193_map_1025_crop_large.wcs).plot(ax=ax_, edgecolor='#1B813E', lw=1.5,alpha=0.8)
ax_.axis(bounds)
with propagate_with_solar_surface():
for fline, seed in zip(east_flines, seed_east):
if get_loop_length(fline) > 100 * u.Mm:
if fline.is_open:
ax2.plot_coord(fline.coords, color='#58B2DC', linewidth=1.5,alpha=0.5)
ax3.plot_coord(fline.coords, color='#58B2DC', linewidth=1.5,alpha=0.5)
else:
ax2.plot_coord(fline.coords, color='#F05E1C', linewidth=1.5,alpha=0.5)
ax3.plot_coord(fline.coords, color='#F05E1C', linewidth=1.5,alpha=0.5)
for fline, seed in zip(west_flines, seed_west):
if get_loop_length(fline) > 100 * u.Mm:
if fline.is_open:
ax2.plot_coord(fline.coords, color='#58B2DC', linewidth=1.5,alpha=0.5)
ax3.plot_coord(fline.coords, color='#58B2DC', linewidth=1.5,alpha=0.5)
else:
ax2.plot_coord(fline.coords, color='#F05E1C', linewidth=1.5,alpha=0.5)
ax3.plot_coord(fline.coords, color='#F05E1C', linewidth=1.5,alpha=0.5)
for fline, seed in zip(flines_ch_boundary, seeds_ch_boundary):
if get_loop_length(fline) > 100 * u.Mm:
if fline.is_open:
ax2.plot_coord(fline.coords, color='#58B2DC', linewidth=1.5,alpha=0.5)
ax3.plot_coord(fline.coords, color='#58B2DC', linewidth=1.5,alpha=0.5)
else:
pass
ax2.coords[0].axislabels.set_visible(False)
clb, clb_ax = plot_colorbar(im3, ax3,
bbox_to_anchor=(0.05,1.05,0.9,0.05), orientation='horizontal',
fontsize=10,scilimits=(0,2),title=r'Signed Squashing Factor $\mathrm{slog} Q$')
clb_ax.tick_params(top=True, bottom=False, labeltop=True, labelbottom=False,
length=3)
clb_ax.xaxis.tick_top()
clb_ax.xaxis.set_label_position('top')
offset_text = clb_ax.xaxis.get_offset_text()
if offset_text is not None:
offset_text.set_visible(False)
ax1.text(0.03, 0.94, r"\textbf{a) HMI PFSS}", transform=ax1.transAxes, ha='left', va='top', fontsize=10,
color='white', path_effects=[path_effects.Stroke(linewidth=1.2, foreground='black'), path_effects.Normal()],
linespacing=1.8)
ax2.text(0.03, 0.96, r"\textbf{b) AIA 19.3\,nm }" + "\n" + \
r'\textbf{{{}}}\,'.format(aia_193_map_1025_crop_large.date.isot[5:-4]) + \
"\n" + r"\textbf{HMI PFSS}", transform=ax2.transAxes, ha='left', va='top', fontsize=10,
color='white', path_effects=[path_effects.Stroke(linewidth=1.2, foreground='black'), path_effects.Normal()],
linespacing=1.8)
ax3.text(0.03, 0.96, r"\textbf{c) HMI Daily} $\boldsymbol{Q}$\textbf{-map }" + "\n" + \
r'\textbf{{{}}}\,'.format(hmi_qmap_1025.date.isot[5:-4]),
transform=ax3.transAxes, ha='left', va='top', fontsize=10,
color='white', path_effects=[path_effects.Stroke(linewidth=1.2, foreground='black'), path_effects.Normal()],
linespacing=1.8)
fig.savefig("../../figs/ms_eis_eui_upflow/upflow_pfss.pdf",dpi=300,bbox_inches='tight')
plt.show()
WARNING: SunpyMetadataWarning: Missing metadata for observer: assuming Earth-based observer. For frame 'heliographic_stonyhurst' the following metadata is missing: dsun_obs,hglt_obs,hgln_obs For frame 'heliographic_carrington' the following metadata is missing: dsun_obs,crln_obs,crlt_obs [sunpy.map.mapbase]
In [33]:
with rc_context(ms_style_dict):
fig = plt.figure(figsize=(4,3),
constrained_layout=True)
ax = fig.add_subplot(111,)
ax.imshow(pyvista_img)
ax.annotate(r'\textbf{East}', xy=(330, 510), xytext=(150, 700),
arrowprops={"arrowstyle":"-|>",
"connectionstyle":"arc3,rad=0.4","color":"#C1328E"},
fontsize=10, ha='center', va='center', color='#C1328E',
bbox=dict(boxstyle="round", fc="w", ec="#91989F", alpha=0.8),
)
ax.annotate(r'\textbf{West}', xy=(650, 700), xytext=(880, 570),
arrowprops={"arrowstyle":"-|>",
"connectionstyle":"arc3,rad=0.4","color":"#1B813E"},
fontsize=10, ha='center', va='center', color='#1B813E',
bbox=dict(boxstyle="round", fc="w", ec="#91989F", alpha=0.8),
)
ax.axis("off")
fig.savefig("../../figs/test_figs/pfss_karbacher.pdf", dpi=300, bbox_inches='tight')
plt.show()
