概要
気象分野には，渦度という物理量があります．
渦度とは，簡単に言うと「局所1領域における流れの回転の強さを表す量」です．
これは風のデータから簡単に計算することができるので，Pythonを使って渦度を計算してみよう，というのが本記事の内容です．
渦度とは
流体力学では渦度ベクトル$\boldsymbol{\omega}$を，流速ベクトル$\boldsymbol{v}=(u,v,w)$の回転で定義します．
\boldsymbol{\omega} \equiv \text{rot}\boldsymbol{v}
= \left(
\frac{\partial w}{\partial y} - \frac{\partial v}{\partial z},
\frac{\partial u}{\partial z} - \frac{\partial w}{\partial x},
\frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}
\right)
一方，気象学で「渦度」と言ったときには，鉛直渦度，すなわち渦度ベクトルの鉛直成分
\frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}\tag{1}
を指すことが多いです．
今回は，この鉛直渦度を計算していきます．
これ以降，鉛直渦度のことを渦度と呼ぶことにします．
計算方法
渦度は式$(1)$のとおりなのですが，偏微分の形をしているため，このままでは具体的に計算することができません．
そこで，式$(1)$を離散化します．
\frac{v_2-v_1}{\Delta x} - \frac{u_2-u_1}{\Delta y}\tag{2}
ここで，$v_2-v_1$は$x$軸方向（東西方向）の$v$の変化量，$u_2-u_1$は$y$軸方向（南北方向）の$u$の変化量です．
こちらの「色と形で気象予報士！」さんのブログ記事2にあった図がとてもわかりやすいかと思います．
あとは，離散化した式$(2)$にデータを代入して計算するだけです．
Pythonによる実装
風のデータは，大気再解析データから取ってきました．
- データ元：京都大学生存圏研究所のサイト内にある，「解析値を中心に再構成したNetCDFデータ」からダウンロード
- ファイル場所：MSM-P/2023/1205.nc
1205.ncは100MB以上あるので，必要な風データのみを抽出して1205_uv.ncとして保存しました
使用ライブラリは以下です．
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import cartopy.crs as ccrs
データの読み込み
再解析データはNetCDFという形式のファイルなので，xarrayというライブラリを使って読み込みます．
data = xr.open_mfdataset('./1205_uv.nc').load()
風の東西成分と南北成分を，それぞれ変数uとvに格納します．
今回は，2023年12月5日9時（UTC）の500hPa面の渦度を計算することにします．
TIME = '2023-12-05T00:00:00'
P = 500
u = data['u'].sel(time=TIME, p=P).values
v = data['v'].sel(time=TIME, p=P).values
あと，緯度経度座標も取得しておきます．
式$(2)$の$\Delta x, \Delta y$の計算に必要だからです3．
lon, lat = np.meshgrid(data['lon'].values, data['lat'].values)
ある2地点の緯度経度が与えられたときの距離の算出方法は，以前私が書いた記事を参考にしてください．
渦度の計算
あとは，式$(2)$に値をぶち込んでいくだけです．
\frac{v_2-v_1}{\Delta x} - \frac{u_2-u_1}{\Delta y}\tag{2}
なお，愚直に二重forループで実装しました4．
zeta = np.zeros_like(u)
for i in range(1,lat.shape[0]-1):
for j in range(1,lat.shape[1]-1):
u1 = u[i+1][j]
u2 = u[i-1][j]
v1 = v[i][j-1]
v2 = v[i][j+1]
dx = dist((lat[i][j+1], lon[i][j+1]), (lat[i][j-1], lon[i][j-1]))
dy = dist((lat[i+1][j], lon[i+1][j]), (lat[i-1][j], lon[i-1][j]))
zeta[i][j] = (v2-v1)/(2*dx) - (u2-u1)/(2*dy)
結果の可視化
では，計算結果を可視化してみましょう．
地図の描画方法については，以前記事を書いていますので参考にしてください．
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([lon.min(), lon.max(), lat.min(), lat.max()], crs=ccrs.PlateCarree())
ax.coastlines()
ax.gridlines(
draw_labels=['left','bottom'],
xlocs=np.arange(120, 150, 2.5),
ylocs=np.arange(20, 50, 2.5),
)
plt.imshow(
zeta*1e6,
cmap='bwr', vmax=600, vmin=-600,
extent=[lon.min(), lon.max(), lat.min(), lat.max()],
transform=ccrs.PlateCarree())
cbar = plt.colorbar()
cbar.set_label(f'Vorticity at {P}hPa ($10^{-6}$ /s)')
plt.title(TIME + '(UTC)')
plt.tight_layout()
plt.savefig('vorticity.png', dpi=300)
plt.show()
いい感じですね．
念のため，答え合わせをしてみましょう．
下図はウェザーニューズ社提供のLabs Ch.から取ってきた，同時刻同圧力面における渦度（など）の値を示しています．
解像度は違いますが，だいたい一致していますね．
黄海にのびる500hPaトラフに対応する高渦度領域を見てもらえるとわかりやすいかと思います．
ソースコード
Qiitaアドベントカレンダーで使ったコードはこちらのレポジトリにまとめています．
本記事のコードは07-vorticityの中に入っています．