Extracting objects from the plot

[MohamedNedal]

Hello, I would like to ask is it possible to extract the Parker spiral line object (dashed black curve) from the plot somehow, to plot it into another plot?
Also when I set the coord_sys to be Stonyhurst, the inner grid circles disappear. Is it possible to show them?
Here's the code to reproduce the same plot below:

from solarmach import SolarMACH
import warnings
warnings.filterwarnings('ignore')
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi'] = 100
plt.rcParams['figure.facecolor'] = 'white'
plt.rcParams['savefig.dpi'] = 300
plt.rcParams['savefig.format'] = 'png'
plt.rcParams['savefig.facecolor'] = 'white'
plt.rcParams['figure.figsize'] = [5,5]


# necessary options
body_list = ['STEREO-A', 'Earth', 'PSP']
vsw_list = [400, 400, 400]
obs_date = '2022-01-03 12:40:00'

# optional parameters
reference_long = -45
reference_lat = 12
reference_vsw = 400
long_offset = 270
plot_spirals = True
numbered_markers = True
return_plot_object = True
plot_sun_body_line = True
show_earth_centered_coord = False

# initialize
sm = SolarMACH(obs_date, body_list, vsw_list, 
               reference_long, reference_lat,
               coord_sys='Stonyhurst'
              )

# make plot
fig, ax = sm.plot(plot_spirals=plot_spirals, 
                  plot_sun_body_line=plot_sun_body_line, 
                  show_earth_centered_coord=show_earth_centered_coord, 
                  reference_vsw=reference_vsw, 
                  numbered_markers=numbered_markers, 
                  return_plot_object=return_plot_object,
                  long_offset=long_offset
                 )

Versions:
solarmach.__version__ is 0.1.5
matplotlib.__version__ is 3.7.1
python --version is 3.9.12

[jgieseler]

Hi!

I would like to ask is it possible to extract the Parker spiral line object (dashed black curve) from the plot somehow, to plot it into another plot?

Yes, but it is not straightforward. You can use the return_plot_object option to extract the whole matplotlib figure and axis (see https://solarmach.readthedocs.io/en/latest/usage.html#6.-Advanced:-edit-the-figure). Then it gets a bit dirty because you'll need to dig the data out from the matplotlib objects. In your above example, ax.get_children() contains all the different elements of the plot:

>>> ax.get_children()
[<matplotlib.lines.Line2D at 0x7f9992322760>,
 Text(-0.6101001930456195, 0.9644639851756787, '1'),
 <matplotlib.lines.Line2D at 0x7f99923228e0>,
 <matplotlib.lines.Line2D at 0x7f9992322ca0>,
 <matplotlib.lines.Line2D at 0x7f9992322e20>,
 Text(4.286249312733589e-08, 0.9817143423391952, '2'),
 <matplotlib.lines.Line2D at 0x7f999233b100>,
 <matplotlib.lines.Line2D at 0x7f999233b280>,
 <matplotlib.lines.Line2D at 0x7f999233b400>,
 Text(-2.4414444198563934, 0.7537240822989946, '3'),
 <matplotlib.lines.Line2D at 0x7f9992313f70>,
 <matplotlib.lines.Line2D at 0x7f999233b2b0>,
 <matplotlib.patches.FancyArrow at 0x7f9992313b50>,
 <matplotlib.lines.Line2D at 0x7f999233b610>,
 Text(18.3, -11.0, '1'),
 Text(18.3, -29.7, '2'),
 Text(18.3, -48.4, '3'),
 <matplotlib.legend.Legend at 0x7f999233b250>,
 <matplotlib.patches.Circle at 0x7f99922d1fa0>,
 Text(0.94, 0.16, 'Solar-MACH'),
 Text(0.94, 0.12, 'https://solar-mach.github.io'),
 <matplotlib.spines.Spine at 0x7f9997ddb7f0>,
 <matplotlib.spines.Spine at 0x7f9997ddb910>,
 <matplotlib.spines.Spine at 0x7f9997ddba30>,
 <matplotlib.spines.Spine at 0x7f9997ddbb50>,
 <matplotlib.projections.polar.ThetaAxis at 0x7f9997ddbcd0>,
 <matplotlib.projections.polar.RadialAxis at 0x7f9992313370>,
 Text(0.5, 1.0, '2022-01-03 12:40:00\n'),
 Text(0.0, 1.0, ''),
 Text(1.0, 1.0, ''),
 <matplotlib.legend.Legend at 0x7f99922d1580>,
 <matplotlib.patches.Wedge at 0x7f99923225e0>]

In this case, the dashed Parker line is the last child object of type matplotlib.lines.Line2D, so number 13. Then you can obtain the x and y data points of that line, and use them in another plot:

>>> x, y = ax.get_children()[13].get_data()
>>> fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
>>> ax.plot(x, y)
>>> ax.set_rmax(1.1)
>>> ax.set_theta_offset(np.deg2rad(270))
>>> plt.show()

The tricky part is to find which of the Line2D objects is the correct one. I would always start from the end because all the spacecraft Parker lines are plotted in the beginning. One quick check is to look at ax.get_children()[i].get_linestyle() for each i because for the dashed line this is '--', so a quick way to rule out the other Parker lines.

Also when I set the coord_sys to be Stonyhurst, the inner grid circles disappear. Is it possible to show them?

That is actually the oldest open bug we have and not connected to the coordinate system: #5 Unfortunately, no solution is in sight for it. 🫤

[MohamedNedal]

Thank you soo much!
I found this way makes it easier to get the index and name of each object:

for i, obj in enumerate(ax.get_children()):
    print(i, obj)

0 Line2D(STEREO A)
1 Annotation(-0.656177, 0.950801, '1')
2 Line2D(_child2)
3 Line2D(Earth)
4 Annotation(4.09706e-08, 0.990129, '2')
5 Line2D(_child5)
6 Line2D(Parker Solar Probe)
7 Annotation(-0.936361, 0.62178, '3')
8 Line2D(_child8)
9 FancyArrow()
10 Line2D(field line connecting to
ref. long. (vsw=400 km/s))
11 Annotation(1, 1, '1')
12 Annotation(1, 1, '2')
13 Annotation(1, 1, '3')
14 Legend
15 Circle(xy=(0, 0), radius=1.28349)
16 Text(0.94, 0.16, 'Solar-MACH')
17 Text(0.94, 0.12, 'https://solar-mach.github.io')
18 Spine
19 Spine
20 Spine
21 Spine
22 ThetaAxis(408.4333333333336,240.00000000000003)
23 RadialAxis(408.4333333333336,240.00000000000003)
24 Text(0.5, 1.0, '2021-10-28 15:15:00\n')
25 Text(0.0, 1.0, '')
26 Text(1.0, 1.0, '')
27 Legend
28 Wedge(center=(0.5, 0.5), r=0.5, theta1=270, theta2=630, width=0.5)

[MohamedNedal]

@jgieseler excuse me one more question. I'd like to know how to plot two (or more) reference arrows (black arrow) with their equivalent parker spirals

[jgieseler]

@jgieseler excuse me one more question. I'd like to know how to plot two (or more) reference arrows (black arrow) with their equivalent parker spirals

Sure. This needs a bit more code to be copied and re-run. For example, after executing your example code from the first post, you can add a second, red reference arrow at 95° with the following code:

import astropy.constants as aconst
import scipy.constants as const
import astropy.units as u

# defining new reference info
new_reference_long = 95
new_reference_lat = 5
new_reference_vsw = 350

# add anther arrow
arrow_dist = min([sm.max_dist/3.2, 2.])
ax.arrow(np.deg2rad(new_reference_long), 0.01, 0, arrow_dist, head_width=0.2, head_length=0.07, edgecolor='red', facecolor='red', lw=1.8, zorder=5, overhang=0.2)

#add another dashed Parker spiral at arrow
AU = const.au / 1000  # km
r_array = np.arange(0.007, (sm.max_dist+0.1)/np.cos(np.deg2rad(sm.max_dist_lat)) + 3.0, 0.001)
omega_ref = sm.solar_diff_rot(new_reference_lat)
alpha_ref = np.deg2rad(new_reference_long) + omega_ref / (reference_vsw / AU) * (sm.target_solar_radius*aconst.R_sun.to(u.AU).value - r_array) * np.cos(np.deg2rad(new_reference_lat))
ax.plot(alpha_ref, r_array * np.cos(np.deg2rad(new_reference_lat)), '--r', label=f'field line connecting to\nref. long. (vsw={new_reference_vsw} km/s)')

# reset max. radius for plot
ax.set_rmax(sm.max_dist + 0.3)

# save figure in current working directory
fig.savefig('solarmach.png', bbox_inches="tight")

But this is of course a very crude solution, as the function for calculating the Parker spiral is copied from the main code; i.e., if there are some changes to it, this example would be out of date.

[MohamedNedal]

Fantastic! Thank you so much :)

[MohamedNedal]

Hi @jgieseler , excuse me I've tried with other dates and times and to implement the solar wind speed for each body according to the recent update, but I found a little misalignment between the field line and the body location.
So I've changed reference_vsw to be new_reference_vsw in this line:

new_reference_vsw = solarmach.get_sw_speed('PSP', date, trange=1)

alpha_ref = np.deg2rad(new_reference_long)+omega_ref / (new_reference_vsw/AU)*(sm.target_solar_radius*aconst.R_sun.to(u.AU).value - r_array)*np.cos(np.deg2rad(new_reference_lat))

This gives a perfect alignment. Thank you!

[jgieseler]

I'm not sure if I completely follow you, as I also don't know what you try to achieve. What do you mean with "a little misalignment between the field line and the body location"? Keep in mind that (new_)reference_vsw defines the solar wind speed for the dashed Parker spiral that starts at the reference_long location at the Sun; it's not intended to connect to any body/spacecraft. And with solarmach.get_sw_speed('PSP', date, trange=1), you obtain the solar wind speed measured by PSP -- that might or might not be accurate to what you're looking for, depending on how far PSP is separated from the field line you're plotting.

[MohamedNedal]

I was trying to make a simple plot and build on it object by object to have more control on the figure's objects, and to use the solar wind speed at each body instead of the default value 400 km/s. When I changed the solar wind speed, I found this misalignment between the field line and the body, but it's corrected by using (new_reference_vsw/AU) instead of (reference_vsw/AU).
The problem of this approach is that if the magnetic footpoint longitude of the body is not available, it won't plot the field line, like in this example for the PSP. Here's the code:

import numpy as np
from solarmach import SolarMACH
import astropy.constants as aconst
import scipy.constants as const
import astropy.units as u

body_list = ['STEREO-A', 'Earth', 'PSP']
date = '2019-07-01 14:15:00'

vsw_sta = solarmach.get_sw_speed('STEREO-A', date, trange=1)
vsw_earth = solarmach.get_sw_speed('Earth', date, trange=1)
vsw_psp = solarmach.get_sw_speed('PSP', date, trange=1)

vsw_list = [vsw_sta, vsw_earth, vsw_psp]

sm = SolarMACH(date, body_list, coord_sys='Stonyhurst')

psp_lon = sm.coord_table[sm.coord_table['Spacecraft/Body']=='PSP']['Magnetic footpoint longitude (Stonyhurst)'].values[0]
psp_lat = sm.coord_table[sm.coord_table['Spacecraft/Body']=='PSP']['Stonyhurst latitude (°)'].values[0]

sta_lon = sm.coord_table[sm.coord_table['Spacecraft/Body']=='STEREO-A']['Magnetic footpoint longitude (Stonyhurst)'].values[0]
sta_lat = sm.coord_table[sm.coord_table['Spacecraft/Body']=='STEREO-A']['Stonyhurst latitude (°)'].values[0]

earth_lon = sm.coord_table[sm.coord_table['Spacecraft/Body']=='Earth']['Magnetic footpoint longitude (Stonyhurst)'].values[0]
earth_lat = sm.coord_table[sm.coord_table['Spacecraft/Body']=='Earth']['Stonyhurst latitude (°)'].values[0]

fig, ax = sm.plot(
                   plot_spirals=False,
                   plot_sun_body_line=True,
                   transparent=False,
                   numbered_markers=False,
                   long_offset=270,
                   return_plot_object=True
                )

# defining new reference info
new_reference_long = -2
new_reference_lat = 12
new_reference_vsw = 400

# add anther arrow
arrow_dist = min([sm.max_dist/3.2, 2.0])
ax.arrow(np.deg2rad(new_reference_long), 0.01, 0, arrow_dist, head_width=0.2, head_length=0.07,
         edgecolor='dodgerblue', facecolor='dodgerblue', lw=1.8, zorder=5, overhang=0.2)

# add another dashed Parker spiral
r_array = np.arange(0.007, (sm.max_dist+0.1)/np.cos(np.deg2rad(sm.max_dist_lat)) + 3.0, 0.001)
omega_ref = sm.solar_diff_rot(new_reference_lat)
alpha_ref = np.deg2rad(new_reference_long)+omega_ref / (new_reference_vsw/AU)*(sm.target_solar_radius*aconst.R_sun.to(u.AU).value - r_array)*np.cos(np.deg2rad(new_reference_lat))
ax.plot(alpha_ref, r_array*np.cos(np.deg2rad(new_reference_lat)), '--', color='dodgerblue',
        label=f'Field line connecting to\nref. long. (vsw={new_reference_vsw:.2f} km/s)')



new_reference_long = psp_lon
new_reference_lat = psp_lat
new_reference_vsw = solarmach.get_sw_speed('PSP', date, trange=1)
AU = const.au/1000  # km
r_array = np.arange(0.007, (sm.max_dist+0.1)/np.cos(np.deg2rad(sm.max_dist_lat)) + 3.0, 0.001)
omega_ref = sm.solar_diff_rot(new_reference_lat)
alpha_ref = np.deg2rad(new_reference_long)+omega_ref / (reference_vsw/AU)*(sm.target_solar_radius*aconst.R_sun.to(u.AU).value - r_array)*np.cos(np.deg2rad(new_reference_lat))
ax.plot(alpha_ref, r_array*np.cos(np.deg2rad(new_reference_lat)), '-', color='purple')



new_reference_long = sta_lon
new_reference_lat = sta_lat
new_reference_vsw = solarmach.get_sw_speed('STEREO-A', date, trange=1)
r_array = np.arange(0.007, (sm.max_dist+0.1)/np.cos(np.deg2rad(sm.max_dist_lat)) + 3.0, 0.001)
omega_ref = sm.solar_diff_rot(new_reference_lat)
alpha_ref = np.deg2rad(new_reference_long)+omega_ref / (reference_vsw/AU)*(sm.target_solar_radius*aconst.R_sun.to(u.AU).value - r_array)*np.cos(np.deg2rad(new_reference_lat))
ax.plot(alpha_ref, r_array*np.cos(np.deg2rad(new_reference_lat)), '-', color='red')



new_reference_long = earth_lon
new_reference_lat = earth_lat
new_reference_vsw = solarmach.get_sw_speed('Earth', date, trange=1)
r_array = np.arange(0.007, (sm.max_dist+0.1)/np.cos(np.deg2rad(sm.max_dist_lat)) + 3.0, 0.001)
omega_ref = sm.solar_diff_rot(new_reference_lat)
alpha_ref = np.deg2rad(new_reference_long)+omega_ref / (reference_vsw/AU)*(sm.target_solar_radius*aconst.R_sun.to(u.AU).value - r_array)*np.cos(np.deg2rad(new_reference_lat))
ax.plot(alpha_ref, r_array*np.cos(np.deg2rad(new_reference_lat)), '-', color='green')


# reset max. radius for plot
ax.set_rmax(sm.max_dist+0.2)

# show the new legend box
ax.legend(loc='center left', bbox_to_anchor=(1.1,0.5), handlelength=3, markerscale=0.7)
ax.grid(True, which='both')
fig

[jgieseler]

The problem is that you are mixing up the longitude of the spacecraft and the longitude of the magnetic footpoint at the Sun, that is, the footpoint of the magnetic field line connecting to that spacecraft.

For example, in the following figure STEREO-A has a longitude of 270°. But the longitude of the magnetic footpoint of the Parker spiral connecting STEREO-A with the Sun is at ~320°.

The code above is written so that it takes the magnetic footpoint longitude as the starting point (because you asked for another reference arrow with Parker spiral). You can't use it straightforward to connect to a given longitude in the interplanetary medium. For that, you have to start at the spacecraft and calculate back to the Sun.

[MohamedNedal]

I see. Thanks a lot