search

Unidata / MetPy

MetPy is a collection of tools in Python for reading, visualizing and performing calculations with weather data.
https://unidata.github.io/MetPy/
BSD 3-Clause "New" or "Revised" License
1.26k stars 415 forks source link

Calculated negative cape values from Buffalo sounding #1336

Open eliteuser26 opened 4 years ago

eliteuser26 commented 4 years ago

The calculation of the CAPE and CIN areas was done from the Buffalo sounding on March 13, 2020 at 12 UTC with this method in Metpy:

cape_value, cin_value = mpcalc.cape_cin(p, T, Td, parcel_path)

Here are the values which resulted from this method:

CAPE: -192.87189038978863 J/kg CIN: -4.246779773060985 J/kg

This is the sample code used to calculate these values above:

from datetime import datetime
from siphon.simplewebservice.wyoming import WyomingUpperAir
from metpy.units import units
import matplotlib.pyplot as plt
import metpy.plots as plots
import numpy as np
import metpy.calc as mpcalc
import pandas as pd
%matplotlib inline

# set date and upper air station name to retrieve in Siphon
date = datetime(2020, 3, 13, 12)
date_string = date.strftime('%d %B %Y %X UTC')
station = 'BUF'
df = WyomingUpperAir.request_data(date, station)

# specify different parameters to display on Skew-T
p = df['pressure'].values * units(df.units['pressure'])
T = df['temperature'].values * units(df.units['temperature'])
Td = df['dewpoint'].values * units(df.units['dewpoint'])
u = df['u_wind'].values * units(df.units['u_wind'])
v = df['v_wind'].values * units(df.units['v_wind'])
# heights = df['height'].values * units(df.units['height'])

# create figure and indicate size
fig = plt.figure(figsize=(10,10))
# create SkewT plot and link to figure
skew = plots.SkewT(fig)
# plot temperature in red and dewpoint in green vs pressure on SkewT
skew.plot(p, T, 'red')
skew.plot(p, Td, 'green')
# create array of log space (default 50 points) using pressure level using base 10
interval = np.logspace(2,3) * units.hPa
# find nearest neighbour for pressure level based on log space
idx = mpcalc.resample_nn_1d(p, interval[1:])
# plot wind barbs vs pressure based on log space on SkewT
skew.plot_barbs(p[idx], u[idx], v[idx])
# for winter set xlimit between -60 and 30C
skew.ax.set_xlim(-60,30)
# plot dry and moist adiabat lines and mixing lines on SkewT
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()
# add MetPy logo and indicate position on figure
plots.add_metpy_logo(fig, x=115, y=95)
# add station and valid datetime to plot titles with locations
plt.title(f'{station} Sounding', loc='left')
plt.title(f'Valid Time: {date_string}', loc='right')
# draw black vertical skew line of 0C on SkewT
skew.ax.axvline(0, color='k')
# calculate parcel path
parcel_path = mpcalc.parcel_profile(p, T[0], Td[0])
# plot parcel path on SkewT
skew.plot(p, parcel_path, color='k')
# shade cape and cin on SkewT between parcel path and temperature
skew.shade_cape(p, T, parcel_path)
skew.shade_cin(p, T, parcel_path)
dopplershift commented 4 years ago

I can confirm on my system with this:

from datetime import datetime
import metpy.calc as mpcalc
from metpy.units import pandas_dataframe_to_unit_arrays
from siphon.simplewebservice.wyoming import WyomingUpperAir

df = WyomingUpperAir.request_data(datetime(2020, 3, 13, 12), 'BUF')
data = pandas_dataframe_to_unit_arrays(df)
mpcalc.surface_based_cape_cin(data['pressure'], data['temperature'], data['dewpoint'])

gives on current metpy master:

(-192.87189038978863 <Unit('joule / kilogram')>,
 -5.391129989300722 <Unit('joule / kilogram')>)
eliteuser26 commented 4 years ago

It is one of those rare cases where the value for Cape is negative. Normally it is giving the proper positive value as intended. This seems to occur when the parcel profile is very close to the measured temperatures in the lower atmosphere.

eliteuser26 commented 4 years ago

Found another case of negative Cape value for Detroit DTX on April 1st at 12 UTC. This isn't an April fool's joke. The calculated Cape value is -27 J/kg. It seems to occur when the parcel profile is crisscrossing the measured temperature profile.

akrherz commented 4 years ago

I've been throwing my entire RAOB archive through MetPy's sounding functions and did a quick check for negative sbcape for the past few days and behold

station valid sbcape_jkg
KDTX 2020-04-01 12:00:00+00 -27.533386
KXMR 2020-04-01 00:00:00+00 -30.603716
KKPP 2020-03-31 12:00:00+00 -95.12861
KMFL 2020-03-31 12:00:00+00 -390.22668
KCUN 2020-03-31 00:00:00+00 -179.90114
KDDC 2020-03-31 00:00:00+00 -81.269775
KJAX 2020-03-31 00:00:00+00 -21.090246
KJSJ 2020-03-31 00:00:00+00 -233.70105
KCUN 2020-03-30 12:00:00+00 -263.11
KJSJ 2020-03-30 12:00:00+00 -397.9066
KMFL 2020-03-30 12:00:00+00 -190.04387
KSDQ 2020-03-30 12:00:00+00 -131.27994
KSLE 2020-03-30 12:00:00+00 -121.92557
KXKF 2020-03-30 12:00:00+00 -352.6673
KBDI 2020-03-29 12:00:00+00 -266.86176
KBMX 2020-03-29 12:00:00+00 -22.72841
KPIT 2020-03-29 12:00:00+00 -12.695718
KTLH 2020-03-29 12:00:00+00 -401.89224
KBMX 2020-03-28 12:00:00+00 -671.2816
KFFC 2020-03-28 12:00:00+00 -775.76196
KJAN 2020-03-28 12:00:00+00 -427.56732

Overall, 435 out of 22,553 soundings with negative sbcape this year :)

dopplershift commented 4 years ago

Well, I'm not thrilled, but at least that's not an overwhelming amount. Thanks for doing that test @akrherz !

eliteuser26 commented 4 years ago

I haven't checked other sites as I am only using 3 locations for my website which are closest to my position. It tries to calculate the area between the parcel profile and the measured temperatures. However it seems to combine both CAPE and CIN values together to calculate the resulting value. When the profile and measured temperatures are close together it can't determine the LFC along the profile.

eliteuser26 commented 4 years ago

Still seeing more negative Cape values from Buffalo when lifted parcel line is very close to the environmental values. This can be fixed easily in the Metpy code by assigning zero for negative values. Same type of code as for CIN but reversed.

This can also occur when there are 2 areas of CIN and 2 areas of Cape in the same profile. In this case it will be calculated as a negative values if CIN values are larger than Cape values.

dopplershift commented 4 years ago

Last time we dug into that code, I didn't think it should be able to return negative values. Setting anything <0 to 0 sounds like a workaround for conditions that break internal assumptions--I'd much rather fix the code itself, assuming my understanding is correct.

eliteuser26 commented 4 years ago

After thinking about it and looking at the parcel path I would agree with you. The same would apply for positive CIN. I was thinking about looking at the code and results to find out why this is happening. In those cases I found that they were multiple areas of Cape and CIN very close to the environmental profile. Otherwise all other cases are being calculated correctly.