FuXi/make_gfs_input.py at main · Conyboy/FuXi · GitHub

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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
import os
import numpy as np
import pandas as pd
import pygrib as pg
import xarray as xr

def make_gfs(src_name):
    assert os.path.exists(src_name)

    levels = [50, 100, 150, 200, 250, 300, 400, 500, 600, 700, 850, 925, 1000]
    pl_names = ['gh', 't', 'u', 'v', 'r']
    sf_names = ['2t', '10u', '10v', 'mslet']

    try:
        ds = pg.open(src_name)
    except:
        print(f"{src_name} not found")
        return

    input = []
    level = []

    for name in pl_names + sf_names + ["tp"]:

        if name in pl_names:
            try:
                data = ds.select(shortName=name, level=levels)
            except:
                print("pl wrong")
                return

            data = data[:len(levels)]

            if len(data) != len(levels):
                print("pl wrong")
                return

            if name == "gh":
                name = "z"

            for v in data:
                init_time = f'{v.date}-{v.time//100:02d}'

                lat = v.distinctLatitudes
                lon = v.distinctLongitudes

                img, _, _ = v.data()

                if name == "z":
                    img = img * 9.8

                input.append(img)
                level.append(f'{name}{v.level}')
                print(f"{v.name}: {v.level}, {img.shape}, {img.min()} ~ {img.max()}")

        if name in sf_names:
            try:
                data = ds.select(shortName=name)
            except:
                print('sfc wrong')
                return

            name_map = {'2t': 't2m', '10u': 'u10', '10v': 'v10', 'mslet': 'msl'}
            name = name_map[name]

            for v in data:
                img, _, _ = v.data()
                input.append(img)
                level.append(name)
                print(f"{v.name}: {img.shape}, {img.min()} ~ {img.max()}")

        if name == "tp":
            tp = img * 0
            input.append(tp)
            level.append("tp")

    input = np.stack(input)
    assert input.shape[-3:] == (70, 721, 1440)
    assert input.max() < 1e10

    times = [pd.to_datetime(init_time)]
    input = xr.DataArray(
        data=input[None],
        dims=['time', 'level', 'lat', 'lon'],
        coords={'time': times, 'level': level, 'lat': lat, 'lon': lon},
    )

    if np.isnan(input).sum() > 0:
        print("Field has nan value")
        return

    return input


def test_make_gfs():
    d1 = make_gfs('30/gfs.t06z.pgrb2.0p25.f000')
    d2 = make_gfs('30/gfs.t12z.pgrb2.0p25.f000')

    if d1 and d2:
        ds = xr.concat([d1, d2], 'time')
        ds = ds.assign_coords(time=ds.time.astype(np.datetime64))
        ds.to_netcdf('input.nc')