Add example for reproducing the "Where is STEREO Today?" plot #6781
Conversation
c9c95d3 to
73218b7
Compare
|
Bleh, the axes label color isn't coming out correct on Matplotlib 3.7.0, but is correct on Matplotlib 3.6.1. I'll have to track that down. |
|
Got to admit, I do like the concept of a showcase. |
|
The example HMI Showcase: Cutout would also be a good candidate to move to this new "Showcase" category |
Lets move it. |
Since the connection with STEREO-B was lost so is it good to show the location of STEREO-B in the plot? |
The goal is to replicate the figure as close as possible and that includes STEREO-B despite that fact its lost to us. |
|
I'll pick this PR back up once matplotlib 3.7.1 is released, which should have my bug fix in it (matplotlib/matplotlib#25237) |
This comment was marked as outdated.
This comment was marked as outdated.
| Showcase | ||
| ======== | ||
|
|
||
| Examples that use an advanced combination of capabilities in sunpy |
There was a problem hiding this comment.
| Examples that use an advanced combination of capabilities in sunpy | |
| Examples that use an advanced combination of capabilities in sunpy. |
There was a problem hiding this comment.
Not actually a sentence...
examples/showcase/where_is_stereo.py
Outdated
| planets = ['Mercury', 'Venus', 'Earth', 'Mars'] | ||
| periods = [88, 225, 365, 687] | ||
| planet_coords = {planet: get_body_heliographic_stonyhurst(planet, | ||
| obstime + np.arange(period) * u.d) | ||
| for planet, period in zip(planets, periods)} | ||
|
|
||
| stereo_a = get_horizons_coord('STEREO-A', obstime) | ||
| stereo_b = get_horizons_coord('STEREO-B', obstime) | ||
|
|
||
| missions = ['Parker Solar Probe', 'Solar Orbiter', 'BepiColombo'] | ||
| mission_labels = {'Parker Solar Probe': 'PSP', 'Solar Orbiter': 'SO', 'BepiColombo': 'BEPICOLOMBO'} | ||
| periods = [100, 180, 120] # may need to be updated as orbits evolve | ||
| mission_coords = {mission: get_horizons_coord(mission, {'start': obstime, | ||
| 'stop': obstime + period * u.d, | ||
| 'step': '1d'}) | ||
| for mission, period in zip(missions, periods)} |
There was a problem hiding this comment.
| planets = ['Mercury', 'Venus', 'Earth', 'Mars'] | |
| periods = [88, 225, 365, 687] | |
| planet_coords = {planet: get_body_heliographic_stonyhurst(planet, | |
| obstime + np.arange(period) * u.d) | |
| for planet, period in zip(planets, periods)} | |
| stereo_a = get_horizons_coord('STEREO-A', obstime) | |
| stereo_b = get_horizons_coord('STEREO-B', obstime) | |
| missions = ['Parker Solar Probe', 'Solar Orbiter', 'BepiColombo'] | |
| mission_labels = {'Parker Solar Probe': 'PSP', 'Solar Orbiter': 'SO', 'BepiColombo': 'BEPICOLOMBO'} | |
| periods = [100, 180, 120] # may need to be updated as orbits evolve | |
| mission_coords = {mission: get_horizons_coord(mission, {'start': obstime, | |
| 'stop': obstime + period * u.d, | |
| 'step': '1d'}) | |
| for mission, period in zip(missions, periods)} | |
| # Times to sample each orbit at from today, chosen to span at least a full | |
| # orbit of every body plotted | |
| times = np.arange(1000) * u.day | |
| planets = ['Mercury', 'Venus', 'Earth', 'Mars'] | |
| planet_coords = {planet: get_body_heliographic_stonyhurst(planet, | |
| obstime + times) | |
| for planet in planets} | |
| stereo_a = get_horizons_coord('STEREO-A', obstime) | |
| stereo_b = get_horizons_coord('STEREO-B', obstime) | |
| missions = ['Parker Solar Probe', 'Solar Orbiter', 'BepiColombo'] | |
| mission_labels = {'Parker Solar Probe': 'PSP', 'Solar Orbiter': 'SO', 'BepiColombo': 'BEPICOLOMBO'} | |
| mission_coords = {mission: get_horizons_coord(mission, {'start': obstime, | |
| 'stop': obstime + times, | |
| 'step': '1d'}) | |
| for mission in missions} |
I got very confused about periods for a while here, so would suggest the above to simplify and clarify things.
There was a problem hiding this comment.
Hmm, I don't like the idea of plotting multiple orbits just to make the code simpler. It won't matter much for planets, but it could look bad for spacecraft.
An alternative approach that'd be more robust than hard-coding orbital periods, but would mean even more code, is to grab 1000 days as you have done here, but then figure out how many elements make up the first orbit (by looking for when the ecliptic longitude passes the starting point) and then plotting just that first orbit.
There was a problem hiding this comment.
Ah fair enough, I'm guessing it would look bad because of the dashed lines? In that case hardcoding is probably best, but maybe with a clarifying comment along the lines of # Orbital period in days above where periods is defined?
There was a problem hiding this comment.
Ah fair enough, I'm guessing it would look bad because of the dashed lines?
Not just that, but orbital precession means that successive orbits of spacecraft – most notably BepiColombo – do not lie on top of each other.
I've added a function to automatically determine the first full orbit in a trajectory.
| # Set Matplotlib settings to the desired appearance and initialize the axes. | ||
|
|
||
|
|
||
| mpl.rcParams.update({'figure.facecolor': 'black', |
There was a problem hiding this comment.
Oh god, why are we making it so ugly! 😆
Although we already closed #6301 with #6771, this PR adds an example that literally reproduces the "Where is STEREO Today?" plot as requested by #6301. Since this example is complex and particular, it is not intended to be a tutorial or a how-to guide. I have put it into a newly created "Showcase" category of examples, which I envision could house other complex and particular examples.