stuff

6 files changed, 141 insertions, 10 deletions

diff --git a/Readme.md b/Readme.md
index 32de887..38ed16c 100644
--- a/Readme.md
+++ b/Readme.md
@@ -3,3 +3,8 @@
 DTKE_CON for convection forecast?
 
 w_ctmax updraft velocity
+
+## Sources
+
+DWD OpenData: opendata.dwd.de
+ECMWF OpenData: https://www.ecmwf.int/en/forecasts/datasets/open-data
diff --git a/grib.py b/grib.py
index c298450..767f3d7 100755
--- a/grib.py
+++ b/grib.py
@@ -16,15 +16,34 @@ target_lon='10.0646952'
 
 grib = pygrib.open('dwd_icon-d2/combined.grib2')
 
+# GRIB-Objekt wichtige attribute:
+# analDate
+# validDate
+# units
+# level
+# typeOfLevel
+# name
+# shortName
 
 for grb in grib:
-    print(grb)
+    #print(grb)
     vals = grb.values
 
     lats, lons = grb.latlons()
 
+    for k in grb.keys():
+        print(k)
+
+    print(grb.analDate)
+    print(grb.validDate)
+    print(grb.units)
+    print(grb.typeOfLevel)
+    break
+
+    # TODO Data in xarray?
+
     print(grb.name)
-    #print(grb.level)
+    print(grb.level)
 
     #print('lats min/max: ', lats.shape, lats.max(), lats.min())
     #print('lons min/max: ', lons.shape, lons.max(), lons.min())
diff --git a/icon_download.sh b/icon_download.sh
index cf00cca..343b10c 100755
--- a/icon_download.sh
+++ b/icon_download.sh
@@ -9,9 +9,10 @@ MODEL=icon-d2
 MODEL_LONG=icon-d2_germany
 BASE="http://opendata.dwd.de/weather/nwp"
 
-RUN="15"
+RUN="00"
 PARAMETERS=( "t" "relhum" "u" "v" "fi" )
-PARAMETERS_SINGLE_LEVEL=( "cape_ml" "cin_ml" "tot_prec" "w_ctmax" )
+# tot_prec and cape_ml/cin_ml is in 15min intervals and screws with xygrib
+PARAMETERS_SINGLE_LEVEL=( "w_ctmax" )
 PRESSURE_LEVELS=( "1000" "975" "950" "850" "700" "600" "500" "400" "300" "250" "200" )
 OFFSETS=( "000" "003" "006" "009" "012" "015" "018" "024" )
 DATE=$(date +%Y%m%d)
@@ -20,24 +21,24 @@ mkdir -p $OUTDIR
 
 echo -n > "$OUTDIR/index.txt"
 
-for OFFSET in ${OFFSETS[@]}; do
-	for PARAMETER in ${PARAMETERS[@]}; do
-		for LEVEL in ${PRESSURE_LEVELS[@]}; do
+for OFFSET in "${OFFSETS[@]}"; do
+	for PARAMETER in "${PARAMETERS[@]}"; do
+		for LEVEL in "${PRESSURE_LEVELS[@]}"; do
 			URL="$BASE/$MODEL/grib/$RUN/$PARAMETER/${MODEL_LONG}_regular-lat-lon_pressure-level_${DATE}${RUN}_${OFFSET}_${LEVEL}_${PARAMETER}.grib2.bz2"
 			BNAME=$(basename "$URL")
 			echo Getting "$URL"
 			echo "${BNAME%.bz2}" >> $OUTDIR/index.txt
-			wget -q --directory-prefix=$OUTDIR  "$URL"
+			wget -q --directory-prefix=$OUTDIR  "$URL" || echo FAILED!
 
 		done
 	done
 
-	for PARAMETER in ${PARAMETERS_SINGLE_LEVEL[@]}; do
+	for PARAMETER in "${PARAMETERS_SINGLE_LEVEL[@]}"; do
 			URL="$BASE/$MODEL/grib/$RUN/$PARAMETER/${MODEL_LONG}_regular-lat-lon_single-level_${DATE}${RUN}_${OFFSET}_2d_${PARAMETER}.grib2.bz2"
 			BNAME=$(basename "$URL")
 			echo Getting "$URL"
 			echo "${BNAME%.bz2}" >> $OUTDIR/index.txt
-			wget -q --directory-prefix=$OUTDIR  "$URL"
+			wget -q --directory-prefix=$OUTDIR  "$URL" || echo FAILED!
 	done
 done
 
diff --git a/requirements.txt b/requirements.txt
index 7ce482f..aef6541 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,2 +1,4 @@
 metpy
 pygrib
+xarray
+cfgrib
diff --git a/skewt.py b/skewt.py
new file mode 100644
index 0000000..43bcc91
--- /dev/null
+++ b/skewt.py
@@ -0,0 +1,70 @@
+import matplotlib.gridspec as gridspec
+import matplotlib.pyplot as plt
+
+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
+
+class Skewt:
+    def __init__(self, p, T, Td, title=None):
+        self._p = p
+        self._T = T
+        self._Td = Td
+
+
+        # Create a new figure. The dimensions here give a good aspect ratio
+        self._fig = plt.figure(figsize=(9, 9))
+        plt.rcParams["font.family"] = "monospace"
+        #self._fig = plt.figure()
+
+        if title is not None:
+            plt.suptitle(title, x=0, y=0, va='bottom', ha='left')
+
+        # Grid for plots
+        self._gs = gridspec.GridSpec(3, 3)
+        self._skew = SkewT(self._fig, rotation=45, subplot=self._gs[:, :2])
+
+        # Plot the data using normal plotting functions, in this case using
+        # log scaling in Y, as dictated by the typical meteorological plot
+        self._skew.plot(p, T, 'r')
+        self._skew.plot(p, Td, 'b')
+
+        plt.xlabel('$T$ $[^o C]$')
+        plt.ylabel('$p$ $[hPa]$')
+
+    def addWindUV(self, u, v):
+        self._u = u
+        self._v = v
+
+        ax = self._fig.add_subplot(self._gs[0, -1])
+        h = Hodograph(ax, component_range=max(u + v).magnitude)
+        h.add_grid(increment=20)
+        h.plot_colormapped(u, v, self._p)
+        plt.tight_layout()
+        plt.xlabel('$m/s$')
+        plt.ylabel('$m/s$')
+        self._skew.plot_barbs(self._p, u, v)
+
+    def addInfo(self, lines):
+        # TODO
+        ax = self._fig.add_subplot(self._gs[1,-1])
+        ax.text(0, 0, lines, ha='left', va='center', size=12)
+        ax.axis("off")
+
+
+    def plot(self, filename=None):
+        # 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()
+
+        # Good bounds for aspect ratio
+        self._skew.ax.set_xlim(-30, 40)
+
+        if filename is not None:
+            plt.savefig(filename)
+        else:
+            plt.show()
diff --git a/xar.py b/xar.py
new file mode 100755
index 0000000..fca3887
--- /dev/null
+++ b/xar.py
@@ -0,0 +1,34 @@
+#!/usr/bin/env python3
+
+import xarray as xr
+from metpy.units import units
+import metpy.calc as mpcalc
+
+import skewt
+
+#grib = pygrib.open('dwd_icon-d2/combined.grib2')
+data = xr.load_dataset('dwd_icon-d2/combined.grib2', engine='cfgrib')
+
+lat = 47.9626
+lon = 11.9964
+
+for_temp = data.sel(latitude=lat, longitude = lon, method='nearest')
+
+for_temp = for_temp[['r', 't', 'u', 'v']]
+
+for step in for_temp.coords['step']:
+    this_step = for_temp.sel(step=step)
+
+    p = this_step.coords['isobaricInhPa'].values * units.hPa
+    T = this_step.t.values * units.K
+    relHum = this_step.r.values * units.percent
+    Td = mpcalc.dewpoint_from_relative_humidity(T, relHum)
+    u = this_step.u.values * (units.m / units.s)
+    v = this_step.v.values * (units.m / units.s)
+
+    skt = skewt.Skewt(p=p, T=T, Td=Td)
+    skt.addWindUV(u, v)
+    skt.addInfo("TEST")
+    skt.plot(filename=f'skewt.png')
+
+    break