aboutsummaryrefslogtreecommitdiff
path: root/xar.py
blob: 759bca44b2fea08a41b1c56931c3b2336878fa8d (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
#!/usr/bin/env python3

import datetime

import xarray as xr
from metpy.units import units
import metpy.calc as mpcalc
import numpy as np

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)

config = {
    'source':'dwd_icon-d2/combined.grib2',
    'plots':[
        {
            'lat':47.9626,
            'lon':11.9964,
            'name':'Antersberg'
        },
    ]
}

def run(config):
    data = xr.load_dataset(config['source'], engine='cfgrib')

    for plot in config['plots']:
        _plot(data, **plot)

def _plot(data, lat, lon, name):
    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)

        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(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}")
        skt.addAnalysis(shade=True, analysis='lcl')
        skt.addInfo("FORECAST DWD ICON-D2")

        init_for_filename = init.strftime('%Y-%m-%d-%HUTC')

        skt.plot(filename=f'skewt_{name}_{init_for_filename}+{hours_since_init_str}.png')

if __name__ == '__main__':
    run(config)