#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy
from epygram import config
from epygram.base import FieldSet
# Constants
Cpd = config.Cpd
Rd = config.Rd
Rv = config.Rv
################
### FORMULAS ###
################
[docs]def hybridpressureAB2fluxpressure(A, B, Psurf):
"""
Computes the pressure at flux levels defined by hybrid-pressure coefficients.
- *A* = table of A coefficients \n
- *B* = table of B coefficients \n
- *Psurf* = surface pressure in Pa. \n
Returns a table of pressure.
"""
if not len(A) == len(B):
raise ValueError("A, B must have the same size.")
# to numpy
A = numpy.array(A)
B = numpy.array(B)
# computation
pi_tilde = A + B*Psurf
return pi_tilde
[docs]def hybridpressureAB2masspressure(A, B, Psurf):
"""
Computes the pressure at mass levels defined by hybrid-pressure coefficients.
- *A* = table of A coefficients \n
- *B* = table of B coefficients \n
- *Psurf* = surface pressure in Pa. \n
Returns a table of pressure.
"""
pi_tilde = hybridpressureAB2fluxpressure(A, B, Psurf)
pi = flux2masspressures(pi_tilde)
return pi
[docs]def flux2masspressures(pi_tilde, Ptop=config.default_Ptop):
"""
Converts pressures at flux levels to mass levels.
"""
L = len(pi_tilde)
if not isinstance(pi_tilde, numpy.ndarray):
pi_tilde = numpy.array(pi_tilde)
pi = numpy.zeros(L)
for k in range(1, L+1):
ik = k-1 # python arranging
if k == 1:
pi[ik] = (pi_tilde[ik] - Ptop) / (1. + Cpd/Rd)
else:
pi[ik] = numpy.sqrt(pi_tilde[ik]*pi_tilde[ik-1])
return pi
[docs]def mass2fluxpressures(pi, Ptop=config.default_Ptop):
"""
Converts pressures at mass levels to flux levels.
"""
L = len(pi)
if not isinstance(pi, numpy.ndarray):
pi = numpy.array(pi)
pi_tilde = numpy.zeros(L)
for k in range(1, L+1):
ik = k-1 # python arranging
if k == 1:
pi_tilde[ik] = (pi[ik] + Ptop) * (1. + Cpd/Rd)
else:
pi_tilde[ik] = pi[ik]**2 / pi_tilde[ik-1]
return pi_tilde
[docs]def pressure2altitude(R, T, pi=None, pi_tilde=None, Pdep=0., Phi_surf=0):
"""
Converts a pressure profile to mass-levels altitude (= height above ground if *Phi_surf == 0*).
The pressure can be given at mass levels (*pi*) or flux levels (*pi_tilde*),
or both (assuming they are consistent), with unit: Pa.
- *R* = table of R = specific gas constant (J/kg/K) \n
- *T* = table of Temperature (K) \n
- *Pdep* = table of Pressure departure: P-Pi, where P is the hydrostatic
pressure and Pi the Non-Hydrostatic pressure (Pa). \n
0. if Hydrostatic (default) \n
- *Phi_surf* = surface geopotential (J/kg)
"""
# get full pressures
if pi == None:
pi = flux2masspressures(pi_tilde)
elif pi_tilde == None:
pi_tilde = mass2fluxpressures(pi)
else:
if not isinstance(pi, numpy.ndarray):
pi = numpy.array(pi)
if not isinstance(pi_tilde, numpy.ndarray):
pi_tilde = numpy.array(pi_tilde)
L = len(pi)
if not len(pi) == len(pi_tilde) == len(R) == len(T):
raise ValueError("R, T, pi, pi_tilde must have the same size.")
if not isinstance(R, numpy.ndarray):
R = numpy.array(R)
if not isinstance(T, numpy.ndarray):
T = numpy.array(T)
if not isinstance(Pdep, numpy.ndarray):
Pdep = numpy.array(Pdep)
# Pre-computation
delta = numpy.zeros(L)
alpha = numpy.zeros(L)
for k in range(1, L+1):
ik = k-1 # python arranging
if k == 1:
delta[ik] = 1. + Cpd/Rd
alpha[ik] = 1.
else:
delta[ik] = (pi_tilde[ik] - pi_tilde[ik-1]) / pi[ik]
alpha[ik] = 1 - pi[ik]/pi_tilde[ik]
# Geopotential
Phi = numpy.zeros(L)
partialsum = numpy.zeros(L)
for k in reversed(range(1, L+1)):
ik = k-1 # python arranging
if k == L:
partialsum[ik] = 0.
else:
partialsum[ik] = partialsum[ik+1] + R[ik+1]*T[ik+1]*delta[ik+1]/(1.+Pdep[ik+1]/pi[ik+1])
Phi[ik] = Phi_surf + partialsum[ik] + R[ik]*T[ik]*alpha[ik]/(1.+Pdep[ik]/pi[ik])
# Altitude
z = Phi / config.g0
return z
[docs]def hybridpressureAB2altitude(A, B, R, T, Psurf, Pdep=0., Phi_surf=0, Ptop=config.default_Ptop):
"""
Computes the altitude of mass levels defined by hybrid-pressure coefficients.
(= height above ground if *Phi_surf == 0*).
- *A* = table of A coefficients \n
- *B* = table of B coefficients \n
- *Psurf* = surface pressure in Pa. \n
- *R* = table of R = specific gas constant (J/kg/K) \n
- *T* = table of Temperature (K) \n
- *Pdep* = table of Pressure departure: P-Pi, where P is the hydrostatic
pressure and Pi the Non-Hydrostatic pressure (Pa). \n
0. if Hydrostatic (default) \n
- *Phi_surf* = surface geopotential (J/kg)
"""
if not len(A) == len(B) == len(R) == len(T):
raise ValueError("A, B, R, T, must have the same size.")
pi_tilde = hybridpressureAB2fluxpressure(A, B, Psurf)
pi = flux2masspressures(pi_tilde, Ptop=Ptop)
z = pressure2altitude(R, T, pi=pi, pi_tilde=pi_tilde, Pdep=Pdep, Phi_surf=Phi_surf)
return z
[docs]def gfl2R(q, ql=0., qi=0., qr=0., qs=0., qg=0.):
"""
Computes air specific gas constant R according to specific humidity,
and hydrometeors if present.
"""
q = numpy.array(q)
ql = numpy.array(ql)
qi = numpy.array(qi)
qr = numpy.array(qr)
qs = numpy.array(qs)
qg = numpy.array(qg)
R = Rd + (Rv-Rd)*q - Rd*(ql+qi+qr+qs+qg)
return R
################
### GRAPHICS ###
################
[docs]def plotprofiles(profiles, labels=None, fidtype=None, Ycoordinate=None,
unit='SI', title=None, logscale=False, ema=False, zoom=None):
"""
To plot a series of profiles. Returns a :mod:`matplotlib` *Figure*.
Args: \n
- *profiles* being a :class:`epygram.base.FieldSet` of :class:`epygram.V1DField`,
or a single :class:`epygram.V1DField`. \n
All profiles are supposed to have the same unit,
and the same vertical coordinate.
- *labels* = a list of labels for the profiles (same length and same order).
- *fidtype* = type of fid for labelling the curve with *fid[fidtype]*;
if *None*, labels with raw fid.
- *Ycoordinate* = label for the Y coordinate.
- *unit* = label for X coordinate.
- *title* = title for the plot.
- *logscale* = to set Y logarithmic scale
- *ema* = to make emagram-like plots of Temperature
- *zoom*: a dict containing optional limits to zoom on the plot. \n
Syntax: e.g. {'ymax':500, ...}.
"""
import matplotlib.pyplot as plt
plt.rc('font', family='serif')
plt.rc('figure', autolayout=True)
colors = ['red', 'blue', 'green', 'orange', 'magenta', 'darkolivegreen', 'yellow', 'salmon', 'black']
linestyles = ['-', '--', '-.', ':']
if not isinstance(profiles, FieldSet):
p = profiles.deepcopy()
profiles = FieldSet()
profiles.append(p)
p0 = profiles[0]
if p0.geometry.coordinate in ('hybrid_pressure', 'pressure'):
reverseY = True
else:
reverseY = False
if p0.geometry.coordinate == 'hybrid_pressure':
Y = range(1, p0.geometry.grid['fieldlevels_number']+1)
if Ycoordinate == None:
if p0.geometry.coordinate == 'hybrid_pressure':
Ycoordinate = 'Level \nHybrid-Pressure \ncoordinate'
elif p0.geometry.coordinate == 'pressure':
Ycoordinate = 'Pressure (hPa)'
elif p0.geometry.coordinate == 'altitude':
Ycoordinate = 'Altitude (m)'
elif p0.geometry.coordinate == 'height':
Ycoordinate = 'Height (m)'
elif p0.geometry.coordinate == 'potential_vortex':
Ycoordinate = 'Potential \nvortex \n(PVU)'
else:
Ycoordinate = 'unknown \ncoordinate'
# Figure
fig, ax = plt.subplots(figsize=(6., 9.))
if logscale:
ax.set_yscale('log')
if reverseY:
plt.gca().invert_yaxis()
i = 0
for p in profiles:
if p.geometry.coordinate == 'pressure':
Y = numpy.array(p.geometry.grid['levels']) / 100
elif p.geometry.coordinate in ('altitude', 'height', 'potential_vortex'):
Y = numpy.array(p.geometry.grid['levels'])
if labels != None:
label = labels[i]
else:
if fidtype != None:
label = p.fid[fidtype]
else:
label = str(p.fid)
data = p.data
if ema:
mindata = numpy.inf
maxdata = -numpy.inf
alpha = 0.75
templines = numpy.arange(round(min(data),-1)-10, round(max(data),-1)+10+50, 10)
for t in templines:
plt.plot([t, t+(max(data)-min(data)) * alpha], [max(Y), min(Y)], color='grey', linestyle=':')
ax.set_yticks(numpy.linspace(0,1000, 11))
data = data + abs(Y-Y[-1]) * (max(data)-min(data))/(max(Y)-min(Y)) * alpha
unit = 'K'
mindata = min(mindata, min(data))
maxdata = max(maxdata, max(data))
plt.plot(data, Y, label=label,
color=colors[i%len(colors)], linestyle=linestyles[i//len(colors)])
i += 1
# Decoration
if reverseY:
ax.set_ylim(bottom=max(Y))
else:
ax.set_ylim(bottom=min(Y))
if zoom != None:
if 'ymin' in zoom.keys():
ax.set_ylim(bottom=zoom['ymin'])
if 'ymax' in zoom.keys():
ax.set_ylim(top=zoom['ymax'])
if title != None:
ax.set_title(title)
legend = ax.legend(loc='upper right', shadow=True)
for label in legend.get_texts():
label.set_fontsize('medium')
ax.set_xlabel(r'$'+unit+'$')
ax.set_ylabel(Ycoordinate)
if ema:
plt.grid(axis='y')
plt.xlim(mindata-10, maxdata+10)
else:
plt.grid()
return fig