| 35 |
michaesp |
1 |
# -*- coding:utf-8 -*-
|
|
|
2 |
|
|
|
3 |
import numpy as np
|
|
|
4 |
from mpl_toolkits.basemap import Basemap
|
|
|
5 |
from matplotlib.collections import LineCollection
|
|
|
6 |
from matplotlib.colors import BoundaryNorm
|
|
|
7 |
from matplotlib.pyplot import get_cmap
|
|
|
8 |
|
|
|
9 |
|
|
|
10 |
'''
|
|
|
11 |
todo :
|
|
|
12 |
improve the plotting of trajectories crossing the 180th meridian
|
|
|
13 |
see how LAGRANTO does it
|
|
|
14 |
test with map centered on -180
|
|
|
15 |
'''
|
|
|
16 |
|
|
|
17 |
|
|
|
18 |
class Mapfigure(Basemap):
|
|
|
19 |
"""
|
|
|
20 |
Class based on Basemap with additional functionallity
|
|
|
21 |
such as plot_trajectories
|
|
|
22 |
"""
|
|
|
23 |
def __init__(self, resolution='i', projection='cyl',
|
|
|
24 |
domain=None, lon=None, lat=None, **kwargs):
|
|
|
25 |
|
|
|
26 |
if (domain is None) & (lon is not None):
|
|
|
27 |
domain = [lon.min(), lon.max(), lat.min(), lat.max()]
|
|
|
28 |
elif (domain is None) & (lon is None):
|
|
|
29 |
raise TypeError('lon, lat or domain need to be specified')
|
|
|
30 |
kwargs['llcrnrlon'] = domain[0]
|
|
|
31 |
kwargs['urcrnrlon'] = domain[1]
|
|
|
32 |
kwargs['llcrnrlat'] = domain[2]
|
|
|
33 |
kwargs['urcrnrlat'] = domain[3]
|
|
|
34 |
kwargs['resolution'] = resolution
|
|
|
35 |
kwargs['projection'] = projection
|
|
|
36 |
super(Mapfigure, self).__init__(**kwargs)
|
|
|
37 |
if (lon is not None):
|
|
|
38 |
self.x, self.y = self(lon, lat)
|
|
|
39 |
|
|
|
40 |
def drawmap(self, continent=False, nbrem=5, nbrep=5,
|
|
|
41 |
coastargs={}, countryargs={},
|
|
|
42 |
meridiansargs={}, parallelsargs={}):
|
|
|
43 |
"""
|
|
|
44 |
draw basic features on the map
|
|
|
45 |
nbrem: interval bewteen meridians
|
|
|
46 |
nbrep: interval between parallels
|
|
|
47 |
"""
|
|
|
48 |
self.drawcoastlines(**coastargs)
|
|
|
49 |
self.drawcountries(**countryargs)
|
|
|
50 |
merid = np.arange(0, 360, nbrem)
|
|
|
51 |
parall = np.arange(0, 180, nbrep)
|
|
|
52 |
self.drawmeridians(merid, labels=[0, 0, 0, 1], **meridiansargs)
|
|
|
53 |
self.drawparallels(parall, labels=[1, 0, 0, 0], **parallelsargs)
|
|
|
54 |
if continent:
|
|
|
55 |
self.fillcontinents(color='lightgrey')
|
|
|
56 |
|
|
|
57 |
def plot_traj(self, trajs, variable, cmap='Spectral', levels=None,
|
|
|
58 |
**kwargs):
|
|
|
59 |
segments, colors = list(zip(*
|
|
|
60 |
[(self._get_segments(traj),
|
|
|
61 |
traj[variable][:-1])
|
|
|
62 |
for traj in trajs])
|
|
|
63 |
)
|
|
|
64 |
segments = np.concatenate(segments)
|
|
|
65 |
colors = np.concatenate(colors)
|
|
|
66 |
cmap = get_cmap(cmap)
|
|
|
67 |
if levels is None:
|
|
|
68 |
levels = np.arange(0, np.nanmax(trajs[variable]))
|
|
|
69 |
norm = BoundaryNorm(levels, cmap.N)
|
|
|
70 |
lc = LineCollection(segments, array=colors, cmap=cmap, norm=norm,
|
|
|
71 |
**kwargs)
|
|
|
72 |
self.ax.add_collection(lc)
|
|
|
73 |
|
|
|
74 |
return lc
|
|
|
75 |
|
|
|
76 |
def _get_segments(self, trajs):
|
|
|
77 |
x, y = self(trajs['lon'], trajs['lat'])
|
|
|
78 |
points = np.array([x, y]).T.reshape(-1, 1, 2)
|
|
|
79 |
segments = np.concatenate([points[:-1], points[1:]], axis=1)
|
|
|
80 |
# remove all segments crossing the 180th meridian !! to be improved
|
|
|
81 |
diff = segments[:, 0, 0] - segments[:, 1, 0]
|
|
|
82 |
index = np.where((diff < 300) & (diff > -300))
|
|
|
83 |
return segments[index[0], :, :]
|
|
|
84 |
# return segments
|
|
|
85 |
|
|
|
86 |
|
|
|
87 |
class SkewT:
|
|
|
88 |
"""
|
|
|
89 |
Create a skewT-logP diagramm from
|
|
|
90 |
give useful function
|
|
|
91 |
"""
|
|
|
92 |
|
|
|
93 |
# Private attributes
|
|
|
94 |
SKEWNESS = 37.5
|
|
|
95 |
T_ZERO = 273.15
|
|
|
96 |
# P_bot is used to define the skewness of the plot
|
|
|
97 |
P_BOT = 100000
|
|
|
98 |
L = 2.501e6 # latent heat of vaporization
|
|
|
99 |
R = 287.04 # gas constant air
|
|
|
100 |
RV = 461.5 # gas constant vapor
|
|
|
101 |
EPS = R/RV
|
|
|
102 |
|
|
|
103 |
CP = 1005.
|
|
|
104 |
CV = 718.
|
|
|
105 |
KAPPA = (CP-CV)/CP
|
|
|
106 |
G = 9.81
|
|
|
107 |
|
|
|
108 |
# constants used to calculate moist adiabatic lapse rate
|
|
|
109 |
# See formula 3.16 in Rogers&Yau
|
|
|
110 |
A = 2./7.
|
|
|
111 |
B = EPS*L*L/(R*CP)
|
|
|
112 |
C = A*L/R
|
|
|
113 |
|
|
|
114 |
def __init__(self, ax, prange={'pbot': 1000., 'ptop': 100., 'dp': 1.}):
|
|
|
115 |
""" initalize a skewT instance """
|
|
|
116 |
|
|
|
117 |
self.pbot = prange['pbot']*100.
|
|
|
118 |
self.ptop = prange['ptop']*100.
|
|
|
119 |
self.dp = prange['dp']*100.
|
|
|
120 |
self.ax = ax
|
|
|
121 |
# Defines the ranges of the plot
|
|
|
122 |
self.plevs = np.arange(self.pbot, self.ptop-1, -self.dp)
|
|
|
123 |
self._isotherms()
|
|
|
124 |
self._isobars()
|
|
|
125 |
self._dry_adiabats()
|
|
|
126 |
self._moist_adiabats()
|
|
|
127 |
# self._mixing_ratio()
|
|
|
128 |
|
|
|
129 |
def _skewnessTerm(self, P):
|
|
|
130 |
return self.SKEWNESS * np.log(self.P_BOT/P)
|
|
|
131 |
|
|
|
132 |
def _isotherms(self):
|
|
|
133 |
for temp in np.arange(-100, 50, 10):
|
|
|
134 |
self.ax.semilogy(temp + self._skewnessTerm(self.plevs), self.plevs,
|
|
|
135 |
basey=np.e, color='blue',
|
|
|
136 |
linestyle=('solid' if temp == 0 else 'dashed'),
|
|
|
137 |
linewidth = .5)
|
|
|
138 |
|
|
|
139 |
def _isobars(self):
|
|
|
140 |
for n in np.arange(self.P_BOT, self.ptop-1, -10**4):
|
|
|
141 |
self.ax.plot([-40, 50], [n, n], color='black', linewidth=.5)
|
|
|
142 |
|
|
|
143 |
def _mixing_ratio(self):
|
|
|
144 |
rdv = 0.622
|
|
|
145 |
B1 = 243.04 # °C
|
|
|
146 |
C1 = 610.94 # Pa
|
|
|
147 |
A1 = 17.625
|
|
|
148 |
t = np.arange(-30, 50, 10)
|
|
|
149 |
m = np.zeros((self.plevs.size, t.size))
|
|
|
150 |
for i, temp in enumerate(t):
|
|
|
151 |
es = C1 * np.exp(A1*temp/(B1+temp))
|
|
|
152 |
m[:, i] = rdv*es/(self.plevs-es)
|
|
|
153 |
t, p = np.meshgrid(t, self.plevs)
|
|
|
154 |
self.ax.contour(t, p, m)
|
|
|
155 |
|
|
|
156 |
def _dry_adiabats(self):
|
|
|
157 |
for tk in self.T_ZERO+np.arange(-30, 210, 10):
|
|
|
158 |
dry_adiabat = tk * (self.plevs/self.P_BOT)**self.KAPPA - (
|
|
|
159 |
self.T_ZERO + self._skewnessTerm(self.plevs))
|
|
|
160 |
self.ax.semilogy(dry_adiabat, self.plevs, basey=np.e,
|
|
|
161 |
color='brown',
|
|
|
162 |
linestyle='dashed', linewidth=.5)
|
|
|
163 |
|
|
|
164 |
def _moist_adiabats(self):
|
|
|
165 |
ps = [p for p in self.plevs if p <= self.P_BOT]
|
|
|
166 |
tlevels = np.concatenate((np.arange(-40., 10.1, 5.),
|
|
|
167 |
np.arange(12.5, 45.1, 2.5)))
|
|
|
168 |
for temp in tlevels:
|
|
|
169 |
moist_adiabat = []
|
|
|
170 |
for p in ps:
|
|
|
171 |
temp -= self.dp*self.gamma_s(temp, p)
|
|
|
172 |
moist_adiabat.append(temp + self._skewnessTerm(p))
|
|
|
173 |
self.ax.semilogy(moist_adiabat, ps, basey=np.e, color='green',
|
|
|
174 |
linestyle='dashed', linewidth=.5)
|
|
|
175 |
|
|
|
176 |
def plot_data(self, p, T, color='black', style='solid'):
|
|
|
177 |
self.ax.semilogy(T + self._skewnessTerm(p*100), p*100,
|
|
|
178 |
basey=np.e, color=(color),
|
|
|
179 |
linestyle=(style), linewidth=1.5)
|
|
|
180 |
self.ax.axis([-40, 50, self.pbot, self.ptop])
|
|
|
181 |
self.ax.set_xlabel('Temperature ($^{\circ}$ C)')
|
|
|
182 |
xticks = np.arange(-40, 51, 5)
|
|
|
183 |
self.ax.set_xticks(xticks, ['' if tick % 10 != 0 else str(tick)
|
|
|
184 |
for tick in xticks])
|
|
|
185 |
self.ax.set_ylabel('Pressure (hPa)')
|
|
|
186 |
yticks = np.arange(self.pbot, self.ptop-1, -10**4)
|
|
|
187 |
self.ax.set_yticks(yticks)
|
|
|
188 |
self.ax.set_yticklabels(['{:2.0f}'.format(label)
|
|
|
189 |
for label in yticks/100])
|
|
|
190 |
|
|
|
191 |
def plot_windsbarbs(self, p, u, v, offset=40):
|
|
|
192 |
x = p.copy()
|
|
|
193 |
x[:] = offset
|
|
|
194 |
ax2 = self.ax.twinx()
|
|
|
195 |
ax2.barbs(x[::2], p[::2]*100, u[::2], v[::2])
|
|
|
196 |
ax2.set_yscale('log', basey=np.e)
|
|
|
197 |
yticks = np.arange(self.pbot, self.ptop-1, -10**4)
|
|
|
198 |
ax2.yaxis.set_ticks(yticks)
|
|
|
199 |
ax2.set_ylim([self.pbot, self.ptop])
|
|
|
200 |
ax2.set_xlim([-40, 50])
|
|
|
201 |
|
|
|
202 |
def es(self, T):
|
|
|
203 |
"""Returns saturation vapor pressure (Pascal) at temperature T (Celsius)
|
|
|
204 |
Formula 2.17 in Rogers&Yau"""
|
|
|
205 |
return 611.2*np.exp(17.67*T/(T+243.5))
|
|
|
206 |
|
|
|
207 |
def gamma_s(self, T, p):
|
|
|
208 |
"""Calculates moist adiabatic lapse rate for T (Celsius) and p (Pa)
|
|
|
209 |
Note: We calculate dT/dp, not dT/dz
|
|
|
210 |
See formula 3.16 in Rogers&Yau for dT/dz,
|
|
|
211 |
but this must be combined with
|
|
|
212 |
the dry adiabatic lapse rate (gamma = g/cp) and the
|
|
|
213 |
inverse of the hydrostatic equation (dz/dp = -RT/pg)"""
|
|
|
214 |
esat = self.es(T)
|
|
|
215 |
wsat = self.EPS*esat/(p-esat) # Rogers&Yau 2.18
|
|
|
216 |
numer = self.A*(T+self.T_ZERO) + self.C*wsat
|
|
|
217 |
denom = p * (1 + self.B*wsat/((T+self.T_ZERO)**2))
|
|
|
218 |
return numer/denom # Rogers&Yau 3.16
|