GSMの500hPa面天気図作成 3
今回は、天気図上に低気圧などのマークを表示する方法について紹介します。コードの説明は、図上の座標取得、GPVから極値とその位置の取得、マークの描画の3点となります。
これで気象庁が作成している500hPa面天気図と同じ内容のカラー版が作成できます。この記事の最後に、Jupyter Notebook用のコードをダウンロードできるようにしています。
図上の座標を取得する
緯度経度で指定したポイントを図の座標系へ変換する関数のコードは次の通りとなります。この関数を使えば、文字列やマークを指定した位置に表示することができます。詳細はコード内のコメントにて説明しています。参考にしてください。
## 緯度経度で指定したポイントの図上の座標などを取得する関数 transform_lonlat_to_figure()
# 図法の座標 => pixel座標 => 図の座標 と3回の変換を行う
# pixel座標: plt.figureで指定した大きさxDPIに合わせ、左下を原点とするpixelで測った座標
# 図の座標: axesで指定した範囲を(0,1)x(0,1)とする座標
# 3つの座標を出力する
# 図の座標, Pixel座標, 図法の座標
def transform_lonlat_to_figure(lonlat, ax, proj):
# lonlat:経度と緯度 (lon, lat)
# ax: Axes図の座標系 ex. fig.add_subplot()の戻り値
# proj: axで指定した図法
#
# 例 緯度経度をpointで与え、ステレオ図法る場合
# point = (140.0,35.0)
# proj= ccrs.Stereographic(central_latitude=60, central_longitude=140)
# fig = plt.figure(figsize=(20,16))
# ax = fig.add_subplot(1, 1, 1, projection=proj)
# ax.set_extent([108, 156, 17, 55], ccrs.PlateCarree())
#
## 図法の変換
# 参照 https://scitools.org.uk/cartopy/docs/v0.14/crs/index.html
point_proj = proj.transform_point(*lonlat, ccrs.PlateCarree())
#
# pixel座標へ変換
# 参照 https://matplotlib.org/stable/tutorials/advanced/transforms_tutorial.html
point_pix = ax.transData.transform(point_proj)
#
# 図の座標へ変換
point_fig = ax.transAxes.inverted().transform(point_pix)
return point_fig, point_pix, point_proj
GPVから低気圧などの極値の位置を求める
GPVから極値の位置を求め際には、注意が必要です。単純に周囲8格子の値より大きい(または小さい)値を選ぶと、多くの極値が近接して存在することがあるため、間引く必要があります。
まずは(n+1)格子 x (n+1)格子 の矩形真ん中の格子が、これら格子点の中で最大値(または最小値)となっている格子をサーチします。こうして得られた複数の格子において、それぞれの格子点間の距離を求め、指定した距離より短い格子においては、どちらか値が大きい方(または小さい方)のみ極値として検出することにします。値が同じ場合は、極値間の距離の合計が最も小さいもの(極値が密集している位置と考えられる)を選択しています。
コードは次の通りです。(まだまだ改善の余地があります。いい改善方法があれば、アドバイスください。)
## 極大/極小ピーク検出関数
def detect_peaks(image, filter_size=3, dist_cut=5.0, flag=0):
# filter_size: この値xこの値 の範囲内の最大値のピークを検出
# dist_cut: この距離内のピークは1つにまとめる
# flag: 0:maximum検出 0以外:minimum検出
if flag==0:
local_max = maximum_filter(image,
footprint=np.ones((filter_size, filter_size)), mode='constant')
detected_peaks = np.ma.array(image, mask=~(image == local_max))
else:
local_min = minimum_filter(image,
footprint=np.ones((filter_size, filter_size)), mode='constant')
detected_peaks = np.ma.array(image, mask=~(image == local_min))
peaks_index = np.where((detected_peaks.mask != True))
# peak間の距離行例を求める
(x,y) = peaks_index
size=y.size
dist=np.full((y.size, y.size), -1.0)
for i in range(size):
for j in range(size):
if i == j:
dist[i][j]=0.0
elif i>j:
d = math.sqrt(((y[i] - y[j])*(y[i] - y[j]))
+((x[i] - x[j])*(x[i] - x[j])))
dist[i][j]= d
dist[j][i]= d
# 距離がdist_cut内のpeaksの距離の和と、そのピーク番号を取得する
Kinrin=[]
dSum=[]
for i in range(size):
tmpA=[]
distSum=0.0
for j in range(size):
if dist[i][j] < dist_cut and dist[i][j] > 0.0:
tmpA.append(j)
distSum=distSum+dist[i][j]
dSum.append(distSum)
Kinrin.append(tmpA)
# Peakから外すPeak番号を求める. peak間の距離和が最も小さいものを残す
cutPoint=[]
for i in range(size):
val = dSum[i]
val_i=image[x[i]][y[i]]
for k in Kinrin[i]:
val_k=image[x[k]][y[k]]
if flag==0 and val_i < val_k:
cutPoint.append(i)
break
if flag!=0 and val_i > val_k:
cutPoint.append(i)
break
if val > dSum[k]:
cutPoint.append(i)
break
if val == dSum[k] and i > k:
cutPoint.append(i)
break
# 戻り値用に外すpeak番号を配列から削除
newx=[]
newy=[]
for i in range(size):
if (i in cutPoint):
continue
newx.append(x[i])
newy.append(y[i])
peaks_index=(np.array(newx),np.array(newy))
return peaks_index
天気図上にマークを描画
天気図上の渦度極大域に+のマークや数値を描画するコードを次に示します。+マークなどのマーカーはplot()関数を利用して、渦度の極値の値や"L"などの文字はtext()関数を利用して描画させます。
## 正渦度中心に+マーク & 渦度の値描画
# 極値の位置を、ある程度間引いて取得
maxid = detect_peaks(ds['vorticity'].values, filter_size=3, dist_cut=4.0)
for i in range(len(maxid[0])):
wlon = ds['lon'][maxid[1][i]]
wlat = ds['lat'][maxid[0][i]]
# 図の範囲内に座標があるか確認
fig_z, _, _ = transform_lonlat_to_figure((wlon,wlat),ax,proj)
if ( fig_z[0] > 0.0 and fig_z[0] < 1.0 and fig_z[1] > 0.0 and fig_z[1] < 1.0):
val = ds['vorticity'].values[maxid[0][i]][maxid[1][i]]
# 表示する数値
ival = int(val * 1000000.0)
if ival > 30:
# マークは、緯度経度で
ax.plot(wlon, wlat, marker='+' , markersize=7, color="red", transform=latlon_proj)
if ival > 50:
# 数値は、位置を中心位置からずらすために、図の座標で指定
ax.text(fig_z[0], fig_z[1]-0.01, str(ival), size=14,
color="black", transform=ax.transAxes,
verticalalignment="top",
horizontalalignment="center")
サンプルコード
これまでのスクリプトを組み合わせると、次の天気図を作成できます。気象庁数値予報天気図の500hPa面天気図と同じ内容ですが、カラー化したことで見やすくなっていると思います。
この図を作成したスクリプト(jupyter Notebook用)を、次からダウンロードできます。作図たい方は、ご利用ください。
次回は、Qベクトルを算出して、その収束・発散の図を作成します。
この記事が気に入ったらサポートをしてみませんか?