Closed ahill818 closed 5 years ago
stefan-hofer
commented
5 years ago Just by eyeballing: I'm surprised that there is an LFC in the second case at 815 hPa. Doesn't look to me that the parcel profile (black) is warmer than (or right of) the ambient air temperature. What is the ambient profile doing near the surface in the second case? It seems there might be some temperature values missing between 1000 hPa and 920 hPa.
ahill818
commented
5 years ago I'm going to guess it's just dry adiabatic 1000-920 hPa in the second case, so you can't actually see the parcel T profile under the black line.
And you are probably right about the LFC; I have a couple try/except statements in the code above that set the LFC to the LCL in the case where mpcalc.lfc() fails, just so I can generate the plot
ahill818
commented
5 years ago OK I've tried to create a minimum working example for our purposes. I've attached a notebook that handles some data from a model sounding and plots the profile and calculates statistics.
# coding: utf-8
# In[1]:
get_ipython().magic('matplotlib inline')
#
# Hodograph Inset
# ===============
#
# Layout a Skew-T plot with a hodograph inset into the plot.
#
#
# In[2]:
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import numpy as np
import pandas as pd
import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.plots import add_metpy_logo, Hodograph, SkewT
from metpy.units import units
# Upper air data can be obtained using the siphon package, but for this example we will use
# some of MetPy's sample data.
#
#
# In[3]:
col_names = ['pressure', 'height', 'temperature', 'dewpoint', 'direction', 'speed']
df = pd.read_fwf(get_test_data('may4_sounding.txt', as_file_obj=False),
skiprows=5, usecols=[0, 1, 2, 3, 6, 7], names=col_names)
df['u_wind'], df['v_wind'] = mpcalc.wind_components(df['speed'],
np.deg2rad(df['direction']))
# Drop any rows with all NaN values for T, Td, winds
df = df.dropna(subset=('temperature', 'dewpoint', 'direction', 'speed',
'u_wind', 'v_wind'), how='all').reset_index(drop=True)
# We will pull the data out of the example dataset into individual variables and
# assign units.
#
#
# In[4]:
p = df['pressure'].values * units.hPa
T = df['temperature'].values * units.degC
Td = df['dewpoint'].values * units.degC
wind_speed = df['speed'].values * units.knots
wind_dir = df['direction'].values * units.degrees
u, v = mpcalc.wind_components(wind_speed, wind_dir)
z = (df['height'].values * units.meters).to(units.km)
# In[44]:
p = [ 902.1554,
897.9034,
893.6506,
889.4047,
883.063,
874.6284,
866.2387,
857.887,
849.5506,
841.2686,
833.0042,
824.7891,
812.5049,
796.2104,
776.0027,
751.9025,
727.9612,
704.1409,
680.4028,
656.7156,
629.077,
597.4286,
565.6315,
533.5961,
501.2452,
468.493,
435.2486,
401.4239,
366.9387,
331.7026,
295.6319,
258.6428,
220.9178,
182.9384,
144.959,
106.9778,
69.00213] * units.hPa
t = [ -3.039381,
-3.703779,
-4.15996,
-4.562574,
-5.131827,
-5.856229,
-6.568434,
-7.276881,
-7.985013,
-8.670911,
-8.958063,
-7.631381,
-6.05927,
-5.083627,
-5.11576,
-5.687552,
-5.453021,
-4.981445,
-5.236665,
-6.324916,
-8.434324,
-11.58795,
-14.99297,
-18.45947,
-21.92021,
-25.40522,
-28.914,
-32.78637,
-37.7179,
-43.56836,
-49.61077,
-54.24449,
-56.16666,
-57.03775,
-58.28041,
-60.86264,
-64.21677] * units.degC
td = [ -22.08774,
-22.18181,
-22.2508,
-22.31323,
-22.4024,
-22.51582,
-22.62526,
-22.72919,
-22.82095,
-22.86173,
-22.49489,
-21.66936,
-21.67332,
-21.94054,
-23.63561,
-27.17466,
-31.87395,
-38.31725,
-44.54717,
-46.99218,
-43.17544,
-37.40019,
-34.3351,
-36.42896,
-42.1396,
-46.95909,
-49.36232,
-48.94634,
-47.90178,
-49.97902,
-55.02753,
-63.06276,
-72.53742,
-88.81377,
-93.54573,
-92.92464,
-91.57479] * units.degC
wind_speed = df['speed'].values * units.knots
wind_dir = df['direction'].values * units.degrees
u = [ -4.548029,
-5.060958,
-5.170549,
-5.172879,
-5.104956,
-4.918171,
-4.668983,
-4.353795,
-3.940358,
-3.283814,
-1.323716,
2.205689,
6.328888,
9.783625,
11.63874,
12.85858,
15.08335,
18.16523,
21.24494,
24.22534,
28.55559,
33.66334,
38.19795,
41.22777,
42.92673,
44.49194,
47.182,
52.06959,
59.53011,
65.812,
71.88536,
83.87706,
88.8103,
88.62223,
86.58784,
73.31921,
49.98649] * units.knots
v = [ 15.69195,
18.00935,
18.82138,
19.23343,
19.59081,
19.77925,
19.87817,
19.98829,
20.2242,
20.90455,
23.37206,
24.52552,
23.4171,
20.71219,
16.74612,
12.51609,
12.70995,
14.9863,
15.97422,
16.03022,
16.64295,
18.01906,
19.80627,
21.36071,
22.69773,
24.44924,
26.87852,
30.4185,
35.47622,
42.26982,
52.09302,
64.79327,
66.51744,
53.17516,
38.22067,
36.1087,
20.14054] * units.knots
z = ([ 18.70802,
56.1215,
93.54421,
130.9969,
187.2617,
262.4107,
337.7491,
413.2926,
489.058,
565.064,
641.3884,
718.3456,
835.6757,
994.3338,
1196.642,
1443.93,
1697.404,
1958.339,
2227.387,
2504.746,
2840.104,
3237.853,
3654.145,
4092.083,
4555.602,
5049.623,
5580.284,
6154.932,
6781.671,
7470.441,
8236.378,
9106.034,
10118.53,
11325.17,
12813.3,
14757.47,
17601.05] * units.meters).to(units.km)
# In[28]:
z[z.magnitude < 8.0]
# In[45]:
# Create a new figure. The dimensions here give a good aspect ratio
fig = plt.figure(figsize=(9, 9))
# Grid for plots
skew = SkewT(fig, rotation=45)
ax1 = fig.add_axes([0.98,0.2,0.1,0.3])
ax2 = fig.add_axes([1.14,0.2,0.1,0.3])
ax1.axis('off')
ax2.axis('off')
ax_hod = fig.add_axes([0.96, 0.58, 0.3, 0.3])
skew.ax.set_ylim(1000, 100)
skew.ax.set_xlim(-40, 60)
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()
# Plot the data using normal plotting functions, in this case using
# log scaling in Y, as dictated by the typical meteorological plot
skew.plot(p, t, 'r')
skew.plot(p, td, 'g')
skew.plot_barbs(p, u, v, xloc=0.97)
# Create a hodograph
h = Hodograph(ax_hod, component_range=80.)
h.add_grid(increment=20)
im = h.plot_colormapped(u.magnitude[z.magnitude < 8.0], v.magnitude[z.magnitude < 8.0], z[z.magnitude < 8.0])
cbpos = fig.add_axes([1.18, 0.48, 0.1, 0.5])
cbpos.axis('off')
fig.colorbar(im)
prof = mpcalc.parcel_profile(p, t[0], td[0]).to('degC')
lcl_p, lcl_t = mpcalc.lcl(p[0], t[0], td[0])
parcel_t_start = prof.T[0]
parcel_p_start = p[0]
lfc_p, lfc_t = mpcalc.lfc(p,t,td,prof)
el_p,el_t = mpcalc.el(p,t,td)
sbcape,sbcin = mpcalc.surface_based_cape_cin(p,t,td)
mucape,mucin = mpcalc.most_unstable_cape_cin(p,t,td)
ax1.axis('off')
ax1.text(0.01,1.0,'Parcel T Start:',size=8)
ax1.text(0.01,0.94,'Parcel P Start:',size=8)
ax1.text(0.01,0.88,'LCL Pressure:',size=8)
ax1.text(0.01,0.82,'LCL Temperature:',size=8)
ax1.text(0.01,0.76,'LFC Pressure:',size=8)
ax1.text(0.01,0.7,'LFC Temperature:',size=8)
ax1.text(0.01,0.64,'EL:',size=8)
ax1.text(0.01,0.58,'SBCAPE:',size=8)
ax1.text(0.01,0.52,'SBCIN:',size=8)
ax1.text(0.01,0.46,'MUCAPE:',size=8)
ax1.text(0.01,0.40,'MUCIN:',size=8)
ax2.text(0.01,1.0,'{0} degC'.format(parcel_t_start,size=8))
ax2.text(0.01,0.94,'{0} hPa'.format(parcel_p_start,size=8))
ax2.text(0.01,0.88,'{0} hPa'.format(lcl_p,size=8))
ax2.text(0.01,0.82,'{0} degC'.format(lcl_t,size=8))
ax2.text(0.01,0.76,'{0} hPa'.format(lfc_p,size=8))
ax2.text(0.01,0.7,'{0} degC'.format(lfc_t,size=8))
ax2.text(0.01,0.64,'{0} hPa'.format(el_p,size=8))
ax2.text(0.01,0.58,'{0} J/kg'.format(sbcape,size=8))
ax2.text(0.01,0.52,'{0} J/kg'.format(sbcin,size=8))
ax2.text(0.01,0.46,'{0} J/kg'.format(mucape,size=8))
ax2.text(0.01,0.40,'{0} J/kg'.format(mucin,size=8))
# Show the plot
plt.show()
dopplershift
commented
5 years ago I haven't had a chance to look at this in detail, but hopefully soon.
Otherwise: 😱
sgdecker
commented
5 years ago The problem is here: https://github.com/Unidata/MetPy/blob/bf1b701dee22c528a22913e96ad5d8f14a2f2784/metpy/calc/thermo.py#L412
The _less_or_close call should only look at the part of the profile above the LCL. In this sounding, the superadiabatic lapse rate at the surface leads to this function seeing the positive area that exists well below the LCL. Hence the "LFC = LCL" branch is erroneously taken. A PR is coming, although it might generate a merge conflict with my other PR.
sgdecker
commented
5 years ago The above fixes the LFC issue, but the EL still is calculated to be 836 mb for this data. I suppose the EL should be NaN if it is closer to the ground than the LCL.
dopplershift
commented
5 years ago With #1055 merged, I think we've fixed all of the issues identified here. I can run the sample code without error (so long as I have ping 0.8 🙄 ) Sorry it's taken so long. 😞
Please do re-open if you find that we haven't actually fixed them all.
ahill818
commented
5 years ago Thank you all for the hard work!
Kyl67899
commented
3 years ago How did you calculate the CAPEV, CINV and profile V? I am having trouble with those and there is no mp.calc for those. I am also using siphon but everything else is similar but some things are different.
HathewayWill
commented
1 year ago I have produced some weird quirks in calculating LFC temperature and pressure, along with EL, and resulting CAPE calculations. I am using WRF output to generate soundings, and calculate corresponding profile statistics (see below examples to see the quirks). In some instances, the LFC temperature, pressure and EL come out at NaNs; correspondingly, this produces LARGE negative CAPE values. Below are snippets of our code that grabs corresponding profiles from WRFOUT files, and inserts that data into MetPy functions to calculate profiles. I don't have a simple way to share this code for testing because these snippets are embedded in python postprocessing for our realtime deterministic modeling system.
import pandas as pd from metpy.units import units import metpy.calc as mpcalc def calc_sounding_stats(t,td,p,z,u,v): dictionary = {'pressure':p,'height':z,'temperature':t,'dewpoint':td,'u_wind':u,'v_wind':v} df = pd.DataFrame.from_dict(dictionary,orient='columns') s_p = df['pressure'].values * units.hPa s_t = df['temperature'].values * units.degC s_td = df['dewpoint'].values * units.degC s_u, s_v = df['u_wind'].values * units.knot, df['v_wind'].values * units.knot s_spd = np.sqrt(s_u**2 + s_v**2) s_z = df['height'].values * units.meter prof = mpcalc.parcel_profile(s_p, s_t[0], s_td[0]).to('degC') lcl_pressure, lcl_temperature = mpcalc.lcl(s_p[0], s_t[0], s_td[0]) parcel_t_start = prof.T[0] parcel_p_start = s_p[0] try: lfc_p,lfc_t = mpcalc.lfc(s_p,s_t,s_td,prof) except IndexError: lfc_p,lfc_t = lcl_pressure,lcl_temperature u_shear_1km,v_shear_1km = mpcalc.bulk_shear(s_p,s_u,s_v,depth=1000 * units.meter) bulk_shear_1km = np.round(np.sqrt(u_shear_1km**2 + v_shear_1km**2)) bulk_shear_1km_dir = np.round(mpcalc.wind_direction(u_shear_1km,v_shear_1km)) u_shear_6km,v_shear_6km = mpcalc.bulk_shear(s_p,s_u,s_v,depth=6000 * units.meter) bulk_shear_6km = np.round(np.sqrt(u_shear_6km**2 + v_shear_6km**2)) bulk_shear_6km_dir = np.round(mpcalc.wind_direction(u_shear_6km,v_shear_6km)) el_p,el_t = mpcalc.el(s_p,s_t,s_td) try: sbcape,sbcin = mpcalc.surface_based_cape_cin(s_p,s_t,s_td) except IndexError: sbcape,sbcin = 0.0*units.joule/units.kilogram,-999.9*units.joule/units.kilogram try: mucape,mucin = mpcalc.most_unstable_cape_cin(s_p,s_t,s_td) except IndexError: mucape,mucin = 0.0*units.joule/units.kilogram,-999.9*units.joule/units.kilogram rmover,lmover,bunk_mean = mpcalc.bunkers_storm_motion(s_p,s_u,s_v,s_z) srh_0to3_pos, srh_0to3_neg, srh_0to3_tot = mpcalc.storm_relative_helicity(s_u,s_v,s_z,3000 * units.meter, bottom=s_z[0]) srh_0to1_pos, srh_0to1_neg, srh_0to1_tot = mpcalc.storm_relative_helicity(s_u,s_v,s_z,1000 * units.meter, bottom=s_z[0]) pwat = mpcalc.precipitable_water(s_td,s_p) bunk_right_dir = np.round(mpcalc.wind_direction(rmover[0],rmover[1])) bunk_left_dir = np.round(mpcalc.wind_direction(lmover[0],lmover[1])) bunk_right_spd = np.round(np.sqrt(rmover[0]**2 + rmover[1]**2)) bunk_left_spd = np.round(np.sqrt(lmover[0]**2 + lmover[1]**2)) # return a list of calculated statistics return [parcel_t_start.magnitude,parcel_p_start.magnitude,lcl_pressure.magnitude,lcl_temperature.magnitude,lfc_p.magnitude,lfc_t.magnitude,el_p.magnitude,sbcape.magnitude,sbcin.magnitude,mucape.magnitude,mucin.magnitude,bulk_shear_1km.magnitude,bulk_shear_1km_dir.magnitude,bulk_shear_6km.magnitude,bulk_shear_6km_dir.magnitude,srh_0to1_tot.magnitude,srh_0to3_tot.magnitude,bunk_right_spd.magnitude,bunk_right_dir.magnitude,bunk_left_spd.magnitude,bunk_left_dir.magnitude,pwat.magnitude] # temp, td, p, z_agl, uc, and vc are numpy arrays from reading WRF files # function to grab model profile at an interpolated i,j location, returns 1d arrays representing vertical profile of t, td, p, z, u, and v snd_t,snd_td,snd_p,snd_z,snd_u,snd_v = grab_interp_sounding(i_loc,j_loc,temp-273.15,td-273.15,p,z_agl,uc*1.94844,vc*1.94844) # Calculate sounding statistics from the vertical profiles of variables next_snding = calc_sounding_stats(snd_t,snd_td,snd_p,snd_z,snd_u,snd_v) # This is the code that plots the sounding in the examples below def sounding_plot(t, td, p, z, u, v, stats, init_label, valid_label, title, time_int, outdir, plot_name): fig, skew, ax_hod, ax1, ax2 = create_fig_sounding() if (len(str(time_int)) == 1): outtime = 'f0' + str(time_int) else: outtime = 'f' + str(time_int) p = p * units.millibar t = t * units.degC td = td * units.degC z = z * units.meter u = u * units.knots v = v * units.knots p1 = skew.plot(p, t, 'r') p2 = skew.plot(p, td, 'g') p3 = skew.plot_barbs(p[::2], u[::2], v[::2]) prof = mpcalc.parcel_profile(p, t[0], td[0]).to('degC') p4 = skew.plot(p, prof, 'k', linewidth=2) p5 = skew.shade_cin(p, t, prof) p6 = skew.shade_cape(p, t, prof) P.sca(ax_hod) h = Hodograph(ax_hod, component_range=80.) h.add_grid(increment=20) wind_speed = np.sqrt(u**2 + v**2) h.plot_colormapped(u, v, wind_speed) P.sca(ax1) ax1.axis('off') ax1.text(0.01,1.0,'Parcel T Start:',size=8) ax1.text(0.01,0.94,'Parcel P Start:',size=8) ax1.text(0.01,0.88,'LCL Pressure:',size=8) ax1.text(0.01,0.82,'LCL Temperature:',size=8) ax1.text(0.01,0.76,'LFC Pressure:',size=8) ax1.text(0.01,0.7,'LFC Temperature:',size=8) ax1.text(0.01,0.64,'EL:',size=8) ax1.text(0.01,0.58,'SBCAPE:',size=8) ax1.text(0.01,0.52,'SBCIN:',size=8) ax1.text(0.01,0.46,'MUCAPE:',size=8) ax1.text(0.01,0.40,'MUCIN:',size=8) ax1.text(0.01,0.34,'0-1km Bulk Shear:',size=8) ax1.text(0.01,0.28,'0-6km Bulk Shear:',size=8) ax1.text(0.01,0.22,'0-1 SRH:',size=8) ax1.text(0.01,0.16,'0-3 SRH:',size=8) ax1.text(0.01,0.10,'Bunkers Right:',size=8) ax1.text(0.01,0.04,'Bunkers Left:',size=8) ax1.text(0.01,-0.02,'PWAT:',size=8) #ax1.text(0.01,-0.08,'SCP',size=8) #ax1.text(0.01,-0.14,'STP',size=8) #ax1.text(0.01,-0.14,'Lapse Rate',size=8) P.sca(ax2) ax2.axis('off') ax2.text(0.01,1.0,'{0} degC'.format(np.round(float(stats[0]),1)),size=8) ax2.text(0.01,0.94,'{0} hPa'.format(np.round(float(stats[1]),1)),size=8) ax2.text(0.01,0.88,'{0} hPa'.format(np.round(float(stats[2]),1)),size=8) ax2.text(0.01,0.82,'{0} degC'.format(np.round(float(stats[3]),1)),size=8) ax2.text(0.01,0.76,'{0} hPa'.format(np.round(float(stats[4]),1)),size=8) ax2.text(0.01,0.7,'{0} degC'.format(np.round(float(stats[5]),1)),size=8) ax2.text(0.01,0.64,'{0} hPa'.format(np.round(float(stats[6]),1)),size=8) ax2.text(0.01,0.58,'{0} J/kg'.format(np.round(float(stats[7]),1)),size=8) ax2.text(0.01,0.52,'{0} J/kg'.format(np.round(float(stats[8]),1)),size=8) ax2.text(0.01,0.46,'{0} J/kg'.format(np.round(float(stats[9]),1)),size=8) ax2.text(0.01,0.40,'{0} J/kg'.format(np.round(float(stats[10]),1)),size=8) ax2.text(0.01,0.34,'{0} kts at {1} deg'.format(float(stats[11]),float(stats[12])),size=8) ax2.text(0.01,0.28,'{0} kts at {1} deg'.format(float(stats[13]),float(stats[14])),size=8) ax2.text(0.01,0.22,'%s m$^{2}$ s$^{-2}$'%np.round(float(stats[15]),1),size=8) ax2.text(0.01,0.16,'%s m$^{2}$ s$^{-2}$'%np.round(float(stats[16]),1),size=8) ax2.text(0.01,0.10,'{0} kts at {1} deg'.format(float(stats[17]),float(stats[18])),size=8) ax2.text(0.01,0.04,'{0} kts at {1} deg'.format(float(stats[19]),float(stats[20])),size=8) ax2.text(0.01,-0.02,'{0} in'.format(np.round(float(stats[21])/25.4,2)),size=8) t1 = fig.text(0.13, 0.895, title, fontsize=10, fontweight='bold') t2 = fig.text(0.70, 0.915, init_label, fontsize=9) t3 = fig.text(0.66, 0.895, outtime + ', ' + valid_label, fontsize=9) fig_name = outdir + plot_name + '_' + outtime + '.png' P.savefig(fig_name, format="png", bbox_inches='tight') fig.clf()I could insert some if/else statements to identify the large negative CAPE and convert them to zero, but I would like to know if something else is going on with the profile calculations.
Using Python 3.6.3 :: Anaconda custom (64-bit) Metpy 0.9.1
@ahill818 So i'm trying to follow your example and i'm having issues. Do you have the working python code somewhere on your github? I'm trying to use wrf output data
I have produced some weird quirks in calculating LFC temperature and pressure, along with EL, and resulting CAPE calculations. I am using WRF output to generate soundings, and calculate corresponding profile statistics (see below examples to see the quirks). In some instances, the LFC temperature, pressure and EL come out at NaNs; correspondingly, this produces LARGE negative CAPE values. Below are snippets of our code that grabs corresponding profiles from WRFOUT files, and inserts that data into MetPy functions to calculate profiles. I don't have a simple way to share this code for testing because these snippets are embedded in python postprocessing for our realtime deterministic modeling system.
I could insert some if/else statements to identify the large negative CAPE and convert them to zero, but I would like to know if something else is going on with the profile calculations.
Using Python 3.6.3 :: Anaconda custom (64-bit) Metpy 0.9.1