From 7a4b8a91ce5b9bbd2ce806d2463e95634ad8b1a1 Mon Sep 17 00:00:00 2001
From: Jonas Gunz <himself@jonasgunz.de>
Date: Fri, 25 Aug 2023 12:23:02 +0200
Subject: skewt with CAPE/CIN and info

---
 skewt.py | 53 ++++++++++++++++++++++++++++++++++++++++++++++-------
 xar.py   | 31 ++++++++++++++++++++++++++++---
 2 files changed, 74 insertions(+), 10 deletions(-)

diff --git a/skewt.py b/skewt.py
index 43bcc91..0b8f4a7 100644
--- a/skewt.py
+++ b/skewt.py
@@ -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)
diff --git a/xar.py b/xar.py
index fca3887..13d5f62 100755
--- a/xar.py
+++ b/xar.py
@@ -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')
-- 
cgit v1.2.3