diff options
Diffstat (limited to 'plotter')
| -rw-r--r-- | plotter/__init__.py | 0 | ||||
| -rwxr-xr-x | plotter/horizontal.py | 113 | ||||
| -rwxr-xr-x | plotter/vertical_from_grib.py | 88 |
3 files changed, 201 insertions, 0 deletions
diff --git a/plotter/__init__.py b/plotter/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/plotter/__init__.py diff --git a/plotter/horizontal.py b/plotter/horizontal.py new file mode 100755 index 0000000..a81b372 --- /dev/null +++ b/plotter/horizontal.py @@ -0,0 +1,113 @@ +#!/usr/bin/env python3 +import os +import json + +import xarray as xr + +import numpy as np +from metpy.plots import MapPanel, PanelContainer, RasterPlot, ContourPlot + +import misc + +config = { + 'source': 'dwd_icon-eu/combined.grib2', + 'plots': [ + { + 'name':'r_t-750', + 'area': None, + 'layers': [ + { + 'layertype': 'raster', + 'field': 'r', + 'level': 750, + }, + { + 'layertype': 'contour', + 'field': 't', + 'level': 750, + 'contours': 5, + 'clabels': True + }, + ] + }, + ] +} + +def run(source, plots, output='.'): + misc.create_output_dir(output) + data = xr.load_dataset(source, engine='cfgrib') + + for plot in plots: + _plot(data, output, **plot) + +def _plot(data, output, name, layers, area = None): + index = [] + + for step in data.coords['step']: + this_step = data.sel(step=step) + + map_layers = [] + + for layer in layers: + map_layers.append(_layer(this_step, **layer)) + + valid = misc.np_time_convert(step.valid_time.values) + init = misc.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) + init_for_filename = init.strftime('%Y-%m-%d-%HUTC') + + panel = MapPanel() + if area is not None: + panel.area = area + panel.projection = 'mer' + panel.layers = ['coastline', 'borders'] + panel.plots = map_layers + panel.left_title = f'{name} VALID: {valid_str} (INIT +{hours_since_init_str}) INIT: {init_str}' + panel.right_title = 'FORECAST DWD ICON-EU' + + pc = PanelContainer() + pc.size = (12.8, 9.6) + pc.panels = [panel] + pc.draw() + #pc.show() + outname = f'{name}_{init_for_filename}+{hours_since_init_str}.png' + pc.save(os.path.join(output, outname)) + + index.append( + { + 'file': outname, + 'init': init_str, + 'valid': valid_str, + 'valid_offset': hours_since_init_str + } + ) + + with open(os.path.join(output, f'{name}.index.json'), 'w') as f: + f.write(json.dumps(index, indent=4)) + +def _layer(data, layertype, **kwargs): + layertypes={ + 'raster': { + 'obj': RasterPlot, + 'defaults': { + 'colorbar': 'vertical', + } + }, + 'contour': { + 'obj': ContourPlot, + 'defaults': {} + } + } + + args = layertypes[layertype]['defaults'] | kwargs + + ret = layertypes[layertype]['obj'](**args) + ret.data = data + + return ret + +if __name__ == '__main__': + run(**config) diff --git a/plotter/vertical_from_grib.py b/plotter/vertical_from_grib.py new file mode 100755 index 0000000..8aa851a --- /dev/null +++ b/plotter/vertical_from_grib.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python3 +import os + +import datetime +import json + +import xarray as xr +from metpy.units import units +import metpy.calc as mpcalc +import numpy as np + +import skewt + +import misc + +config = { + 'source':'dwd_icon-eu/combined.grib2', + 'plots':[ + { + 'lat':47.9626, + 'lon':11.9964, + 'name':'Antersberg', + 'analysis':'lcl' + }, + ] +} + +def run(source, plots, output='.'): + misc.create_output_dir(output) + data = xr.load_dataset(source, engine='cfgrib') + + for plot in plots: + _plot(data, output, **plot) + +def _plot(data, output, lat, lon, name, analysis=None): + for_temp = data.sel(latitude=lat, longitude = lon, method='nearest') + for_temp = for_temp[['r', 't', 'u', 'v']] + index = [] + + 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) + + valid = misc.np_time_convert(step.valid_time.values) + init = misc.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(f'{name} INIT+' + hours_since_init_str) + skt.addInfo(f"VALID: {valid_str}") + skt.addInfo(f"INIT : {init_str}") + skt.addInfo(f"LAT {lat} LON {lon}") + + if analysis is not None: + skt.addAnalysis(shade=True, analysis=analysis) + + # TODO get from source! + skt.addInfo("FORECAST DWD ICON-D2") + + init_for_filename = init.strftime('%Y-%m-%d-%HUTC') + + outname = f'skewt_{name}_{init_for_filename}+{hours_since_init_str}.png' + skt.plot(filename=os.path.join(output, outname)) + + index.append( + { + 'file': outname, + 'init': init_str, + 'valid': valid_str, + 'valid_offset': hours_since_init_str + } + ) + + with open(os.path.join(output, f'skewt_{name}.index.json'), 'w') as f: + f.write(json.dumps(index, indent=4)) + +if __name__ == '__main__': + run(**config) |