make_era5_input.py

import os
import numpy as np
import pandas as pd
import xarray as xr


def make_era5(init_time, data_dir):

    init_time = pd.to_datetime(init_time)
    print(f"process {init_time} ...")

    pl_file = os.path.join(data_dir, init_time.strftime('P%Y%m%d%H.nc'))
    pl = xr.open_dataset(pl_file)

    sfc_file = os.path.join(data_dir, init_time.strftime('S%Y%m%d%H.nc'))
    sfc = xr.open_dataset(sfc_file)

    tp_file = os.path.join(data_dir, init_time.strftime('R%Y%m%d.nc'))
    tp = xr.open_dataarray(tp_file).fillna(0)
    tp = tp.rolling(time=6).sum() * 1000
    tp = tp.sel(time=tp.time[::6])
    tp = tp.clip(min=0, max=1000)
    sfc['tp'] = tp

    pl_names = ['z', 't', 'u', 'v', 'r']
    sfc_names = ['t2m', 'u10', 'v10', 'msl', 'tp']
    levels = [50, 100, 150, 200, 250, 300, 400, 500, 600, 700, 850, 925, 1000]

    channel = [f'{n.upper()}{l}' for n in pl_names for l in levels]
    channel +=[n.upper() for n in sfc_names]

    ds = []
    for name in pl_names + sfc_names:
        if name in ['z', 't', 'u', 'v', 'r']:
            v = pl[name]

        if name in ['t2m', 'u10', 'v10', 'msl', 'tp']:
            v = sfc[name]
            level = xr.DataArray([1], coords={'level': [1]}, dims=['level'])
            v = v.expand_dims({'level': level}, axis=1)             

        if np.isnan(v).sum() > 0:
            print(f"{name} has nan value")
            raise ValueError

        v.name = "data"
        v.attrs = {}                
        print(f"{name}: {v.shape}, {v.min().values} ~ {v.max().values}")
        ds.append(v)
     
    ds = xr.concat(ds, 'level')
    ds = ds.assign_coords(level=channel)
    ds = ds.rename({'longitude': 'lon', 'latitude': 'lat'})
    ds = ds.astype(np.float32)
    return ds
    
ds12 = make_era5('20230725-12', 'ERA520230725')
ds18 = make_era5('20230725-18', 'ERA520230725')
ds = xr.concat([ds12, ds18], 'time')
ds.to_netcdf('input.nc')