grib2をpython(matplotlib)で地図上で可視化

前回はgrib2ファイルをpythonから操作する環境の構築について説明しました。
今回は読み込んだデータを。matplotlib+basemapで可視化する手法について紹介します。
特にJupyterについては説明しませんが、基本的にJupyter上で試しています(普通にpythonコマンドからでも大丈夫だとは思います)。

インストール

Anaconda環境であれば

apt-get -y install libgl-dev
conda install matplotlib
conda install basemap

でmatplotlibとBasemapがインストールできます。
apt-getのところは環境によって変わりますが、ないとmatplotlibのimport時にlibGL.soがないと言われてしまいます。
dockerの場合は、前回のDockerfileに上記の各行の頭ににRUNを付けて追記すればOKです。

Dockerfile

...以下を追記...
RUN apt-get -y install libgl-dev
RUN conda install matplotlib
RUB conda install basemap

grib2データ読み込み

前回説明したpygribを使って、gribデータを読み込みます。前回と同じサンプルを使います。

wget -P somepath/ http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2017/01/02/Z__C_RJTD_20170102000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin

ここからはpython内で作業していきます。
dockerで作業する場合には前回と同様、以下でcontainerを立ち上げて、その中でpythonを起動してください。

cd path_of_Dokcerfile
docker build -t pygrib_ubuntu .
docker run -it -v /somepath:/temp pygrib_ubuntu

import matplotlib
matplotlib.use('Agg') # サーバーでx環境がない場合などはこれがないとmatplotlibが動かない
import pygrib
import matplotlib.pyplot as plt
import numpy as np
import math
from mpl_toolkits.basemap import Basemap


grbs = pygrib.open('/temp/Z__C_RJTD_20170102000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin')

# selectメソッドで指定した要素(forecastTime=0)にマッチする要素の配列が返ってくる。
# forecastTime=0の要素は複数あるが先頭([0])は海面更生気圧(mslp)が入っている
grb = grbs.select(forecastTime=0)[0]

# lats,lonsは二次元配列で緯度経度が入っている
lats, lons = grb.latlons()

# 中に含まれているデータ(海面更生気圧)を二次元のndarrayで取り出す
mslp = grb.values

これで可視化に必要な、緯度経度とその地点のデータ(海面更生気圧)をgribから取り出すことができました。

可視化

続けて、天気図っぽく等圧線を描いていきます。

# 等圧線をひく間隔をlevelsにリストとして入れる
# MSLPの単位はPaなので、2hPaごとに線をひく場合、200おきに線をひくことになる。
levels = range(math.floor(mslp.min()/100)*100, math.ceil(mslp.max()/100)*100+1,200)

# lat,lonを一次元に変換
# 何故か二次元のままでは描画できない
flat_lats= np.ravel(lats)
flat_lons= np.ravel(lons)

fig = plt.figure()

# マップを作成。4つのパラメータは描画する範囲の指定
m = Basemap(llcrnrlat=lats.min(),urcrnrlat=lats.max(), llcrnrlon=lons.min(),urcrnrlon=lons.max())

# 等圧線をひく
m.contour(flat_lons, flat_lats,mslp,latlon=True,tri=True,levels=levels)

# 海岸線をひく
m.drawcoastlines() 

fig.savefig('/temp/mslp.png')

こんな感じで絵が出来上がります。

出力の際に以下のようなwarningsが出ますが、とりあえずは影響ないので気にしないことにしましょう。

/root/.pyenv/versions/anaconda3-4.1.1/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:3505: MatplotlibDeprecationWarning: The ishold function was deprecated in version 2.0.

  b = ax.ishold()

/root/.pyenv/versions/anaconda3-4.1.1/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py:3572: MatplotlibDeprecationWarning: axes.hold is deprecated.

    See the API Changes document (http://matplotlib.org/api/api_changes.html)

    for more details.

  ax.hold(b)