diff options
| -rw-r--r-- | Readme.md | 5 | ||||
| -rwxr-xr-x | grib.py | 23 | ||||
| -rwxr-xr-x | icon_download.sh | 17 | ||||
| -rw-r--r-- | requirements.txt | 2 | ||||
| -rw-r--r-- | skewt.py | 70 | ||||
| -rwxr-xr-x | xar.py | 34 |
6 files changed, 141 insertions, 10 deletions
@@ -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 @@ -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() @@ -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 |