GSMの500hPa面天気図作成 3

2021年10月3日 00:32

今回は、天気図上に低気圧などのマークを表示する方法について紹介します。コードの説明は、図上の座標取得、GPVから極値とその位置の取得、マークの描画の３点となります。

これで気象庁が作成している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ベクトルを算出して、その収束・発散の図を作成します。

この記事が気に入ったらサポートをしてみませんか？