diff options
Diffstat (limited to 'metchart/plotter')
| -rw-r--r-- | metchart/plotter/__init__.py | 0 | ||||
| -rw-r--r-- | metchart/plotter/debug_data.py | 4 | ||||
| -rwxr-xr-x | metchart/plotter/horizontal.py | 123 | ||||
| -rwxr-xr-x | metchart/plotter/meteogram.py | 164 | ||||
| -rwxr-xr-x | metchart/plotter/vertical_from_grib.py | 101 |
5 files changed, 392 insertions, 0 deletions
diff --git a/metchart/plotter/__init__.py b/metchart/plotter/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/metchart/plotter/__init__.py diff --git a/metchart/plotter/debug_data.py b/metchart/plotter/debug_data.py new file mode 100644 index 0000000..560d13c --- /dev/null +++ b/metchart/plotter/debug_data.py @@ -0,0 +1,4 @@ +def run(data, **kwargs): + print(data) + + return [] diff --git a/metchart/plotter/horizontal.py b/metchart/plotter/horizontal.py new file mode 100755 index 0000000..89945c3 --- /dev/null +++ b/metchart/plotter/horizontal.py @@ -0,0 +1,123 @@ +#!/usr/bin/env python3 +import os +import json + +import xarray as xr + +import numpy as np +import matplotlib.pyplot as plt +from metpy.plots import MapPanel, PanelContainer, RasterPlot, ContourPlot + +from . 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(data, plots, output='.'): + misc.create_output_dir(output) + index = [] + + for plot in plots: + index.append(_plot(data, output, **plot)) + + return index + +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}' + if '_description' in data.attrs: + panel.right_title = data.attrs['_description'] + + pc = PanelContainer() + pc.size = (12.8, 9.6) + #pc.figure.layout='constrained' + 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)) + plt.close('all') + + index.append( + { + 'file': outname, + 'init': init_str, + 'valid': valid_str, + 'valid_offset': hours_since_init_str, + 'display_name': name, + 'id': name + } + ) + + with open(os.path.join(output, f'{name}.index.json'), 'w') as f: + f.write(json.dumps(index, indent=4)) + + return { 'name': name, 'indexfile': f'{name}.index.json', 'list_title': 'INIT+' } + +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/metchart/plotter/meteogram.py b/metchart/plotter/meteogram.py new file mode 100755 index 0000000..8335247 --- /dev/null +++ b/metchart/plotter/meteogram.py @@ -0,0 +1,164 @@ +#!/usr/bin/env python3 +import os +import json + +import numpy as np + +import matplotlib.pyplot as plt + +import metpy.calc as mpcalc + +from .. import misc + +HEIGHT = 13 + +def run(data, plots, output='.', name='meteogram'): + misc.create_output_dir(output) + index = [] + + for plot in plots: + index.append(_plot(data, output, **plot)) + + with open(os.path.join(output, f'{name}.index.json'), 'w') as f: + f.write(json.dumps(index, indent=4)) + + return [{ 'name': name, 'indexfile': f'{name}.index.json', 'list_title': 'Location' }] + +def _get_next_subplot(size, counter=0): + ret = (counter + 1, counter + size) + counter += size + return counter, ret + +def _add_cloudcov(ax, data): + ax.set_ylabel('Pressure level [hPa]') + + clc = ax.contourf(data.valid_time, data.isobaricInhPa, data.ccl.transpose(), cmap='clcov', vmin=0, vmax=100, levels=9) + + # use Format parameter for n/8 + plt.colorbar(clc, label='cloudcov', extendfrac=None, ticks=[100*n/8 for n in range(9)], format=lambda x,_: f'{int(x/12.5)}/8', pad=0.0, fraction=0.015) + + cf = ax.contour(data.valid_time, data.isobaricInhPa, data.t.metpy.convert_units('degC').transpose()) + ax.clabel(cf, inline=True, fontsize=10) + + ax.barbs( + data.valid_time, + data.isobaricInhPa, + data.u.metpy.convert_units('kt').transpose(), + data.v.metpy.convert_units('kt').transpose() + ) + + ax.invert_yaxis() + +def _add_temp_dewpoint(ax, data): + ### Temp + Dewpoint + ax.plot(data.valid_time, data.t2m.metpy.convert_units('degC').transpose(), color='red', label='Temperature (2m)') + ax.plot(data.valid_time, mpcalc.dewpoint_from_relative_humidity(data.t2m, data.r2).transpose(), color='blue', label='Dewpoint (2m)') + ax.plot(data.valid_time, data.sel(isobaricInhPa=850.0).t.metpy.convert_units('degC').transpose(), color='grey', label='Temperature (850hPa)') + ax.set_ylabel('Temperature [degC]') + ax.legend(loc='lower right') + +def _add_mslp(ax, data): + ax.plot(data.valid_time, data.prmsl.metpy.convert_units('hPa').transpose(), color='black', label='Temperature (2m)') + ax.set_ylabel('Mean Sea Level Pressure [hPa]') + +def _add_convective_clouds(ax, data): + # TODO: ADD HBAS_CON, HTOP_CON + # If none: -500m + ax.set_ylim(0, 14) + ax.set_ylabel('Convective Clouds Height [km]') + ax.bar(data.valid_time, alpha=0.5, + bottom=data.HBAS_CON.metpy.convert_units('km').transpose(), + height=(data.HTOP_CON.metpy.convert_units('km')-data.HBAS_CON.metpy.convert_units('km')).transpose(), + align='edge', width=np.timedelta64(3, 'h')) + +def _add_precip(ax, data): + ax.set_ylabel('Total precipitation [mm]') + ax.set_ylim(0, 30) + ax.bar(data.valid_time[:-1], data.tp.diff('step').transpose(), width=np.timedelta64(3, 'h'), + align='edge', alpha=0.7, color='green') + + ax_p = ax.twinx() + ax_p.set_ylabel('Snow depth [m]') + ax_p.set_ylim(bottom=0) + ax_p.plot(data.valid_time, data.sde.transpose(), color='blue') + +def _add_surface_wind(ax, data): + ax.plot(data.valid_time, mpcalc.wind_speed(data.u10.transpose(), data.v10.transpose()), color='black', label='Wind (10m)') + ax.plot(data.valid_time, data.fg10.transpose(), color='red', label='Gust (10m)') + + ax_b = ax.twinx() + ax_b.barbs( + data.valid_time, + [1 for _ in data.valid_time], + data.u10.metpy.convert_units('kt').transpose(), + data.v10.metpy.convert_units('kt').transpose() + ) + ax_b.axis('off') + + ax.set_ylabel('Wind Speed [m/s]') + ax.legend(loc='lower right') + +def _plot(data, output, name, lat, lon): + data = data.sel(latitude=lat, longitude = lon, method='nearest') + + fig = plt.figure(figsize=(12, 12), layout="constrained") + + sp_cnt, spec = _get_next_subplot(4) + ax = fig.add_subplot(HEIGHT,1,spec) + _add_cloudcov(ax, data) + + sp_cnt, spec2 = _get_next_subplot(2,sp_cnt) + ax2 = fig.add_subplot(HEIGHT,1,spec2,sharex=ax) + _add_temp_dewpoint(ax2, data) + + sp_cnt, spec3 = _get_next_subplot(2,sp_cnt) + ax3 = fig.add_subplot(HEIGHT,1,spec3,sharex=ax) + #ax3.legend(loc='lower right') + _add_mslp(ax3, data) + + ax4 = ax3.twinx() + _add_convective_clouds(ax4, data) + + sp_cnt, spec4 = _get_next_subplot(2,sp_cnt) + ax5 = fig.add_subplot(HEIGHT,1,spec4,sharex=ax) + _add_precip(ax5, data) + + sp_cnt, spec5 = _get_next_subplot(2,sp_cnt) + ax6 = fig.add_subplot(HEIGHT,1,spec5,sharex=ax) + _add_surface_wind(ax6, data) + + ### Info Lines + sp_cnt, spec5 = _get_next_subplot(1,sp_cnt) + ax_text = fig.add_subplot(HEIGHT, 1, spec5) + + info_lines = [] + init = misc.np_time_convert(data.time.values) + init_str = init.strftime('%d %b %Y - %HUTC') + init_for_filename = init.strftime('%Y-%m-%d-%HUTC') + + info_lines.append(f'{name}') + info_lines.append(f"INIT : {init_str}") + info_lines.append(f"LAT {lat} LON {lon}") + + if '_description' in data.attrs: + info_lines.append(data.attrs['_description']) + + ax_text.text(0, 0, '\n'.join(info_lines), ha='left', va='center', + size=10, fontfamily='monospace') + ax_text.axis("off") + + ### Output + + outname = f'{name}_{init_for_filename}.png' + plt.savefig(os.path.join(output, outname)) + plt.close('all') + + return ( + { + 'file': outname, + 'init': init_str, + 'valid': init_str, + 'valid_offset': '00', + 'display_name': name, + 'id': name + }) diff --git a/metchart/plotter/vertical_from_grib.py b/metchart/plotter/vertical_from_grib.py new file mode 100755 index 0000000..c632b0b --- /dev/null +++ b/metchart/plotter/vertical_from_grib.py @@ -0,0 +1,101 @@ +#!/usr/bin/env python3 +import os + +import datetime +import json + +import matplotlib.pyplot as plt + +import xarray as xr +from metpy.units import units +import metpy.calc as mpcalc +import numpy as np + +import skewt + +from . import misc + +config = { + 'source':'dwd_icon-eu/combined.grib2', + 'plots':[ + { + 'lat':47.9626, + 'lon':11.9964, + 'name':'Antersberg', + 'analysis':'lcl' + }, + ] +} + +def run(data, plots, output='.'): + misc.create_output_dir(output) + index = [] + + for plot in plots: + index.append(_plot(data, output, **plot)) + + return index + +def _plot(data, output, name, lat=None, lon=None, analysis=None): + for_temp = data[['r', 't', 'u', 'v']] + + if not (lat is None and lon is None): + for_temp = for_temp.sel(latitude=lat, longitude=lon, method='nearest') + + 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) + + if '_description' in data.attrs: + skt.addInfo(data.attrs['_description']) + + 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)) + + plt.close('all') + + index.append( + { + 'file': outname, + 'init': init_str, + 'valid': valid_str, + 'valid_offset': hours_since_init_str, + 'display_name': name, + 'id': name + } + ) + + with open(os.path.join(output, f'skewt_{name}.index.json'), 'w') as f: + f.write(json.dumps(index, indent=4)) + + return {'name': name, 'indexfile': f'skewt_{name}.index.json', 'list_title': 'Location'} + +if __name__ == '__main__': + run(**config) |