数値予報データ

数値予報データ#

気象庁の数値予報データを描いてみましょう。気象庁メソモデルのデータを使います。まず，京都大学生存圏研究所から提供されている，解析値を中心に再構成し，扱いやすいNetCDF形式に変換したデータで3時間毎の解析値を描きます。ここでは例として2022年1月6日のデータを使います。curlでデータをダウンロードして，NetCDFファイルをxarrayで読みます。

import numpy as np
import xarray as xr

!curl -O http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/netcdf/MSM-P/2022/0106.nc

  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100  178M  100  178M    0     0  77.0M      0  0:00:02  0:00:02 --:--:-- 77.0M

ds = xr.open_dataset("0106.nc")

変数名などを確認するために，メタデータを表示します。

3時間毎に8時刻分のデータが格納されていることが分かります。1日分を1行に4時刻ずつ並べて描くことにします。850 hPa面で気温をカラーの陰影で，風をベクトルで描きます。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs


plt.rcParams["font.size"] = 18
fig = plt.figure(figsize=[20, 8], constrained_layout=True)

skip = 12
plev = 850
qx, qy, qk = 1.2, 0.05, 20
axs = []
tlevs = np.arange(252, 289, 3)
for i in range(ds.time.size):
    ax = fig.add_subplot(2, 4, i+1, projection=ccrs.PlateCarree())
    T = ds.temp.sel(p=plev).isel(time=i)
    u = ds.u.sel(p=plev).isel(time=i)[::skip, ::skip]
    v = ds.v.sel(p=plev).isel(time=i)[::skip, ::skip]
    p = ax.contourf(ds.lon, ds.lat, T,
                    tlevs, transform=ccrs.PlateCarree(), cmap="coolwarm")
    q = ax.quiver(ds.lon[::skip], ds.lat[::skip], u, v, scale=400, transform=ccrs.PlateCarree())
    if i+1==ds.time.size:
        ax.quiverkey(q, qx, qy, qk, f"{qk} m/s")
    ax.set_title(ds.time[i].dt.strftime("%Y-%m-%d %HUTC").values)
    ax.coastlines()
    ax.gridlines()
    axs.append(ax)
fig.colorbar(p, ax=axs, shrink=0.8, label="T K")
fig.suptitle(f"MSM Analysis {plev}hPa")
plt.show()

_images/bf2338ba2dfe9ef090e683f331f8feb5ba640ceb65dec4ececa568008416f213.png

5行目: フォントサイズを一括して18ポイントにします。
6行目: プロットの間の間隔を詰めるためにconstrained_layout=Trueとしています。
8行目: ベクトルは22行目でquiver()を使って描きますが，そのままだとデータが多すぎるので12点毎にひとつ描くようにします。
11行目: colorbar()で使うaxを格納する空のリストを作ります。
12行目: 気温は3K毎に色を変えます。
15〜17行目: データセットから気温と風を変数に格納します。sel()で等圧面を気圧の値で，isel()で時刻を番号で選択します。
16, 17行目: skipを使って風を水平方向に間引いています。
18〜19行目: 気温を陰影で描きます。colorbar()に使うために，結果をpに格納します。
20行目: 風をベクトルで描きます。quiverkey()に使うために，結果をqに格納します。
23行目: タイトルとして時刻を設定します。
26行目: axをリストaxsに追加します。
27行目: カラーバーを描きます。
28行目: 全体のタイトルとしてモデル名と気圧面を設定します。

南岸低気圧が東進して寒気が日本に吹き出しています。予報データと比較してみましょう。最新の予報データはMSM-P.ncなどにリンクされています。また，最近1週間のデータが提供されています。執筆時点取得できた2022年1月6日00 UTCからの予報データを用います。

!curl -O http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/latest/MSM-P.nc

  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100  178M  100  178M    0     0  84.1M      0  0:00:02  0:00:02 --:--:-- 84.1M

!curl -O http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/latest/20220106/MSM2022010600P.nc

  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100  178M  100  178M    0     0  75.9M      0  0:00:02  0:00:02 --:--:-- 75.9M

ds1 = xr.open_dataset("MSM2022010600P.nc")

import matplotlib.pyplot as plt
import cartopy.crs as ccrs


plt.rcParams["font.size"] = 18
fig = plt.figure(figsize=[20, 12], constrained_layout=True)

skip = 12
plev = 850
qx, qy, qk = 1.2, 0.05, 20
axs = []
tmax = ds1.temp.sel(p=plev).max()
tlevs = np.arange(252, 289, 3)
for i in range(ds1.time.size):
    T = ds1.temp.sel(p=plev).isel(time=i)
    u = ds1.u.sel(p=plev).isel(time=i)[::skip, ::skip]
    v = ds1.v.sel(p=plev).isel(time=i)[::skip, ::skip]
    p = ax.contourf(ds.lon, ds.lat, T,
                    tlevs, transform=ccrs.PlateCarree(), cmap="coolwarm")
    q = ax.quiver(ds.lon[::skip], ds.lat[::skip], u, v, scale=400, transform=ccrs.PlateCarree())
    if i+1==ds1.time.size:
        ax.quiverkey(q, qx, qy, qk, f"{qk} m/s")
    ax.set_title(ds1.time[i].dt.strftime("%Y-%m-%d %HUTC").values)
    ax.coastlines()
    ax.gridlines()
    axs.append(ax)
fig.colorbar(p, ax=axs, shrink=0.8, label="T K")
fig.suptitle(f"MSM  Forecast {plev}hPa")
plt.show()

_images/8e181959bc28daf3dd2e384eabade6841c8785acc6e6abda9bbf9291c4c3940b.png

6行目: 3行になるので図の高さを大きくします。
28行目: 全体のタイトルのAnalysisをForecastにしました。

予報ではやや寒気の吹き出しが弱い傾向がありますが，南岸低気圧とその後の冬型の強化が精度良く予報されています。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs


plt.rcParams["font.size"] = 18
fig = plt.figure(figsize=[20, 8], constrained_layout=True)

skip = 12
plev = 850
qx, qy, qk = 1.2, -0.02, 5
for i in range(1, ds.time.size):
    ax = fig.add_subplot(2, 4, i+1, projection=ccrs.PlateCarree())
    T = ds.temp.sel(p=plev).isel(time=i)
    dT = ds1.temp.sel(p=plev).isel(time=i) - T
    u = ds.u.sel(p=plev).isel(time=i)[::skip, ::skip]
    v = ds.v.sel(p=plev).isel(time=i)[::skip, ::skip]
    du = ds1.u.sel(p=plev).isel(time=i)[::skip, ::skip] - u
    dv = ds1.v.sel(p=plev).isel(time=i)[::skip, ::skip] - v
    p = ax.contourf(ds.lon, ds.lat, dT, transform=ccrs.PlateCarree(), cmap="coolwarm", vmin=dT.min(), vmax=-dT.min())
    q = ax.quiver(ds.lon[::skip], ds.lat[::skip], du, dv, transform=ccrs.PlateCarree())
    ax.quiverkey(q, qx, qy, qk, f"{qk} m/s")
    ax.set_title(ds.time[i].dt.strftime("%Y-%m-%d %HUTC").values)
    ax.coastlines()
    ax.gridlines()
    fig.colorbar(p, ax=ax, shrink=0.6, label="T K")
    fig.suptitle(f"MSM Forecast-Analysis {plev}hPa")
plt.show()

_images/690f35ef2275415661f60b8f72ce3668c7d5a403ea8bee26d4186f5b3c7ca680.png