pythonでGPVデータを元に天気図を作成する

はじめに

さて、本日12月1日は、なんの日かご存知でしょうか？

特に記念日というわけではないですが、気象に関連した話題で言うと、気象業務法が施行された日というのが1952年の12月1日らしいです。

そんな気象にまつわる日ということで、今日は気象データGPVを扱ったネタをお送りしようと思います。

GPVデータとは、Grid Point Valueの略で、気象庁が数値予報モデルで計算した結果を格子点ごとに表したデータです。気温、湿度、風速、気圧などの気象要素が含まれ、これを元に天気図が作成され、天気予報が行われます。

今回はpythonでこのGPVを使って実際に天気図を作図してみようと思います

作成を試みるのは、こちらの850hPa天気図（AUPQ78）とします。

（出典：気象庁）

この天気図には主に以下の要素が含まれます：

• 等高度線（実線）：850hPa面の高度分布
• 等温線（破線）：850hPa面の気温分布
• 湿域（ドット）：相対湿度が高い領域

それでは、これらの要素をGPVデータから取得して作図していきましょう。

GPVの取得

GPVデータは京都大学生存圏研究所が運営する生存圏データベースや、NOAA Operational Model Archives (NOMADAS)（最新〜10日前くらいのデータ）、National Centers for Environmental Informationのサイトからダウンロードすることができます。

今回は2025年11月30日0時UTCのGSMのGPVを使用します。

GPVデータの読み込み

ではここから実際にGPVデータを扱っていきます。

import pygrib
import numpy as np

# read gpv file 
fname_gfs = "xxxxxxxx.bin"

with pygrib.open(fname_gfs) as gribin:
  for g in gribin:
    print(g)

• GPVの読み込みは　pygribライブラリを使用します
• xxxxxxxx.binにはダウンロードしてきたGPVデータのパスを指定します
• 中身を一通り出力してみます
  以下のようにいろんな要素が格納されていることがわかります。

1:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 0 hrs:from 202510280000
2:Surface pressure:Pa (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202510280000
3:10 metre U wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 10 m:fcst time 0 hrs:from 202510280000
4:10 metre V wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 10 m:fcst time 0 hrs:from 202510280000
5:2 metre temperature:K (instant):regular_ll:heightAboveGround:level 2 m:fcst time 0 hrs:from 202510280000
6:2 metre relative humidity:% (instant):regular_ll:heightAboveGround:level 2 m:fcst time 0 hrs:from 202510280000
7:Low cloud cover:% (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202510280000
8:Medium cloud cover:% (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202510280000
9:High cloud cover:% (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202510280000
10:Total cloud cover:% (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202510280000
...

特定の気象要素、気圧面、予想時刻のデータを取り出す

fcst_hrs = 0
level = 850

with pygrib.open(fname_gfs) as gribin:
  dat = gribin.select(name='Temperature',level = level, forecastTime=fcst_hrs)
  for d in dat:
      print(d)
  
  
  #初期時刻と予報時刻
  ini_time = gribin.select(level = level, forecastTime=fcst_hrs)[0].analDate
  fcst_time = gribin.select(level = level, forecastTime=fcst_hrs)[0].validDate

selectで条件を与えてデータを絞り込みます。nameで気象要素、levelで気圧面、forecastTimeで初期時刻からの予想時刻を指定します

今回、forecastTimeは0なので初期時刻のデータを取り出しているので、初期時刻（ini_time）と予報時刻（fcst_time）は同じ時刻を示します

以下ではselectを使って、850面の気温、高度を取得します。湿数はGPVに要素としてありませんが、気温と相対湿度から求めることができるので、相対湿度を取り出しておきます。
この時必要に応じて、緯度経度（latとlon）の範囲を指定して、取り出す範囲も限定することも可能です。

lon_min = 70
lon_max = 180
lat_min= 0
lat_max =90

with pygrib.open(fname_gfs) as gribin:
  # 850 hpa temp
  temp,lat,lon = gribin.select(name='Temperature', level =level,forecastTime=fcst_hrs)[0].data(
    lat1=lat_min, lat2= lat_max, lon1=lon_min, lon2 = lon_max
  )
  temp = temp - 273.15  # K to C
  height,_,_ = gribin.select(name='Geopotential height', level =level,forecastTime=fcst_hrs)[0].data(
    lat1=lat_min, lat2= lat_max, lon1=lon_min, lon2 = lon_max
  )

  rh,_,_ = gribin.select(name='Relative humidity', level =level,forecastTime=fcst_hrs)[0].data(
    lat1=lat_min, lat2= lat_max, lon1=lon_min, lon2 = lon_max
  )

平滑化

取り出した気温データからそのまま等温線を書くこともできますが、そのままだとデータが細かくて等温線がギザギザになってしまうかもしれません。
そんな時は平滑化して等温線が滑らかになるようにします。

#平滑化
from scipy.ndimage import gaussian_filter
temp = gaussian_filter(temp, sigma=1.0)

sigmaの値を大きくすれば、より平滑化されますが、その分細かな情報が削ぎ落とされてしまいますので、適度な数で調整してください。

作図

作図の基本

さて、必要なデータを取り出したのでここからようやく作図していきます。
作図にはmatplotlibおよびcartopyライブラリを使用します。
以下では高度の等高線を書いています。

import matplotlib.pyplot as plt
import matplotlib.patheffects as path_effects
import cartopy.crs as ccrs
import cartopy.feature as cfeature

fig = plt.figure()
#plotエリアの設定
ax01 = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())

#等高度線の追加
cs01 = ax01.contour(lon, lat, height,
    levels= np.arange(0,3000,60),
    colors='black',
    linewidths=2,
    linestyles='-',)
cs02 = ax01.contour(lon, lat, height,
    levels= np.arange(0,3000,300),
    colors='black',
    linewidths=3,
    linestyles='-',)

contourでデータの等値線を書きます。
levelsで等高線の間隔を指定しています。60mごとのラインと合わせて300mごとのラインを書いて後者を少し太くしています。

等値線のラベルをつける

先ほどのコードの続きに以下を追加します。

fig = plt.figure()
〜〜
省略
〜〜

#ラベルの追加
clbls01 = ax01.clabel(cs01,fmt='%d')
for label in clbls01:
    label.set_rotation(0)
#ラベルバッファ
plt.setp(clbls01,path_effects=[path_effects.withStroke(linewidth=2,foreground='w')])

ただ単にラベルをつけるだけならclabelだけで大丈夫です。
set_rotation(0)でラベルを水平にしています。
path_effectsでラベルの周囲に白色のバッファをつけて、等高線が重なっていても数値がわかりやすいようになっています。

大陸や国境線の描画

さらに以下を追加します。

fig = plt.figure()
〜〜
省略
〜〜

#国境などの追加
ax01.add_feature(cfeature.COASTLINE, linewidths=2, color='gray')
ax01.add_feature(cfeature.BORDERS, linewidths=2, linestyle=':', color='gray')

ax01.add_feature(cfeature....)をつけることで海岸線や国境を追加することができます。cfeaturesは他にも、LAND：陸地、OCEAN：海洋、LAKES：湖沼、RIVERS：河川などがあります。

湿数の計算

湿数とは、簡単にいうと大気の湿り具合の指数です。GPVには直接湿数の値は入っていませんが、計算で求めることができます。
湿数＝気温- 露点温度
です。
露点温度は、気温と相対湿度から求めることができます

from metpy.units import units
import metpy.calc as mpcalc


# 湿数の算出
dp = mpcalc.dewpoint_from_relative_humidity(temp * units.degC , rh * units.percent)
wet = temp - dp.magnitude

求めた湿数をプロットします

fig = plt.figure()
#plotエリアの設定
〜〜
省略
〜〜

#湿数の追加
cs03 = ax01.contourf(lon, lat, wet,
    levels= np.arange(0,3,1),
    hatches = '......',
    colors='none'
    )

うまくドットで湿域が描画できたように思いますが、見本の天気図よりもドットの間隔が粗いような気がします。
これは画像のサイズが小さすぎるためなので、画像サイズを大きくします。

画像サイズを大きくする

fig = plt.figure(figsize=(20, 20))
〜〜
省略
〜〜

plt.figureにfigsize=(20, 20)などでサイズを指定します。

気温の追加

等高度線を追加したようにcontourで温度を追加します。

fig = plt.figure(figsize=(20, 20))
〜〜
省略
〜〜

#等温線の追加
cs04 = ax01.contour(lon, lat, temp,
    levels= np.arange(-30,30,3),
    colors='gray',
    linewidth=0.5,
    linestyles='--',
    ).clabel(fmt='%d',fontsize=15)

グリッドの追加

fig = plt.figure(figsize=(20, 20))
〜〜
省略
〜〜

#グリッドの追加
ax01.gridlines(
    draw_labels=True,
    xlocs = np.arange(80,190,10),
    ylocs = np.arange(0,70,10),
    linestyle='-',
    color='gray',
    linewidth=0.5,
)

• gridlines()で経緯度線を追加します。
• xlocs、ylocsでそれぞれ経度、緯度の描画範囲と間隔を指定します。

投影法の変更

元の天気図に比べると、現在の図面はどうも地図の見え方が違います。
これは元の天気図は、緯度60度、経度140度線を中心とした北極平射図法（NorthPolarStereo）という図法であるのに対して、現状の出力では正距円筒図法（PlateCarree）という図法で描画されているためです。
地図の投影法を変更します。

#マップの投影法の定義
crs_map = ccrs.NorthPolarStereo(
    central_longitude=140,
    true_scale_latitude=60.0
)

#データの座標系の定義
crs_data = ccrs.PlateCarree()

fig = plt.figure(figsize=(20, 20))

#plotエリアの設定
ax01 = fig.add_subplot(1,1,1, projection=crs_map)
#等高度線の追加
cs01 = ax01.contour(lon, lat, height,
    levels= np.arange(0,3000,60),
    colors='black',
    linewidths=2,
    transform=crs_data,
    linestyles='-',)
cs02 = ax01.contour(lon, lat, height,
    levels= np.arange(0,3000,300),
    colors='black',
    linewidths=3,
    transform=crs_data,
    linestyles='-',)

〜〜
省略
〜〜

#湿数の追加
cs03 = ax01.contourf(lon, lat, wet,
    levels= np.arange(0,3,1),
    hatches = '......',
    transform=crs_data,
    colors='none'
    )

#等温線の追加
cs04 = ax01.contour(lon, lat, temp,
    levels= np.arange(-30,30,3),
    colors='gray',
    linewidths=2,
    transform=crs_data,
    linestyles='--',
    ).clabel(fmt='%d',fontsize=15)

〜〜
省略
〜〜

• crs_mapではマップの描画の投影法を定義しています。今回は北極平射図法（NorthPolarStereo）としています。
  • central_longitude：基準となる経度を指定します。
  • true_scale_latitude：基準となる緯度を指定します。
  • plotエリアの設定の際にprojection=crs_mapでマップの投影法を指定しています。
• 一方で、crs_dataでは、描画するGPVのデータは緯度経度の座標系のデータであることを定義しています。
  • 等高度線や等温線などのデータを描画する際にtransform=crs_dataというパラメータでデータの座標系を指定する必要があります。

描画範囲を限定する

このままでは広範囲すぎるので描画範囲を絞り込みます。

fig = plt.figure(figsize=(20, 20))
〜〜
省略
〜〜

#描画範囲の限定
ax01.set_extent([90.0, 155.0,11.0, 50.0],crs_data)

set_extentで描画範囲を限定します。

テキストの追加

マップの余白に解析時刻やなんの気象要素が描画されているかなど、テキストを追加します。

fig = plt.figure(figsize=(20, 20))
〜〜
省略
〜〜

text01 =f"{level}hPa HEIGHT (M) , TEMP (℃) , WET AREA::(T-TD <3℃)"
fig.text(0.20, 0.16, text01, fontsize=20)

text02 =ini_time.strftime("%Y/%m/%d %H%M%SUTC + ") + str(fcst_hrs) + ":00"
fig.text(0.20, 0.14, text02, fontsize=20)

完成

元の天気図では高層気象観測の観測値が記載されていたりと、細かい点では相違がありますが、それらしい図ができました。

おまけ

contourだとラインで描画しますが、contourfにすると塗りになります。
ついでなので気温を塗りつぶしてみます。

無事天気図が完成しました。

最終的なコードは以下になります。

#投影法の定義
crs_map = ccrs.NorthPolarStereo(
    central_longitude=140,
    true_scale_latitude=60.0
)

#データの座標系を定義
crs_data = ccrs.PlateCarree()


fig = plt.figure(figsize=(20, 20))
#plotエリアの設定
ax01 = fig.add_subplot(1,1,1, projection=crs_map)
cs01 = ax01.contour(lon, lat, height,
    levels= np.arange(0,3000,60),
    colors='black',
    transform=crs_data,
    linewidths=2,
    linestyles='-',)
cs02 = ax01.contour(lon, lat, height,
    levels= np.arange(0,3000,300),
    colors='black',
    transform=crs_data,
    linewidths=3,
    linestyles='-',)


#add label
clbls01 = ax01.clabel(cs01,fmt='%d',fontsize =20)
for label in clbls01:
    label.set_rotation(0)
#ラベルバッファ
plt.setp(clbls01,path_effects=[path_effects.withStroke(linewidth=0.4,foreground='w')])

ax01.add_feature(cfeature.COASTLINE, linewidths=2, color='gray')
ax01.add_feature(cfeature.BORDERS, linewidths=2, linestyle=':', color='gray')



cs03 = ax01.contourf(lon, lat, wet,
    levels= np.arange(0,3,1),
    hatches = '......',
    transform=crs_data,
    colors='none'
    )
cs04 = ax01.contour(lon, lat, temp,
    levels= np.arange(-30,30,3),
    colors='gray',
    linewidths=2,
    transform=crs_data,
    linestyles='--',
    ).clabel(fmt='%d',fontsize=25)


cs05 = ax01.contourf(lon, lat, temp,
    levels= np.arange(-30,30,3),
    cmap='RdYlBu_r',
    alpha=0.7,
    transform=crs_data,
    )


#グリッドの追加
ax01.gridlines(
    draw_labels=True,
    xlocs = np.arange(80,190,10),
    ylocs = np.arange(0,70,10),
    linestyle='-',
    color='gray',
    linewidth=0.5,
)

#描画範囲の限定
ax01.set_extent([90.0, 155.0,11.0, 50.0],crs_data)


text01 =f"{level}hPa HEIGHT (M) , TEMP (℃) , WET AREA::(T-TD <3℃)"

fig.text(0.20, 0.16, text01, fontsize=20)

text02 =ini_time.strftime("%Y/%m/%d %H%M%SUTC + ") + str(fcst_hrs) + ":00"
fig.text(0.20, 0.14, text02, fontsize=20)
        

GPVデータの扱いについては、@LabCodeさんのこちらの動画も参考にさせていただきました。ありがとうございます。