気象業界では、数値予報や解析値などの面的なデータにGRIB2というバイナリフォーマットが使われます。@mhangyoさんのこちらの記事ではMatplotlib Basemap Toolkitを用いてGRIB2を可視化していますが、Basemapは開発終了となるため、ここではCartopyを用いた可視化に変えてみます。
（2019-02-28に書いた記事です。）
はじめに：Matplotlib Basemap ToolkitからCartopyへ
pandas等で読み込んだデータを地図にプロットしたいとき、Matplotlib Basemap Toolkitが広く使われていると思うのですが、最近はCartopyへの移行が強く推奨されています。以下に、公式サイトのアナウンスの訳（一部意訳）を載せておきます。
Cartopyとプロジェクト管理の新たなフェーズ、そしてサポート終了のお知らせ
2016年以降、Basemapのソフトウェアプロジェクト管理は新たなフェーズに入りました。CartopyというプロジェクトがBasemapを置き換えていく予定です。CartopyはまだBasemapのすべての機能を実装しているわけではありません。しかし、新しいソフトウェア開発においては可能な限りCartopyを使ったほうがよいですし、Basemapを使っている既存のソフトウェアはCartopyへの移行プロセスを開始したほうがよいでしょう。すべての保守や開発の取り組みはCartopyに集中させるべきです。
Ben Rootさんが2020年までのBasemapの保守を進んで引き受けてくれましたので、プルリクエストの査読やバグ修正は行われる予定です。また、NumPyやMatplotlibなどのパッケージとの互換性を保証するための保守もなされます。新機能の追加も受け入れます。しかし、強調しておきますが、新たな開発の取り組みはCartopyに集中させるべきです。というのも、Python 2.7の公式サポートが2020年に終了するときに、Ben RootさんはBasemapのリリースをしてサポートを終えるつもりだからです。
環境準備
まず、pygribとCartopyの入った環境を用意しましょう。
自分の場合、最近、環境作りではDockerを用いることが多く、特にJupyterを用いたデータサイエンス環境を作りたい場合はたいていscipy-notebookを拡張しています。今回は、以下のようなDockerfileで行ないました。ハマったところもあるので、おまけとして記事の最後に書いてあります。
FROM jupyter/scipy-notebook
USER root
RUN conda install --quiet --yes \
'pygrib' \
'cartopy' && \
conda clean -tipsy && \
fix-permissions $CONDA_DIR
RUN pip install --upgrade numpy
データは、@mhangyoさんの記事と同じMSM（Meso-Scale Model, 気象庁メソモデル）のデータファイルを用意しておきます。
GRIB2の読み込みから可視化まで
Pythonでの読み込みと可視化に入っていきます。
以下のコードは、先程のDockerfileからビルドしたDockerイメージのJupyterで動作を確認していますが、matplotlibの出力先などを適切に設定してあげれば他の環境でも動くはずです。
import
まずはimportです。pygribとCartopy以外はよく見るimport文です。
%matplotlib inline
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import pygrib
import seaborn as sns
読み込み
GRIB2を読み込んで、緯度、経度、格子点値を変数に格納します。読み込みなので、当然ながらBasemap使用時からの変化は特にありません。
grbs = pygrib.open('Z__C_RJTD_20170102000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin')
grb = grbs.select(forecastTime=0)[0]
lats, lons = grb.latlons()
mslp = grb.values
可視化の準備
続いて、可視化のためにデータを調整していきます。
今回は、等圧線にラベルも入れることにします。気圧の値は、GRIB2内にはPa単位で入っていますが、気象業界ではhPa単位で表すのが通例です。今回もラベルはhPa単位で表示したいので、Pa単位の格子点値をすべてhPa単位に変換しておきましょう。
mslp_hpa = mslp / 100
また、等高線はhPaが2の倍数となるところに引きたいので、2の倍数のみをレベルとするように調整します。
levels = np.arange(mslp_hpa.min()//2*2, mslp_hpa.max()//2*2+2, 2)
正距円筒図法でのプロット
まずは、Basemapのデフォルト（projection='cyl'）と同じ正距円筒図法でプロットしてみましょう。正距円筒図法は、緯度と経度が直角かつ等間隔に交差する最もシンプルな図法で、緯度と経度の格子データとして表される気象データとは相性がよい図法でもあります。
lat_s = lats.min()
lat_n = lats.max()
lon_w = lons.min()
lon_e = lons.max()
proj = ccrs.PlateCarree() # 正距円筒図法を指定
ax = plt.axes(projection=proj)
ax.gridlines()
ax.coastlines(resolution='10m')
ax.set_extent((lon_w, lon_e, lat_s, lat_n), proj) # 緯度経度の範囲を指定
CS = ax.contour(lons, lats, mslp_hpa, levels, transform=proj)
ax.clabel(CS, fmt='%.f') # 気圧を整数で表示するよう、fmtを指定
正距円筒図法は気象データとの相性はよいのですが、高緯度の地域で潰れた感じになるのが少し気になりますね。
正射図法でのプロット
続いて、正射図法で表示してみます。正射図法は、地球の表面をそのまま真横の平面へと投影した図法です。遥か彼方から地球を眺めたときの半球面を描くことになるため、自然な感じになります。
正距円筒図法は、すべての緯線と経線を平等に扱うため、特にパラメータの指定は必要ありませんでした。しかし正射図法では、「どこを中心に投影するか」、すなわち「どこを中心に眺めるか」で見たい地域の形が変わってくるため、投影の中心の緯度経度を指定してあげる必要があります。例えば、日本の中心を投影の中心とすれば日本はきれいに描けますが、太平洋の中心を投影の中心とすると日本はかなり歪んでしまいます。
ここでは、描画範囲の中心を投影の中心にしてプロットしてみましょう。
lat_s = lats.min()
lat_n = lats.max()
lon_w = lons.min()
lon_e = lons.max()
proj = ccrs.Orthographic((lon_w + lon_e) / 2, (lat_s + lat_n) / 2) # 対象緯度範囲、経度範囲の中心を投影の中心とした正射図法を指定
latlon_proj = ccrs.PlateCarree() # 緯度経度の処理用に正距円筒図法も使う
ax = plt.axes(projection=proj) # 地図は正射図法で表示
ax.gridlines()
ax.coastlines(resolution='10m')
ax.set_extent((lon_w, lon_e, lat_s, lat_n), latlon_proj) # 緯度経度で指定した描画範囲は正距円筒図法で処理
CS = ax.contour(lons, lats, mslp_hpa, levels, transform=latlon_proj) # 緯度経度の座標値は正距円筒図法で処理
ax.clabel(CS, fmt='%.f')
一つの描画用のソースコードで図法が2つ登場するので、最初は混乱しやすいかもしれません。
**「緯度経度を指定する箇所では、正距円筒図法（cartopy.crs.PlateCarree）を必ず使う」**と考えるとよいかと思います。
おまけ：Dockerでのpygrib環境作りでハマったところの共有
今回、pygribを入れるため、最初は以下のようなDockerfileを用意しました。
FROM jupyter/scipy-notebook
USER root
RUN conda install --quiet --yes \
'pygrib' \
'cartopy' && \
conda clean -tipsy && \
fix-permissions $CONDA_DIR
ところが、イメージのビルドは成功したように見えるものの、立ち上げると、以下のようにpygribのimportができません。
import pygrib
print(pygrib.__version__)
出力：
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
RuntimeError: module compiled against API version 0xc but this version of numpy is 0xb
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
<ipython-input-1-7b90afe9041e> in <module>
----> 1 import pygrib
2
3 print(pygrib.__version__)
ImportError: numpy.core.multiarray failed to import
Anacondaでこのようなことが発生した場合、NumPyを入れ直してあげるとよいようです。本記事冒頭のDockerfileを再掲します。
FROM jupyter/scipy-notebook
USER root
RUN conda install --quiet --yes \
'pygrib' \
'cartopy' && \
conda clean -tipsy && \
fix-permissions $CONDA_DIR
RUN pip install --upgrade numpy
このイメージだと、きちんとimportが通ります。
import pygrib
print(pygrib.__version__)
出力：
2.0.3
一つの記事にするほど深く調べたわけでもないので、ついでとしてこちらに書いておきました。