Subversion Repositories lagranto.arpege

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 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