はじめに

こんにちは。筆者は大学院で気象データの解析・研究に取り組み始め、はや２年が経とうとしているのですが、今回は先月に福岡や熊本で発生した線状降水帯発生時の大気の様子をPythonのmatplotlibを使って捉えてみたいと思います！

↓ 当時の雨雲レーダーの様子 (https://tenki.jp/ より引用)

長崎県から熊本県にかけて線状の雨雲エコーがくっきりと見えていますね

このとき日本周辺の大気はどんな様子だったのしょうか？
今回はPythonによる気象データの可視化手順について紹介したいと思います！

対象読者

Python matplotlib の基本的な使い方を理解している
matplotlib を使った気象データの可視化の流れを知りたい
Xarrayを用いた GPV データ(後述) の基本的な扱い方を知りたい
自分で大気の状態を可視化してみたい

完成イメージ

ぱっと見たい方向けに先に完成した図を載せます。

上空約1,500m の水蒸気量の様子
上空約1,500m の気温の様子

どうでしょうか？テレビのニュースで見たことがあるようなないような、、そんな図が出来上がりました。日本列島、特に九州に向かって中国大陸から多量の水蒸気が流入している様子が一目瞭然ですね。

また2枚目の図から、九州に向かって来る大気は周囲より幾分か暖かいことが分かります。

九州で線状降水帯が発生していたとみられる時間帯前後では、日本の西側から大雨の元となる暖湿流が流入するような大気場となっていたようです、、！

実装

ここからは具体的な可視化の手順および実装について見ていきます。

使用データ

今回用いたデータは、京都大学生存圏データベースのメソ数値予報モデルGPV です (以後 MSM とよびます) 。MSMは、気象庁が提供する日本周辺の大気環境場の解析値データセットであり、GPVとは格子点データのことで、3次元の格子点座標に気温や気圧などの各パラメータが格納されています。

ファイル形式はNetCDFとよばれるもので、気象学や海洋学の分野でよく用いられるバイナリファイルとなっています。詳しい説明は割愛しますが、データを配列として読み書きできることやデータの説明を合わせて格納する自己記述性をもっています。

イメージは「この階のこのエリアにあるこの車種の車を取り出して！」と頼むと、対応する車を取り出してくれる立体駐車場のようなものです。座標範囲やデータ変数名を指定すると、容易に対応する配列を取得することができます。

可視化の流れ

簡単に実装の流れを見ておきます。

Xarray を使ってNetCDFファイルの内容を xr.Dataset としてインスタンスに保持し、任意の配列を取得します。取得した配列を、cartopy で描画した白地図上に matplotlibを用いてplotしていきます。

実行環境

$ python -V
Python 3.11.9

xarray     : 2023.6.0
netCDF4    : 1.6.2
matplotlib : 3.8.4
cartopy    : 0.22.0

実装ステップ

NetCDFファイルの読み込み

京大生存圏サーバー上のNetCDFファイルにアクセスし、xr.Datasetとして読み込んだ結果を返す関数を用意しておきます。今回はローカルにファイルをダウンロードしたくないので tempfile に書き込んでおいて、それを xr.open_datasetメソッドで読み込みます。

fetcher.py

from tempfile import NamedTemporaryFile

import requests
import xarray as xr


def fetch_content(url: str) -> bytes:
    try:
        response = requests.get(url)
        response.raise_for_status()
        return response.content
    except requests.RequestException as e:
        print(f"Error fetching the file: {e}")
        raise


def read_nc_dataset(nc_content: bytes) -> xr.Dataset:
    with NamedTemporaryFile(suffix=".nc", delete=False) as tmp:
        tmp.write(nc_content)
        tmp.flush()
        tmpfile_path = tmp.name
    return xr.open_dataset(tmpfile_path)

取得したデータセットの確認

試しに、2025年8月10日の気圧面データにアクセスしてみます。

main.py

from fetcher import fetch_content, read_nc_dataset

url = "https://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/netcdf/MSM-P/2025/0810.nc"
contents = fetch_content(url)
ds = read_nc_dataset(contents)
print(ds)

$ python main.py 

<xarray.Dataset>
Dimensions:  (lon: 241, lat: 253, p: 16, time: 8)
Coordinates:
  * lon      (lon) float32 120.0 120.1 120.2 120.4 ... 149.6 149.8 149.9 150.0
  * lat      (lat) float32 47.6 47.5 47.4 47.3 47.2 ... 22.8 22.7 22.6 22.5 22.4
  * p        (p) float32 1e+03 975.0 950.0 925.0 ... 250.0 200.0 150.0 100.0
  * time     (time) datetime64[ns] 2025-08-10 ... 2025-08-10T21:00:00
Data variables:
    z        (time, p, lat, lon) float64 ...
    w        (time, p, lat, lon) float64 ...
    u        (time, p, lat, lon) float64 ...
    v        (time, p, lat, lon) float64 ...
    temp     (time, p, lat, lon) float64 ...
    rh       (time, p, lat, lon) float64 ...
Attributes:
    Conventions:  CF-1.0
    history:      created by create_1daync_msm_p.rb  2025-08-11

xr.Datasetとして上手く読み込めていますね。

Xarray について

Xarray は多次元配列を扱う Python のライブラリで、numpy 配列に次元(dimension)・座標(coordinate)・属性(attribute) などのメタデータを加えて扱うことができます。numpy や netcdf4 といったライブラリを用いて配列のみを扱うのに比べて、以下のような利点があります。

座標軸名でアクセスできる (次元の順番を気にしなくてよい)
sel メソッドにより、直感的に配列をスライスできる (indexアクセスが避けられる)
単位などの付加情報を保持できる

上記のような利点により、次元や配列のインデックスを気にする必要がなくなるため、配列の直感的な操作が可能になり、次元の取り間違えや、意図とは異なった範囲の配列取得といったミスを避けることができます。
その他、 Xarray の詳細な説明については、詳しい記事が豊富にあるのでそちらに譲りたいと思います。 ( 参考：xarray を用いたデータ解析 )

データセットの構造

先述したように、xr.Dataset は次元(dimension)と座標(coordinate)をもちます。少しややこしいのですが、dimensionは配列の軸そのものを表し、coordinateはそのdimensionの軸に対応する実際のラベル値(座標)になります。基本的には1つの次元(軸)に対し、1つの座標(ラベル値)がつくことが多いためどちらも同名となっていますが、time 次元に unix_time 座標と datetime座標を定義するというように、1つの次元に複数の座標を割り当てることも可能です。
これをもとに先ほどのデータセットの構造を把握しておきます。

$ python main.py

<xarray.Dataset>
# 次元は 経度緯度方向に241×253格子、高度方向に16層、時間軸方向に8時刻分
Dimensions:  (lon: 241, lat: 253, p: 16, time: 8)
# それぞれの次元と同じ名前で、座標が定義されている
Coordinates:
 * lon      (lon) float32 120.0 120.1 120.2 120.4 ... 149.6 149.8 149.9 150.0
 * lat      (lat) float32 47.6 47.5 47.4 47.3 47.2 ... 22.8 22.7 22.6 22.5 22.4
 * p        (p) float32 1e+03 975.0 950.0 925.0 ... 250.0 200.0 150.0 100.0
 * time     (time) datetime64[ns] 2025-08-10 ... 2025-08-10T21:00:00
# 格納されている変数名の情報
Data variables:
   z        (time, p, lat, lon) float64 ...
   w        (time, p, lat, lon) float64 ...
   u        (time, p, lat, lon) float64 ...
   v        (time, p, lat, lon) float64 ...
   temp     (time, p, lat, lon) float64 ...
   rh       (time, p, lat, lon) float64 ...
# 属性情報
Attributes:
   Conventions:  CF-1.0
   history:      created by create_1daync_msm_p.rb  2025-08-11

データ配列の抽出方法

先述したとおり、Xarrayでは次元名と座標値を指定してindexによらないアクセスが可能です。配列の抽出は非常に簡単で、ds.sel(次元名＝座標値)で任意の範囲の配列を取得することができます。

ds.sel(
    time="2025-08-10T18:00:00",
    p=850.0,
    lon=slice(135, 145),
    lat=slice(40, 30),    # lat座標は降順で格納されているので注意
)

> 
<xarray.Dataset>
Dimensions:  (lon: 81, lat: 101)
Coordinates:
  * lon      (lon) float32 135.0 135.1 135.2 135.4 ... 144.6 144.8 144.9 145.0
  * lat      (lat) float32 40.0 39.9 39.8 39.7 39.6 ... 30.4 30.3 30.2 30.1 30.0
    p        float32 850.0
    time     datetime64[ns] 2025-08-10T18:00:00
Data variables:
    z        (lat, lon) float64 ...
    w        (lat, lon) float64 ...
    u        (lat, lon) float64 ...
    v        (lat, lon) float64 ...
    temp     (lat, lon) float64 ...
    rh       (lat, lon) float64 ...

このように selメソッド及び slice を活用することで、任意の配列範囲を切り取ることができます。
では、一気に図の作成までいってみましょう。

データ配列の plot

main.py

from typing import cast

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from cartopy.mpl.geoaxes import GeoAxes
from fetcher import fetch_content, read_nc_dataset

# MSM-S にアクセスすると地表面データを取得可能
url = "https://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/netcdf/MSM-P/2025/0810.nc"
contents = fetch_content(url)
ds = read_nc_dataset(contents)

# 8/11 0300JSTの850hPa面値を指定 (格納時刻はUTC)
ds = ds.sel(time="2025-08-10T18:00:00", p=850.0)

# plotする変数を抽出
temperature = ds["temp"]
geo_height = ds["z"]
wind_u = ds["u"]
wind_v = ds["v"]

# 座標軸配列を取得
lons = ds["lon"]
lats = ds["lat"]

# figureを用意
fig = plt.figure(figsize=(10, 8))
proj = ccrs.PlateCarree()
ax = cast(GeoAxes, fig.add_subplot(1, 1, 1, projection=proj))

# 海岸線を描画
ax.coastlines(linewidth=1, resolution="10m")

# 日本付近を描画範囲に設定
ax.set_extent([120, 150, 22.4, 47.6], crs=proj)

# 陰影plot
shade = ax.contourf(
    lons,
    lats,
    temperature,
    levels=20,
    cmap="jet",
    extend="both",
    transform=proj,
)

# 等高線plot
contour = ax.contour(
    lons,
    lats,
    geo_height,
    colors="black",
    linewidths=0.5,
    transform=proj,
)

# ベクトルplot
vector = ax.quiver(
    lons,
    lats,
    wind_u,
    wind_v,
    transform=proj,
    regrid_shape=22,
    scale=24,
    color="grey",
    angles="xy",
    scale_units="xy",
)
plt.show()

↓ できた図

それっぽくなってきました。このままだと少し物足りないので、タイトル・緯度経度ラベル・カラーバー・等値線ラベル・凡例ベクトル等さらにプロットして完成を目指します。

完成版

main.py

from typing import cast

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
from cartopy.mpl.geoaxes import GeoAxes
from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter
from fetcher import fetch_content, read_nc_dataset
from matplotlib.axes import Axes
from mpl_toolkits.axes_grid1 import make_axes_locatable

# MSM-S にアクセスすると地表面データを取得可能
url = "https://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/netcdf/MSM-P/2025/0810.nc"
contents = fetch_content(url)
ds = read_nc_dataset(contents)

# 8/11 0300JSTの850hPa面値を指定 (格納時刻はUTC)
ds = ds.sel(time="2025-08-10T18:00:00", p=850.0)

# plotする変数を抽出
temperature = ds["temp"]
geo_height = ds["z"]
wind_u = ds["u"]
wind_v = ds["v"]

# 座標軸配列を取得
lons = ds["lon"]
lats = ds["lat"]

# figureを用意
fig = plt.figure(figsize=(10, 8))
proj = ccrs.PlateCarree()
ax = cast(GeoAxes, fig.add_subplot(1, 1, 1, projection=proj))

# 海岸線を描画
ax.coastlines(linewidth=1, resolution="10m")

# 陰影plot
shade_levels = np.arange(13.5, 24.5, 0.5)
temperature_celcius = temperature - 273.15
shade = ax.contourf(
    lons,
    lats,
    temperature_celcius,
    levels=shade_levels,
    cmap="jet",
    extend="both",
    transform=proj,
)

# カラーバー
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.2, axes_class=Axes)
fig.add_axes(cax)
cbar = fig.colorbar(shade, cax=cax, orientation="vertical")

# カラーバーのラベル
cbar.set_label("[℃]", labelpad=0.9, y=1.08, rotation=0, fontsize=10)

# 等高線plot
contour_levels = np.arange(1400, 1600, 10)
contour = ax.contour(
    lons,
    lats,
    geo_height,
    levels=contour_levels,
    colors="black",
    linewidths=0.3,
    transform=proj,
)

# 等高線ラベル
ax.clabel(contour, fmt="%d", fontsize=8)

# ベクトルplot
vector = ax.quiver(
    lons,
    lats,
    wind_u,
    wind_v,
    transform=proj,
    regrid_shape=22,
    scale=20,
    color="grey",
    angles="xy",
    scale_units="xy",
)

# ベクトルの凡例
legend_value = 20
ax.quiverkey(
    vector,
    X=0.92,
    Y=-0.08,
    U=legend_value,
    label=f"{legend_value} " + r"[$\mathrm{m\,s^{-1}}$]",
    labelpos="E",
    transform=proj,
)

# 緯度経度ラベル
lon_ticks = np.arange(120, 150 + 1, 5)
lat_ticks = np.arange(20, 50 + 1, 5)
ax.set_xticks(lon_ticks, crs=proj)
ax.set_yticks(lat_ticks, crs=proj)
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())

# 日本付近を描画範囲に設定
ax.set_extent([120, 150, 22.4, 47.6], crs=proj)

# タイトル
ax.set_title("2025/08/10 1800UTC 850hPa wind & temperature", fontsize=14)

plt.show()

plt.clf()
plt.close()

↓ 完成

<補足>

make_axes_locatable メソッドを用いることで、カラーバーをキャンパスのサイズに合わせてきれいに配置することができます。
ax.quiver メソッドの引数 regrid_shape や scale を調整することで、プロットするベクトルをいいかんじに間引いてくれたり、どれくらいの長さスケールで値を表現するかを決めることができます。regrid_shape ではどうやら内部で scipy が動いていて、単に間引くだけでなく周囲の値を重み付けしてプロットしてくれるようです。
1枚目の水蒸気混合比をplotするには近似式から算出する必要があります。

おわりに

今回は NetCDF 形式の GPV 気象データを Python matplotlibで可視化してみました。ちなみにもう少し上の高度を見てみると偏西風蛇行の様子も見えたりします。
↓ 上空約6000mの気温と風

たまには自分で手を動かして大気の状態を眺めてみるのもおもしろいかもしれません。
以上！！

参考記事

Xarray の基本： xarray を用いたデータ解析
matplotlibを用いた気象GPVデータの描画 :
PythonでGSMのGPVデータから気温の水平分布を描画する
NetCDFファイルの仕様 : PythonでnetCDFを扱う -気象データの読み取り-