diff options
| -rw-r--r-- | skewt.py | 53 | ||||
| -rwxr-xr-x | xar.py | 31 |
2 files changed, 74 insertions, 10 deletions
@@ -1,6 +1,8 @@ import matplotlib.gridspec as gridspec import matplotlib.pyplot as plt +import numpy as np + import metpy.calc as mpcalc from metpy.cbook import get_test_data from metpy.plots import add_metpy_logo, Hodograph, SkewT @@ -12,6 +14,8 @@ class Skewt: self._T = T self._Td = Td + self._info_lines = [] + # Create a new figure. The dimensions here give a good aspect ratio self._fig = plt.figure(figsize=(9, 9)) @@ -46,20 +50,55 @@ class Skewt: plt.ylabel('$m/s$') self._skew.plot_barbs(self._p, u, v) - def addInfo(self, lines): - # TODO + def addInfo(self, line): + self._info_lines.append(line) + + def addAnalysis(self, analysis='ccl', shade=False): + f = {'ccl': self._cclAnalysis, 'lcl': self._lclAnalysis} + + lvl, parcel = f[analysis]() + + self._skew.plot(self._p, parcel, 'y') + self._skew.plot(lvl[0], lvl[1], 'o', markerfacecolor='red', linewidth=1) + + cape, cin = mpcalc.cape_cin(self._p, self._T, self._Td, parcel) + + self.addInfo(f'CAPE {int(cape.magnitude)} $J/kg$ CIN {int(cin.magnitude)} $J/kg$') + + if shade: + self._skew.shade_cape(self._p,self._T,parcel) + self._skew.shade_cin(self._p,self._T,parcel) + + def _cclAnalysis(self): + #p = np.arange(max(self._p).magnitude, min(self._p).magnitude, -50) * units.hPa + + ccl = mpcalc.ccl(self._p,self._T,self._Td) + ccl_ground=mpcalc.dry_lapse(self._p[:1], ccl[1], ccl[0]) + ccl_ground_parcel= mpcalc.parcel_profile(self._p, ccl_ground[0], self._Td[0]) + + return (ccl, ccl_ground_parcel) + + def _lclAnalysis(self): + ground_parcel= mpcalc.parcel_profile(self._p, self._T[0], self._Td[0]) + lcl = mpcalc.lcl(self._p[0],self._T[0],self._Td[0]) + + return (lcl, ground_parcel) + + def _buildInfoBox(self): ax = self._fig.add_subplot(self._gs[1,-1]) - ax.text(0, 0, lines, ha='left', va='center', size=12) + ax.text(0, 0, '\n'.join(self._info_lines), ha='left', va='center', + size=10, fontfamily='monospace') ax.axis("off") - def plot(self, filename=None): + self._buildInfoBox() + # Add the relevant special lines #self._skew.ax.set_ylim(max(self._p), min(self._p)) self._skew.ax.set_ylim(1000, 100) - self._skew.plot_dry_adiabats() - self._skew.plot_moist_adiabats() - self._skew.plot_mixing_lines() + self._skew.plot_dry_adiabats(linewidth=1) + self._skew.plot_moist_adiabats(linewidth=1) + self._skew.plot_mixing_lines(linewidth=1) # Good bounds for aspect ratio self._skew.ax.set_xlim(-30, 40) @@ -4,8 +4,19 @@ import xarray as xr from metpy.units import units import metpy.calc as mpcalc +import numpy as np + +import datetime + import skewt +def np_time_convert(dt64): + unix_epoch = np.datetime64(0, 's') + one_second = np.timedelta64(1, 's') + seconds_since_epoch = (dt64 - unix_epoch) / one_second + + return datetime.datetime.utcfromtimestamp(seconds_since_epoch) + #grib = pygrib.open('dwd_icon-d2/combined.grib2') data = xr.load_dataset('dwd_icon-d2/combined.grib2', engine='cfgrib') @@ -26,9 +37,23 @@ for step in for_temp.coords['step']: u = this_step.u.values * (units.m / units.s) v = this_step.v.values * (units.m / units.s) + valid = np_time_convert(step.valid_time.values) + init = np_time_convert(step.time.values) + + valid_str = valid.strftime('%d %b %Y - %HUTC') + init_str = init.strftime('%d %b %Y - %HUTC') + hours_since_init_str = str(int(this_step.step.values / np.timedelta64(1,'h'))).zfill(2) + skt = skewt.Skewt(p=p, T=T, Td=Td) skt.addWindUV(u, v) - skt.addInfo("TEST") - skt.plot(filename=f'skewt.png') + skt.addInfo(f"VALID: {valid_str}") + skt.addInfo(f"INIT : {init_str}") + skt.addInfo(f"LAT {lat} LON {lon}") + skt.addInfo('INIT+' + hours_since_init_str) + skt.addInfo(f'') + skt.addAnalysis(shade=True, analysis='lcl') + skt.addInfo("DWD ICON-D2") + + init_for_filename = init.strftime('%Y-%m-%d-%H') - break + skt.plot(filename=f'skewt_{init_for_filename}+{hours_since_init_str}.png') |