In [4]:
%matplotlib inline
from numpy import load, array, radians, sin, cos, linspace, mean, log, isnan, nan, nanmin, nanmax, nanmean, abs, zeros, exp, where,\
concatenate, diff
from numpy.ma import masked_array
from scipy.interpolate import interp1d
In [5]:
# Exception handling, with line number and stuff
import linecache
import sys
def PrintException():
exc_type, exc_obj, tb = sys.exc_info()
f = tb.tb_frame
lineno = tb.tb_lineno
filename = f.f_code.co_filename
linecache.checkcache(filename)
line = linecache.getline(filename, lineno, f.f_globals)
print 'EXCEPTION IN ({}, LINE {} "{}"): {}'.format(filename, lineno, line.strip(), exc_obj)
In [6]:
Rs_da=287.05 # Specific gas const for dry air, J kg^{-1} K^{-1}
Rs_v=461.51 # Specific gas const for water vapour, J kg^{-1} K^{-1}
Cp_da=1004.6 # Specific heat at constant pressure for dry air
Cv_da=719. # Specific heat at constant volume for dry air
Cp_v=1870. # Specific heat at constant pressure for water vapour
Cv_v=1410. # Specific heat at constant volume for water vapour
Cp_lw=4218 # Specific heat at constant pressure for liquid water
Epsilon=0.622 # Epsilon=R_s_da/R_s_v; The ratio of the gas constants
degCtoK=273.15 # Temperature offset between K and C (deg C)
rho_w=1000. # Liquid Water density kg m^{-3}
grav=9.81 # Gravity, m s^{-2}
Lv=2.5e6 # Latent Heat of vaporisation
boltzmann=5.67e-8 # Stefan-Boltzmann constant
mv=18.0153 # Mean molar mass of water vapor(g/mol)
In [7]:
import re
station_list_search='/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/igra-stations.txt'
station_metadata=[]
f = open(station_list_search,'r')
for line in f:
line = line.strip()
line=re.sub(r'([A-Z])\s([A-Z])', r'\1_\2',line)
line=re.sub(r'([A-Z])\s\s([A-Z])', r'\1_\2',line)
station_metadata.append(line.split())
f.close()
def station_info_search(stat):
for line in station_metadata:
if "%s" % stat in line:
st = line[2].lower().title().replace('_',' ')
lo = float(line[3])
la = float(line[4])
st_height = float(line[5])
return st,la,lo, st_height
In [8]:
# Parse the data
stat='43003'
station_name,la,lo, st_height=station_info_search(stat)
load_file = load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section_'
'IND_SOUNDING_INTERP_MEAN_Climat_19600501_20111001_relativedelta(days=+7)_%s.npz' % stat)
# % (date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), delta, stat))
variable_list={'pressures': 0, 'temps':1, 'dewpoints':2, 'winddirs':3, 'windspeeds':4, 'pot_temp':5,
'sat_vap_pres':6, 'vap_press':7, 'rel_hum':8, 'wvmr':9, 'sp_hum':10, 'sat_temp':11,
'theta_e':12, 'theta_e_sat':13, 'theta_e_minus_theta_e_sat':14, 'geop_height':15}
data=load_file['date_bin_mean_all_dates_one_station']
p=data[1,0,:]/100.
T=data[1,1,:]-273.15
Tk=data[1,1,:]
Td=T-data[1,2,:]
Tdk=Tk-data[1,2,:]
h=data[1,15,:]*10.
def u_v_winds(wind_direction, wind_speed):
wind_rad = radians(wind_direction)
u_wind=-((wind_speed)*sin(wind_rad))
v_wind=-((wind_speed)*cos(wind_rad))
return u_wind,v_wind
#print T
#print p
#print Td
u_wind,v_wind = u_v_winds(data[1,3,:], data[1,4,:])
In [9]:
def variable_name_index_match(variable, variable_list):
for key, value in variable_list.iteritems(): # iter on both keys and values
if key.startswith('%s' % variable) and key.endswith('%s' % variable):
arr_index_var=value
return arr_index_var
variable_list_line={'lcl_temp': 0, 'lcl_vpt':1, 'pbl_pressure':2, 'surface_pressure':3, 'T_eq_0':4}
In [1]:
# A few basic checks and possible conversions
def PressureinhPaCheck(pres):
try:
maxpres=max(pres)
except TypeError:
maxpres=pres
if maxpres>2000.:
print "WARNING: P>2000 hPa; did you input a value in Pa?"
# If necessary, convert to celsius
def TempInCelsiusCheck(tc):
if all(tc[~isnan(tc)]<100.):
tc += 273.15
return tc
# If necessary, convert to kelvin
def TempInKelvinCheck(tk):
if all(tk[~isnan(tk)]<100.):
tk += 273.15
return tk
# Routines to find enivoronmental temperature/dewpoint at pressure level, based on
# linear interpolation of sounding
def TempC500mb(tempc, pres):
"""Temperature in Celsius at 500 mb
INPUTS:
Temperature profile from sounding - tempc (C)
Pressure profile from sounding pres (hPa)
OUTPUTS: Temp(C)
Source: http://weather.uwyo.edu/upperair/indices.html
Prints a warning if a pressure value below 2000 Pa input, to ensure
that the units were input correctly.
"""
PressureinhPaCheck(pres)
tempc = TempInCelsiusCheck(tempc)
return interp_sounding(tempc,pres,500.)
def TempC850mb(tempc, pres):
"""Temperature in Celsius at 500 mb
INPUTS:
Temperature profile from sounding - tempc (C)
Pressure profile from sounding pres (hPa)
OUTPUTS: Temp(C)
Source: http://weather.uwyo.edu/upperair/indices.html
Prints a warning if a pressure value below 2000 Pa input, to ensure
that the units were input correctly.
"""
PressureinhPaCheck(pres)
tempc = TempInCelsiusCheck(tempc)
return interp_sounding(tempc,pres,850.)
def DewTempC850mb(dew_tempc, pres):
"""Temperature in Celsius at 500 mb
INPUTS:
Temperature profile from sounding - tempc (C)
Pressure profile from sounding pres (hPa)
OUTPUTS: Temp(C)
Source: http://weather.uwyo.edu/upperair/indices.html
Prints a warning if a pressure value below 2000 Pa input, to ensure
that the units were input correctly.
"""
PressureinhPaCheck(pres)
dew_tempc = TempInCelsiusCheck(dew_tempc)
return interp_sounding(dew_tempc,pres,850.)
def TempC700mb(tempc, pres):
"""Temperature in Celsius at 500 mb
INPUTS:
Temperature profile from sounding - tempc (C)
Pressure profile from sounding pres (hPa)
OUTPUTS: Temp(C)
Source: http://weather.uwyo.edu/upperair/indices.html
Prints a warning if a pressure value below 2000 Pa input, to ensure
that the units were input correctly.
"""
PressureinhPaCheck(pres)
tempc = TempInCelsiusCheck(tempc)
return interp_sounding(tempc,pres,700.)
def DewTempC700mb(dew_tempc, pres):
"""Temperature in Celsius at 500 mb
INPUTS:
Temperature profile from sounding - tempc (C)
Pressure profile from sounding pres (hPa)
OUTPUTS: Temp(C)
Source: http://weather.uwyo.edu/upperair/indices.html
Prints a warning if a pressure value below 2000 Pa input, to ensure
that the units were input correctly.
"""
PressureinhPaCheck(pres)
dew_tempc = TempInCelsiusCheck(dew_tempc)
return interp_sounding(dew_tempc,pres,700.)
def MeanFirst500m(vr, height, st_height):
"""Average of variable for first 500m above surface
INPUTS:
Sounding variable - vr
Heights of input sounding = height
Station height - st_height
OUTPUTS: Variable average for first 500m above surface
Source: http://weather.uwyo.edu/upperair/indices.html
Prints a warning if a pressure value below 2000 Pa input, to ensure
that the units were input correctly.
"""
# Calculate average for first 500m
fifty_m_above_surface=st_height+50.
fivehundred_m_above_surface=st_height+500.
y_points = linspace(fifty_m_above_surface, fivehundred_m_above_surface, 200) # Points for height interpolation in first 500m
#print y_points
#print vr
#print height
vr_500m = interp_sounding(vr,height,y_points)
mean_vr = nanmean(vr_500m)
#print fifty_m_above_surface
#print fivehundred_m_above_surface
return mean_vr
# Routines to lift parcel from various starting characteristics
def TCParcelLiftedFrom850To500(tempc, dew_tempc, pres):
"""Temperature in Celsius at 500 mb of a parcel lifted from 850 mb
INPUTS:
Temperature profile from sounding - tempc (C)
Dewpoint temperature profile from sounding - tempc (C)
Pressure profile from sounding pres (hPa)
OUTPUTS: Temp(C)
Source: http://weather.uwyo.edu/upperair/indices.html
Prints a warning if a pressure value below 2000 Pa input, to ensure
that the units were input correctly.
"""
try:
maxpres=max(pres)
except TypeError:
maxpres=pres
if maxpres>2000:
print "WARNING: P>2000 hPa; did you input a value in Pa?"
tempc = TempInCelsiusCheck(tempc)
dew_tempc = TempInCelsiusCheck(dew_tempc)
t850c = interp_sounding(tempc,pres,850.)
td850c = interp_sounding(dew_tempc,pres,850.)
#print t850c
#print td850c
parcel_profile = ParcelAscentDryToLCLThenMoistC(850.,t850c,td850c, pres)
t500c_lift_from_850 = interp_sounding(parcel_profile,pres,500.)
return t500c_lift_from_850
def TCParcelLiftedFromFirst500mTo500(tempc, dew_tempc, pres, heights, st_height):
"""Temperature in Celsius at 500 mb of a parcel lifted from 850 mb
INPUTS:
Temperature profile from sounding - tempc (C)
Dewpoint temperature profile from sounding - tempc (C)
Pressure profile from sounding pres (hPa)
OUTPUTS: Temp(C)
Source: http://weather.uwyo.edu/upperair/indices.html
Prints a warning if a pressure value below 2000 Pa input, to ensure
that the units were input correctly.
"""
try:
maxpres=max(pres)
except TypeError:
maxpres=pres
if maxpres>2000:
print "WARNING: P>2000 hPa; did you input a value in Pa?"
tempc = TempInCelsiusCheck(tempc)
dew_tempc= TempInCelsiusCheck(dew_tempc)
tc_first_500m = MeanFirst500m(tempc, heights, st_height)
tdc_first_500m = MeanFirst500m(dew_tempc, heights, st_height)
p_first_500m = MeanFirst500m(pres, heights, st_height)
parcel_profile = ParcelAscentDryToLCLThenMoistC(p_first_500m,tc_first_500m,tdc_first_500m, pres)
t500c_lift_from_first_500m = interp_sounding(parcel_profile,pres,500.)
print t500c_lift_from_first_500m
return t500c_lift_from_first_500m
def LiftDry(startp,startt,startdp,y_points):
"""Lift a parcel to discover certain properties.
INPUTS:
startp: Pressure (hPa)
startt: Temperature (C)
startdp: Dew Point Temperature (C)
"""
assert startt>startdp,"Not a valid parcel. Check Td<Tc %s %s" % (startt, startdp)
#Pres=linspace(startp,100,100)
if startt>100.:
startt=startt-273.15
# Lift the dry parcel
T_dry=(startt+273.15)*(y_points/startp)**(Rs_da/Cp_da)-273.15
return T_dry
# Routines from https://github.com/tchubb/SkewT/blob/master/build/lib.linux-x86_64-2.7/skewt/thermodynamics.py
def LiftWet(startt,pres):
#--------------------------------------------------------------------
# Lift a parcel moist adiabatically from startp to endp.
# Init temp is startt in C, pressure levels are in hPa
#--------------------------------------------------------------------
if startt>100.:
startt=startt-273.15
if pres[0]<pres[-1]:
pres = pres[::-1]
temp=startt
t_out=zeros(pres.shape);t_out[0]=startt
for ii in range(pres.shape[0]-1):
delp=pres[ii]-pres[ii+1]
temp=temp-100*delp*GammaW(temp+273.15,(pres[ii]-delp/2)*100)
t_out[ii+1]=temp
return t_out[::-1]
def Theta(tempk,pres,pref=100000.):
"""Potential Temperature
INPUTS:
tempk (K)
pres (Pa)
pref:
OUTPUTS: Theta (K)
Source: Wikipedia
Prints a warning if a pressure value below 2000 Pa input, to ensure
that the units were input correctly.
"""
return tempk*(pref/pres)**(Rs_da/Cp_da)
def TempK(theta,pres,pref=100000.):
"""Inverts Theta function."""
try:
minpres=min(pres)
except TypeError:
minpres=pres
if minpres<2000.:
print "WARNING: P<2000 Pa; did you input a value in hPa?"
return theta*(pres/pref)**(Rs_da/Cp_da)
def ThetaE():
"""Equivalent potential temperature"""
raise NotImplementedError
def ThetaV(tempk,pres,e):
"""Virtual Potential Temperature
INPUTS
tempk (K)
pres (Pa)
e: Water vapour pressure (Pa) (Optional)
"""
mixr=MixRatio(e,pres)
theta=Theta(tempk,pres)
return theta*(1.+mixr/Epsilon)/(1.+mixr)
def GammaW(tempk,pres,e=None):
"""Function to calculate the moist adiabatic lapse rate (deg C/Pa) based
on the temperature, pressure, and rh of the environment.
INPUTS:
tempk (K)
pres (Pa)
RH (%)
RETURNS:
GammaW: The moist adiabatic lapse rate (Dec C/Pa)
"""
tempc=tempk-degCtoK
es=SatVap(tempc)
ws=MixRatio(es,pres)
if e is None:
# assume saturated
e=es
w=MixRatio(e,pres)
tempv=VirtualTempFromMixR(tempk,w)
latent=Latentc(tempc)
A=1.0+latent*ws/(Rs_da*tempk)
B=1.0+Epsilon*latent*latent*ws/(Cp_da*Rs_da*tempk*tempk)
Rho=pres/(Rs_da*tempv)
Gamma=(A/B)/(Cp_da*Rho)
return Gamma
def Density(tempk,pres,mixr):
"""Density of moist air"""
virtualT=VirtualTempFromMixR(tempk,mixr)
return pres/(Rs_da*virtualT)
def VirtualTemp(tempk,pres,e):
"""Virtual Temperature
INPUTS:
tempk: Temperature (K)
e: vapour pressure (Pa)
p: static pressure (Pa)
OUTPUTS:
tempv: Virtual temperature (K)
SOURCE: hmmmm (Wikipedia)."""
tempvk=tempk/(1-(e/pres)*(1-Epsilon))
return tempvk
def VirtualTempFromMixR(tempk,mixr):
"""Virtual Temperature
INPUTS:
tempk: Temperature (K)
mixr: Mixing Ratio (kg/kg)
OUTPUTS:
tempv: Virtual temperature (K)
SOURCE: hmmmm (Wikipedia). This is an approximation
based on a m
"""
return tempk*(1.0+0.6*mixr)
def Latentc(tempc):
"""Latent heat of condensation (vapourisation)
INPUTS:
tempc (C)
OUTPUTS:
L_w (J/kg)
SOURCE:
http://en.wikipedia.org/wiki/Latent_heat#Latent_heat_for_condensation_of_water
"""
return 1000.*(2500.8 - 2.36*tempc + 0.0016*tempc**2. - 0.00006*tempc**3.)
def SatVap(tempc,phase="liquid"):
"""Calculate saturation vapour pressure over liquid water and/or ice.
INPUTS:
tempc: (C)
phase: ['liquid'],'ice'. If 'liquid', do simple dew point. If 'ice',
return saturation vapour pressure as follows:
Tc>=0: es = es_liquid
Tc <0: es = es_ice
RETURNS: e_sat (Pa)
SOURCE: http://cires.colorado.edu/~voemel/vp.html (#2:
CIMO guide (WMO 2008), modified to return values in Pa)
This formulation is chosen because of its appealing simplicity,
but it performs very well with respect to the reference forms
at temperatures above -40 C. At some point I'll implement Goff-Gratch
(from the same resource).
"""
over_liquid=6.112*exp(17.67*tempc/(tempc+243.12))*100.
over_ice=6.112*exp(22.46*tempc/(tempc+272.62))*100.
# return where(tempc<0,over_ice,over_liquid)
if phase=="liquid":
# return 6.112*exp(17.67*tempc/(tempc+243.12))*100.
return over_liquid
elif phase=="ice":
# return 6.112*exp(22.46*tempc/(tempc+272.62))*100.
return where(tempc<0,over_ice,over_liquid)
else:
raise NotImplementedError
def MixRatio(e,p):
"""Mixing ratio of water vapour
INPUTS
e (Pa) Water vapor pressure
p (Pa) Ambient pressure
RETURNS
qv (kg kg^-1) Water vapor mixing ratio`
"""
return Epsilon*e/(p-e)
def MixR2VaporPress(qv,p):
"""Return Vapor Pressure given Mixing Ratio and Pressure
INPUTS
qv (kg kg^-1) Water vapor mixing ratio`
p (Pa) Ambient pressure
RETURNS
e (Pa) Water vapor pressure
"""
return qv*p/(Epsilon+qv)
def VaporPressure(dwpt):
"""Water vapor pressure
INPUTS
dwpt (C) Dew Point Temperature (for SATURATION vapor
pressure use tempc)
RETURNS
e (Pa) Water Vapor Pressure
SOURCE:
Bolton, Monthly Weather Review, 1980, p 1047, eq. (10)
"""
return 611.2*exp(17.67*dwpt/(243.5+dwpt))
def DewPoint(e):
""" Use Bolton's (1980, MWR, p1047) formulae to find tdew.
INPUTS:
e (Pa) Water Vapor Pressure
OUTPUTS:
Td (C)
"""
ln_ratio=log(e/611.2)
Td=((17.67-ln_ratio)*degCtoK+243.5*ln_ratio)/(17.67-ln_ratio)
return Td-degCtoK
In [2]:
# Several indices
# http://weather.uwyo.edu/upperair/indices.html
def KIndex(tempc, dew_tempc, pres):
T850c = TempC850mb(tempc, pres)
T500c = TempC500mb(tempc, pres)
TD850c = DewTempC850mb(dew_tempc,pres)
T700c = TempC700mb(tempc,pres)
TD700c = DewTempC700mb(dew_tempc,pres)
return (T850c - T500c) + TD850c - (T700c - TD700c)
def CrossTotalsIndex(tempc,dew_tempc, pres):
TD850c = DewTempC850mb(dew_tempc,pres)
T500c = TempC500mb(tempc, pres)
return TD850c - T500c
def VerticalTotalsIndex(tempc,pres):
T850c = TempC850mb(tempc,pres)
T500c = TempC500mb(tempc, pres)
return T850c - T500c
def TotalTotalsIndex(tempc, dew_tempc, pres):
T850c = TempC850mb(tempc,pres)
T500c = TempC500mb(tempc, pres)
TD850c = DewTempC850mb(tempc,pres)
return (T850c - T500c) + (TD850c - T500c)
def LiftedIndex(tempc, dew_tempc, pres, heights, st_height):
"""LIFT = T500 - Tparcel
T500 = temperature in Celsius of the environment at 500 mb
Tparcel = 500 mb temperature in Celsius of a lifted parcel with
the average pressure, temperature, and dewpoint of the layer 500 m
above the surface
INPUTS:
500mb temp of parcel lifted from 850mb (C)
500mb temp (C)
pref:
OUTPUTS: Temp(C)
Source: http://glossary.ametsoc.org/wiki/Stability_index
http://weather.uwyo.edu/upperair/indices.html
Prints a warning if a pressure value below 2000 Pa input, to ensure
that the units were input correctly.
"""
T500c = TempC500mb(tempc, pres)
t500c_lift_from_first_500m = TCParcelLiftedFromFirst500mTo500(tempc, dew_tempc, pres, heights, st_height)
print t500c_lift_from_first_500m
print T500c
lift = T500c - t500c_lift_from_first_500m
return lift
def ShowalterIndex(tempc, dew_tempc, pres):
"""SHOW = T500 - Tparcel
T500 = Temperature in Celsius at 500 mb
Tparcel = Temperature in Celsius at 500 mb of a parcel lifted from 850 mb
INPUTS:
500mb temp of parcel lifted from 850mb (C)
500mb temp (C)
pref:
OUTPUTS: Temp(C)
Source: http://glossary.ametsoc.org/wiki/Stability_index
Prints a warning if a pressure value below 2000 Pa input, to ensure
that the units were input correctly.
"""
t500c = TempC500mb(tempc,pres)
t500c_lift_from_850 = TCParcelLiftedFrom850To500(tempc, dew_tempc, pres)
show = t500c - t500c_lift_from_850
return show
In [13]:
KIndex(T,Tdk,p)
Out[13]:
In [11]:
#Interpolate Sounding
#y_points=linspace(5000, 100000, 200) # Points for pressure interpolation
def interp_sounding(variable, pressures_s,y_points):
#print y_points
#print variable
#nan_mask = masked_array(array(variable, dtype=float), isnan(array(variable, dtype=float)))
#nan_mask_p = masked_array(array(pressures_s, dtype=float), isnan(array(variable, dtype=float)))
variable = [x for (y,x) in sorted(zip(pressures_s, variable), key=lambda pair: pair[0])]
pressures_s = [y for (y,x) in sorted(zip(pressures_s, variable), key=lambda pair: pair[0])]
interp = interp1d(pressures_s, variable, bounds_error=False, fill_value=nan)
y_interp = interp(y_points)
return y_interp
In [12]:
def LiftedCondensationLevelTemp(init_temp_k, dew_init_temp_k):
if (init_temp_k<100.):
init_temp_k = init_temp_k +273.15
if (dew_init_temp_k<100.):
dew_init_temp_k = dew_init_temp_k +273.15
return (1./(1./(dew_init_temp_k-56.) + log(init_temp_k/dew_init_temp_k)/800.)) + 56.
In [435]:
def LiftedCondensationLevelPres(mean_pres, lcl_temp, mean_temp_k, kappa):
return mean_pres * ((lcl_temp/mean_temp_k)**(1/kappa))
In [32]:
# This needs trying
def LCLMethod2(temp_k, dew_temp_k, pres, height, st_height):
temp_k = TempInKelvinCheck(temp_k)
dew_temp_k = TempInKelvinCheck(dew_temp_k)
dew_temp_c = dew_temp_k-273.15
dew_mean_temp_c = MeanFirst500m(dew_temp_c, height, st_height)
mean_temp_k = MeanFirst500m(temp_k, height, st_height)
vap_press = VaporPressure(dew_mean_temp_c)
mean_pres = MeanFirst500m(pres, height, st_height)
wvmr = MixRatio(vap_press,mean_pres*100.)
sat_temp_k = 55.+2840./(3.5*log(mean_temp_k)-log(vap_press/100.)-4.805)
lcl = mean_pres*(sat_temp_k/mean_temp_k)**((Cp_da+Cp_v*wvmr)/(Rs_da*(1+wvmr/Epsilon)))
return lcl
In [33]:
LCLMethod2(Tk, Tdk, p, h, st_height)
Out[33]:
In [27]:
p
Out[27]:
LCLT Temperature (K) at the LCL, the lifting condensation level, from an average of the lowest 500 meters. LCLT = [1 / ( 1 / ( DWPK - 56 ) + LN ( TMPK / DWPK ) / 800 )] + 56 LCLP Pressure (hPa) at the LCL, the lifting condensation level, from an average of the lowest 500 meters. LCLP = PRES * ( LCLT / ( TMPC + 273.15 ) ) ** ( 1 / KAPPA ) Poisson's equation
In [437]:
def PoissonConstant(wvmr):
"""http://glossary.ametsoc.org/wiki/Poisson_constant
May need to tweak low limit for dry air (=0.2854)"""
return where(wvmr>0., 0.2854*(1-0.24*wvmr), 0.2854)
In [20]:
def LFCParcelAscent(parcel_profile, temp_k, pres):
lfc_idx = nanmax(where((parcel_profile>temp_k) & (pres<(nanmax(pres)-50)))[0])
#lfc_temp = temp_k[lfc_idx]
lfc_pres = pres[lfc_idx]
# If parcel unstable throughout sounding set LFC to LCL
if all(parcel_profile>temp_k):
#lfc_temp = lcl_temp
lfc_pres = lcl_pres
return lfc_pres
In [19]:
def EQLVParcelAscent(parcel_profile, temp_k, pres):
idx_saddles=[]
for i in (where((parcel_profile>temp_k))[0]):
if i in (where(parcel_profile<temp_k)[0]+1):
#print i
idx_saddles.append(i)
#print idx_saddles
idx_eqlv = nanmin(idx_saddles)
#print idx_eqlv
# Interpolate around index of highest saddle point to zero difference
y_points_zero=0 # Points from original sounding interpolation
eqlv_p = interp_sounding(p[idx_eqlv-1:idx_eqlv+1],(parcel_profile-temp_k)[idx_eqlv-1:idx_eqlv+1],y_points_zero)
eqlv_t = interp_sounding(temp_k[idx_eqlv-1:idx_eqlv+1],(parcel_profile-temp_k)[idx_eqlv-1:idx_eqlv+1],y_points_zero)
return eqlv_p, eqlv_t
In [18]:
def ParcelAscentDryToLCLThenMoistC(init_parcel_pres,init_parcel_temp_c,init_parcel_dew_temp_c, pres):
dry_parcel_TC = LiftDry(init_parcel_pres,init_parcel_temp_c,init_parcel_dew_temp_c, pres)
dry_parcel_TK = dry_parcel_TC+273.15
if (init_parcel_temp_c>100.):
init_parcel_temp_c = init_parcel_temp_c - 273.15
if (init_parcel_dew_temp_c>100.):
init_parcel_dew_temp_c = init_parcel_dew_temp_c - 273.15
lcl_temp = LiftedCondensationLevelTemp(init_parcel_temp_c+273.15, init_parcel_dew_temp_c+273.15)
i_idx = (abs(dry_parcel_TK-lcl_temp)).argmin()
# Find initial temp at LCL for moist parcel ascent
y_points_zero=0
moist_t_init = dry_parcel_TC[i_idx]
temp_moist_adi_parcel_above_lcl_c = LiftWet(moist_t_init,pres[0:i_idx])
#temp_moist_adi_parcel_above_lcl_k = temp_moist_adi_parcel_above_lcl_c + 273.15
# Combine dry parcel ascent below LCL and moist above
parcel_profile=concatenate((temp_moist_adi_parcel_above_lcl_c, dry_parcel_TC[i_idx::]))
return parcel_profile
In [17]:
def ParcelAscentDryToLCLThenMoistK(init_parcel_pres,init_parcel_temp_c,init_parcel_dew_temp_c, pres):
dry_parcel_TC = LiftDry(init_parcel_pres,init_parcel_temp_c,init_parcel_dew_temp_c, pres)
dry_parcel_TK = dry_parcel_TC+273.15
if (init_parcel_temp_c>100.):
init_parcel_temp_c = init_parcel_temp_c - 273.15
if (init_parcel_dew_temp_c>100.):
init_parcel_dew_temp_c = init_parcel_dew_temp_c - 273.15
lcl_temp = LiftedCondensationLevelTemp(init_parcel_temp_c+273.15, init_parcel_dew_temp_c+273.15)
i_idx = (abs(dry_parcel_TK-lcl_temp)).argmin()
# Find initial temp at LCL for moist parcel ascent
y_points_zero=0
moist_t_init = dry_parcel_TC[i_idx]
temp_moist_adi_parcel_above_lcl_c = LiftWet(moist_t_init,pres[0:i_idx])
temp_moist_adi_parcel_above_lcl_k = temp_moist_adi_parcel_above_lcl_c + 273.15
# Combine dry parcel ascent below LCL and moist above
parcel_profile=concatenate((temp_moist_adi_parcel_above_lcl_k, dry_parcel_TK[i_idx::]))
return parcel_profile
In [16]:
# CAPE and CIN
def CapeCin(pres, temp_k, dew_temp_k, height, st_height):
# Calculate pressure, temperature and dew point temperature averages for first 500m
temp_k = TempInKelvinCheck(temp_k)
dew_temp_k = TempInKelvinCheck(dew_temp_k)
mean_temp_k = MeanFirst500m(temp_k, height, st_height)
dew_mean_temp_k = MeanFirst500m(dew_temp_k, height, st_height)
mean_pres = MeanFirst500m(pres, height, st_height)
#temp_c=KelvinToCelsius(temp_k)
#dew_temp_c=KelvinToCelsius(dew_temp_k)
mean_temp_c = KelvinToCelsius(mean_temp_k)
dew_mean_temp_c = KelvinToCelsius(dew_mean_temp_k)
# Find LCL temp and pressure
lcl_t = LiftedCondensationLevelTemp(mean_temp_k, dew_mean_temp_k)
vap_press = VaporPressure(dew_mean_temp_c)
wvmr = MixRatio(vap_press,mean_pres*100.)
kappa=PoissonConstant(wvmr)
lcl_p = LiftedCondensationLevelPres(mean_pres, lcl_t, mean_temp_k, kappa)
# Calculate dry parcel ascent
parcel_profile = ParcelAscentDryToLCLThenMoistK(mean_pres,mean_temp_c,dew_mean_temp_c, pres)
# Find equilibrium level
eqlv_p, eqlv_t = EQLVParcelAscent(parcel_profile, temp_k, pres)
# Find LFC
lfc_p = LFCParcelAscent(parcel_profile, temp_k, pres)
# Calculate CAPE and CIN
delta_z=diff(height) # Find delta height in sounding
# Take all but lowest pressure level (so length matches delta_z)
pp_diff = parcel_profile[1::]
Tk_diff = temp_k[1::]
p_diff=pres[1::]
sum_ascent = abs(delta_z)*(pp_diff-Tk_diff)/Tk_diff
CAPE = grav*sum(sum_ascent[((pp_diff-Tk_diff)>0) & (p_diff>eqlv_p) & (p_diff<lfc_p)])
#Taking levels above lcl but uwyo specifies top of mixed layer
CIN = grav*sum(sum_ascent[((pp_diff-Tk_diff)<0) & (p_diff>lfc_p) & (p_diff<lcl_p)])
#print "LCL_Test %s" % lcl_test
#print "Mean Pressure 500m %s" % mean_pres
#print "Dew Mean Temp 500m C %s" %dew_mean_temp_c
#print "Mean Temp 500m C %s" %mean_temp_c
#print "Vapour Pressure %s" %vap_press
#print "LCL Pressure %s" % lcl_p
#print "LFC Pressure %s" % lfc_p
#print "EQLV Pressure %s" % eqlv_p
#print "CAPE %s" % CAPE
#print "CIN %s" % CIN
return eqlv_p, parcel_profile,lcl_p, lfc_p, lcl_t, delta_z, CAPE, CIN
In [15]:
def CapeCinPBLInput(pres, temp_k, dew_temp_k, height, st_height, pbl_pressure): # PBL_pressure input temporary ?
# Calculate pressure, temperature and dew point temperature averages for first 500m
temp_k = TempInKelvinCheck(temp_k)
dew_temp_k = TempInKelvinCheck(dew_temp_k)
mean_temp_k = MeanFirst500m(temp_k, height, st_height)
dew_mean_temp_k = MeanFirst500m(dew_temp_k, height, st_height)
mean_pres = MeanFirst500m(pres, height, st_height)
#temp_c=TempInKelvinCheck(temp_k)
#dew_temp_c=TempInKelvinCheck(dew_temp_k)
mean_temp_c = mean_temp_k - 273.15
dew_mean_temp_c = dew_mean_temp_k - 273.15
# Find LCL temp and pressure
lcl_t = LiftedCondensationLevelTemp(mean_temp_k, dew_mean_temp_k)
vap_press = VaporPressure(dew_mean_temp_c)
wvmr = MixRatio(vap_press,mean_pres*100.)
kappa=PoissonConstant(wvmr)
lcl_p = LiftedCondensationLevelPres(mean_pres, lcl_t, mean_temp_k, kappa)
# Calculate dry parcel ascent
parcel_profile = ParcelAscentDryToLCLThenMoistC(mean_pres,mean_temp_c,dew_mean_temp_c, pres)
parcel_profile = parcel_profile+273.15
# Find equilibrium level
eqlv_p, eqlv_t = EQLVParcelAscent(parcel_profile, temp_k, pres)
# Find LFC
lfc_p = LFCParcelAscent(parcel_profile, temp_k, pres)
# Calculate CAPE and CIN
delta_z=diff(height) # Find delta height in sounding
# Take all but lowest pressure level (so length matches delta_z)
pp_diff = parcel_profile[1::]
Tk_diff = temp_k[1::]
p_diff=pres[1::]
sum_ascent = abs(delta_z)*(pp_diff-Tk_diff)/Tk_diff
CAPE = grav*sum(sum_ascent[((pp_diff-Tk_diff)>.) & (p_diff>eqlv_p) & (p_diff<lfc_p)])
#Taking levels above lcl but uwyo specifies top of mixed layer
#CIN = grav*sum(sum_ascent[((pp_diff-Tk_diff)<0) & (p_diff>lfc_p) & (p_diff<lcl_p)])
CIN = grav*sum(sum_ascent[((pp_diff-Tk_diff)<0.) & (p_diff>lfc_p) & (p_diff<pbl_pressure)])
#print "LCL_Test %s" % lcl_test
#print "Mean Pressure 500m %s" % mean_pres
#print "Dew Mean Temp 500m C %s" %dew_mean_temp_c
#print "Mean Temp 500m C %s" %mean_temp_c
#print "Vapour Pressure %s" %vap_press
#print "LCL Pressure %s" % lcl_p
#print "LFC Pressure %s" % lfc_p
#print "EQLV Pressure %s" % eqlv_p
#print "CAPE %s" % CAPE
#print "CIN %s" % CIN
return eqlv_p, parcel_profile,lcl_p, lfc_p, lcl_t, delta_z, CAPE, CIN
In [444]:
variable='pbl_pressure'
var_index = variable_name_index_match(variable, variable_list_line)
pbl_pressure = load_file['date_bin_mean_all_dates_one_station_single'][bin,0,var_index]
print pbl_pressure
In [ ]:
EQLV, pp, lclp,lfcp, lclt, delta_z, CAPE, CIN=CapeCin(p, Tk, Tdk, h, st_height)
You could compare the potential temperature of the low level parcel with that of the environment at every level, and where it becomes bigger than a threshold (I used 0.7 K), that's the CBL top
Need to compare LCL calc in interp_climat and CapeCin, 920 vs 890
Calculate PBL for CIN
Want to graph CAPE CIN etc for station by time
Want to graph time of sounding histogram and frequency
In [490]:
# Parse the data
import datetime
from dateutil.relativedelta import relativedelta
import re
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
def station_info_search(stat):
for line in station_metadata:
if "%s" % stat in line:
st = line[2].lower().title().replace('_',' ')
lo = float(line[3])
la = float(line[4])
st_height = float(line[5])
return st,la,lo, st_height
station_list_search='/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/igra-stations.txt'
station_metadata=[]
f = open(station_list_search,'r')
for line in f:
line = line.strip()
line=re.sub(r'([A-Z])\s([A-Z])', r'\1_\2',line)
line=re.sub(r'([A-Z])\s\s([A-Z])', r'\1_\2',line)
station_metadata.append(line.split())
f.close()
station_list_cs=[43003, 43014, 42867, 43371, 43353, 43285, 43192, 43150, 42339, 40990, 40948]
#station_list_cs=[43014]
grid=[]
date_min_month_bin=datetime.datetime(2011,5,1,0,0,0)
date_max_month_bin=datetime.datetime(2011,10,1,0,0,0)
delta = relativedelta(weeks=+1)
d = date_min_month_bin
while (d <= date_max_month_bin):
grid.append(d.timetuple().tm_yday)
d += delta
#years = mdates.YearLocator() # every year
months = mdates.MonthLocator() # every month
dateFmt = mdates.DateFormatter('%B %d')
variable='pbl_pressure'
var_index = variable_name_index_match(variable, variable_list_line)
for stat in station_list_cs:
load_file = load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section_'
'IND_SOUNDING_INTERP_MEAN_Climat_19600501_20111001_relativedelta(days=+7)_%s.npz' % stat)
cape_time_var=[]
cin_time_var=[]
date_plot_no_empty=[]
date_plot = [datetime.datetime.fromordinal(d).replace(year=2011) for d in grid]
for i in range(load_file['date_bin_mean_all_dates_one_station'].shape[0]):
pbl_pressure = load_file['date_bin_mean_all_dates_one_station_single'][i,0,var_index]
try:
data=load_file['date_bin_mean_all_dates_one_station']
p=data[i,0,:]/100
T=data[i,1,:]-273.15
Td=T-data[i,2,:]
h=data[i,15,:]*10
#print T
#print p
#print Td
station_name,la,lo, st_height=station_info_search(stat)
#u_wind,v_wind = u_v_winds(data[i,3,:], data[i,4,:])
print
EQLV, pp, lclp,lfcp, lclt, delta_z, CAPE, CIN=CapeCinPBLInput(p, T, Td, h, st_height, pbl_pressure/100)
cape_time_var.append(CAPE)
cin_time_var.append(CIN)
date_plot_no_empty.append(date_plot[i])
except Exception,e:
PrintException()
#print stat
#print p
#print Tk
#print Tdk
#print h
#print st_height
#print cape_time_var
#print cin_time_var
try:
print "Station Height % s" % st_height
station_name,la,lo, st_height=station_info_search(stat)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot_date(date_plot_no_empty, cape_time_var, fmt="-", label='CAPE')
ax.plot_date(date_plot_no_empty, cin_time_var, fmt="-", label='CIN')
plt.legend()
plt.title('%s' % station_name)
# format the ticks
#ax.xaxis.set_major_locator(months)
ax.xaxis.set_major_formatter(dateFmt)
#ax.xaxis.set_minor_locator(months)
fig.autofmt_xdate()
#fig.show()
fig.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_CAPE_CIN_Climatology_Station_Time_Plot.png'\
% (stat, date_min_month_bin.strftime('%Y%m%d'), date_max_month_bin.strftime('%Y%m%d')),\
format='png', bbox_inches='tight')
except Exception, e:
PrintException()
In [489]:
grid=[]
date_min=datetime.datetime(1960,5,1,0,0,0)
date_max=datetime.datetime(2014,10,1,0,0,0)
delta = relativedelta(weeks=+1)
date_min_month_bin=datetime.datetime(2011,5,1,0,0,0)
date_max_month_bin=datetime.datetime(2011,10,1,0,0,0)
delta = relativedelta(weeks=+1)
d = date_min_month_bin
while (d <= date_max_month_bin):
grid.append(d.timetuple().tm_yday)
d += delta
for stat in station_list_cs:
load_file = np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section_'
'IND_SOUNDING_INTERP_MEAN_Climat_%s_%s_%s_%s.npz'
% (date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), delta, stat))
load_file['date_bin_mean_all_dates_one_station_single'][:,0,1]
date_plot = [datetime.datetime.fromordinal(d).replace(year=2011) for d in grid]
station_name,la,lo, st_height=station_info_search(stat)
fig = plt.figure()
ax = fig.add_subplot(111)
variable='pbl_pressure'
var_index = variable_name_index_match(variable, variable_list_line)
ax.plot_date(date_plot[1:], load_file['date_bin_mean_all_dates_one_station_single'][1:,0,var_index], fmt="-", label='PBL')
variable='lcl_vpt'
var_index = variable_name_index_match(variable, variable_list_line)
ax.plot_date(date_plot[1:], load_file['date_bin_mean_all_dates_one_station_single'][1:,0,var_index], fmt="-", label='LCL')
variable='surface_pressure'
var_index = variable_name_index_match(variable, variable_list_line)
ax.plot_date(date_plot[1:], load_file['date_bin_mean_all_dates_one_station_single'][1:,0,var_index], fmt="-", label='SP')
variable='T_eq_0'
var_index = variable_name_index_match(variable, variable_list_line)
ax.plot_date(date_plot[1:], load_file['date_bin_mean_all_dates_one_station_single'][1:,0,var_index], fmt="-", label='T_eq_0')
plt.legend()
plt.title('%s' % station_name)
plt.ylim(50000,100500)
plt.gca().invert_yaxis()
# format the ticks
#ax.xaxis.set_major_locator(months)
ax.xaxis.set_major_formatter(dateFmt)
#ax.xaxis.set_minor_locator(months)
fig.autofmt_xdate()
#fig.show()
fig.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_PBL_etc_Climatology_Station_Time_Plot.png'\
% (stat, date_min_month_bin.strftime('%Y%m%d'), date_max_month_bin.strftime('%Y%m%d')),\
format='png', bbox_inches='tight')