-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathplt-bfly.py
135 lines (122 loc) · 4.43 KB
/
plt-bfly.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
# Script to generate a solar magnetic field butterfly diagram
# Import libraries
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import sunpy
import sunpy.io
import os.path
import datetime
# light_figs()
# Specify data directories and data range
cr0 = 1911
cr1 = 2269
mdidat = os.path.expanduser('~/data/mdi.Synoptic_Mr.polfil/')
mdicrs = np.arange(cr0, 2096)
hmidat = os.path.expanduser('~/data/hmi.Synoptic_Mr.polfil/')
hmicrs = np.arange(2097, cr1)
# Specify any parameters needed for the output map
latres = 1440
lat = np.linspace(-1,1,latres)
# Generate a time array
tarr = []
# Generate any output arrays
bfly = np.zeros([latres, len(mdicrs)+len(hmicrs)], dtype=np.double)
i = 0
# Loop through MDI files
for cr in mdicrs:
# Read and process a given synoptic chart into a profile
filenm = mdidat+'synop_Mr_0.polfil.'+str(cr).zfill(4)+'.fits'
if os.path.isfile(filenm):
d = sunpy.io.read_file(filenm)
tarr.append(d[0].header['T_ROT'])
crdat = d[0].data
crdat[np.where(abs(crdat) > 3000)] = 0
bfly[:,i] = np.interp(lat,np.linspace(-1,1,crdat.shape[0]),
crdat.sum(axis=1)/sum(np.isfinite(crdat[540,:])))
else:
bfly[:,i] = np.NaN
i += 1
# Loop through HMI files
for cr in hmicrs:
# Read and process a given synoptic chart into a profile
filenm = hmidat+'hmi.synoptic_mr_polfil_720s.'+str(cr).zfill(4)+'.Mr_polfil.fits'
if os.path.isfile(filenm):
d = sunpy.io.read_file(filenm)
tarr.append(d[1].header['T_ROT'])
crdat = d[1].data
bfly[:,i] = crdat.sum(axis=1)/crdat.shape[1]
else:
bfly[:,i] = np.NaN
i += 1
# Convert the timing array
tarr_cr = [ datetime.datetime.strptime((tarr[x])[0:16], "%Y.%m.%d_%H:%M")
for x in range(len(tarr)) ]
tmin, tmax = mdates.date2num([tarr_cr[0], tarr_cr[-1]])
# Plot the resulting diagram
f, (ax1) = plt.subplots(1, figsize=[12,6])
im = ax1.imshow(bfly, vmin=-10, vmax=10, extent=[tmin,tmax,-1,1],
aspect='auto', cmap='Greys_r')
ax1.set_yticks([-1,-0.5,0,0.5,1])
ax1a = ax1.twinx()
ax1b = ax1.twiny()
latticks = np.array([-90,-60,-30,30,60,90])
ax1a.set_yticks(np.sin(latticks*np.pi/180))
ax1a.set_yticklabels(latticks)
ax1b.set_xlim(cr0, cr1)
#ax1b.vlines([1911.5, 2104], ymin=-1, ymax=1, color='r', linestyles='dashed', alpha=0.5)
#ax1b.vlines([2096, 2252], ymin=-1, ymax=1, color='b', linestyles='dashed', alpha=0.5)
ax1.xaxis_date()
ax1.set_xlabel('Year')
ax1.set_ylabel('Sine latitude')
ax1a.set_ylabel('Latitude (degrees)')
ax1b.set_xlabel('Carrington rotation')
cb = plt.colorbar(im, ax=ax1a, label='Mean magnetic flux density (G)',
use_gridspec=True, fraction=0.05, pad=0.1, extend='both')
#f.tight_layout()
plt.savefig('bfly.pdf', dpi=300)
plt.savefig('bfly.png', dpi=300)
# Plot the resulting diagram vertically
f, (ax1) = plt.subplots(1, figsize=[5,7])
im = ax1.imshow(np.rot90(bfly), vmin=-10, vmax=10, extent=[-1,1,tmax,tmin],
aspect='auto', cmap='Greys_r')
ax1.set_xticks([-1,-0.5,0,0.5,1])
ax1a = ax1.twiny()
ax1b = ax1.twinx()
latticks = np.array([-90,-60,-30,30,60,90])
ax1a.set_xticks(np.sin(latticks*np.pi/180))
ax1a.set_xticklabels(latticks)
ax1b.set_ylim(cr1, cr0)
#ax1b.vlines([2193], ymin=-1, ymax=1)
ax1.yaxis_date()
ax1.set_ylabel('Year')
ax1.set_xlabel('Sine latitude')
ax1a.set_xlabel('Latitude (degrees)')
ax1b.set_ylabel('Carrington rotation')
cb = plt.colorbar(im, ax=ax1b, label='Mean magnetic flux density (G)',
use_gridspec=True, fraction=0.1, pad=0.1, extend='both',
location='bottom')
f.tight_layout()
plt.savefig('bfly_v.pdf', dpi=300)
plt.savefig('bfly_v.png', dpi=300)
# keynote_figs()
f, (ax1) = plt.subplots(1, figsize=[6,3.5])
im = ax1.imshow(bfly, vmin=-10, vmax=10, extent=[tmin,tmax,-1,1],
aspect='auto', cmap='Greys_r')
ax1.set_yticks([-1,-0.5,0,0.5,1])
ax1a = ax1.twinx()
ax1b = ax1.twiny()
latticks = np.array([-90,-60,-30,30,60,90])
ax1a.set_yticks(np.sin(latticks*np.pi/180))
ax1a.set_yticklabels(latticks)
ax1b.set_xlim(cr0, cr1)
ax1.xaxis_date()
ax1.set_xlabel('Year')
ax1.set_ylabel('Sine latitude')
ax1a.set_ylabel('Latitude (degrees)')
ax1b.set_xlabel('Carrington rotation')
cb = plt.colorbar(im, ax=ax1a, label='Mean magnetic flux density (G)',
use_gridspec=True, fraction=0.1, pad=0.15, extend='both')
#f.tight_layout()
plt.savefig('bfly_small.pdf', dpi=300)
plt.savefig('bfly_small.png', dpi=300)