-
Notifications
You must be signed in to change notification settings - Fork 3
/
hmi_rot2d.py
112 lines (91 loc) · 3.59 KB
/
hmi_rot2d.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
#!/usr/bin/env python
import numpy as np
from matplotlib import pyplot as pl
from argparse import ArgumentParser
from matplotlib.projections import PolarAxes
from mpl_toolkits.axisartist import floating_axes
from mpl_toolkits.axisartist.grid_finder import (FixedLocator, MaxNLocator,
DictFormatter)
# data from SDO/HMI webpage
# http://jsoc.stanford.edu/HMI/Global_products.html
parser = ArgumentParser()
parser.add_argument('-o', '--output', type=str,
help="save plot to file")
parser.add_argument('--dpi', type=float,
help="resolution in DPI for saved image, passed to "
"rcParams['figure.dpi']")
parser.add_argument('--figsize', type=float, nargs=2,
help="figure size, passed to rcParams['figure.figsize']")
parser.add_argument('--cmap', type=str, default='viridis',
help="name of matplotlib colourmap to use (default='viridis')")
args = parser.parse_args()
if args.figsize:
pl.rcParams['figure.figsize'] = args.figsize
if args.dpi:
pl.rcParams['figure.dpi'] = args.dpi
try:
rot2d = np.load('data/hmi_rot2d.npy')
err2d = np.load('data/hmi_err2d.npy')
rmesh = np.load('data/hmi_rmesh.npy')
except IOError:
try:
from urllib2 import urlopen
except ImportError:
from urllib.request import urlopen
response = urlopen('http://jsoc.stanford.edu/SUM86/D917240671/S00000/rot.2d')
rot2d = np.loadtxt(response.readlines())
response.close()
np.save('data/hmi_rot2d.npy', rot2d)
response = urlopen('http://jsoc.stanford.edu/SUM86/D917240671/S00000/err.2d')
err2d = np.loadtxt(response.readlines())
response.close()
np.save('data/hmi_err2d.npy', err2d)
response = urlopen('http://jsoc.stanford.edu/SUM86/D917240671/S00000/rmesh.orig')
rmesh = np.loadtxt(response.readlines())[::4]
response.close()
np.save('data/hmi_rmesh.npy', rmesh)
# rot2d has 49 columns, latitudes are 90-i*15/8; i starts at 0
lat = np.array([15./8.*i for i in np.arange(49)])/180.*np.pi
r, th = np.meshgrid(rmesh, lat)
# adapted from
# http://matplotlib.org/examples/axes_grid/demo_floating_axes.html
fig = pl.gcf()
tr = PolarAxes.PolarTransform()
angle_ticks = [(0.0, r"$0^\circ$"),
(0.25*np.pi, r"$45^\circ$"),
(0.5*np.pi, r"$90^\circ$")]
grid_locator1 = FixedLocator([v for v, s in angle_ticks])
tick_formatter1 = DictFormatter(dict(angle_ticks))
grid_locator2 = MaxNLocator(2)
grid_helper = floating_axes.GridHelperCurveLinear(
tr, extremes=(0, 0.5*np.pi, 0, 1),
grid_locator1=grid_locator1,
grid_locator2=grid_locator2,
tick_formatter1=tick_formatter1,
tick_formatter2=None)
ax1 = floating_axes.FloatingSubplot(fig, 111, grid_helper=grid_helper)
fig.add_subplot(ax1)
ax1.axis["left"].set_axis_direction("bottom")
ax1.axis["right"].set_axis_direction("top")
ax1.axis["bottom"].set_visible(False)
ax1.axis["top"].set_axis_direction("bottom")
ax1.axis["top"].toggle(ticklabels=True, label=True)
ax1.axis["top"].major_ticklabels.set_axis_direction("top")
ax1.axis["top"].label.set_axis_direction("top")
ax1.axis["left"].label.set_text(r'$r/R_\odot$')
aux_ax = ax1.get_aux_axes(tr)
aux_ax.patch = ax1.patch
ax1.patch.zorder = 0.9
Ncontours = 20
data = rot2d.T[::-1]
data[err2d.T[::-1]/data>0.01] = np.nan
c = aux_ax.contourf(th, r, data, Ncontours, cmap=args.cmap)
pl.colorbar(c, label='rotation rate (nHz)')
# plot base of convection zone
th = np.linspace(0., np.pi/2, 101)
r = np.ones(len(th))*0.713
aux_ax.plot(th, r, 'k--')
if args.output:
pl.savefig(args.output)
else:
pl.show()