pythonで格子点データを扱う（入門編）

はじめに

前回は、wgrib2というコマンドを使って、GRIB2形式データをcsv形式に変換したり、Fortranで読み込む方法を解説しました。
今回はpythonで格子点データを扱う方法についてのお話です。
というのも、pythonは、データサイエンス分野では欠かせない言語であり、豊富なモジュールで統計手法や作図を行うのに広く使われていて、気象データと何か（例えば、商品の販売数など）との関係を分析する（関係式を作る）のに強力なツールだと思います。
pythonでデータの読み方を理解できれば、あとは様々な応用ができると思います。csv形式等テキストの読み方は、他のサイトでも情報が得られると思いますので、特に気象データとして特徴的なバイナリ形式の格子点データを扱う方法について解説したいと思います。

環境として必要なのは、今回もwgrib2。また、python3はインストール済みであることとします。
pythonのモジュールとして、numpyを使います。matplotlib、Basemapも出てきますので、pipなどでインストールしておきましょう。

• 前提条件
  • wgrib2インストール済み（方法は前回記事を参照）
  • python3インストール済み
  • numpyインストール済み
  • matplotlib、Basemapインストール済み（作図する場合）

wgrib2で単精度浮動小数点バイナリに変換したデータを格子点配列で読み込む

せっかくpythonを使うのだから、GRIB2形式を直接pythonで読み込む方法をお話した方が良いのかもしれませんが、wgrib2は要素毎や時刻毎に処理したり、GRIB2形式ファイルを扱うのがとても簡単なので、wgrib2で2次元データを一つ抽出して単精度浮動小数点バイナリに変換したデータを使うことにします。大量のデータを処理する場合はシェルスクリプトを組むと良いと思います。（もっと簡単に格子配列を読み込む方法をご存知でしたら教えてください）
今回もGRIB2形式データはMSM-GPV(地上:Z__C_RJTD_yyyymmddhh0000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin)を用います。
単精度バイナリファイルを出力するので、コマンドは、

$ wgrib2 Z__C_RJTD_20200130000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin | grep "TMP" | grep "anl" | wgrib2 -i -no_header Z__C_RJTD_20200130000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin -no_append -ieee msm.bin

こんな感じです。地上気温の初期値を出力しています。この出力ファイルmsm.binは、4バイト毎に東西481格子、南北505格子の順に南西端から格納されています。

$ ls -la msm.bin
-rw-rw-r-- 1 user user 971620 1月 30 17:31 msm.bin

バイナリのエンコードはビッグエンディアンです。
これをpythonで読み込んでみます。

import numpy as np

# 設定
imax = 481 # 東西格子数
jmax = 505 # 南北格子数
datasize = imax*jmax*4 # データサイズ

# バイナリファイルの読み込み
binfile = open('msm.bin','rb')

# 2次元配列に変換
array = np.fromfile(binfile, dtype='>f', sep='',count=datasize).reshape(jmax,imax)

print(array)

openで'rb'を指定してバイナリを読み込みます。
numpyのfromfileでバイナリデータから配列に変換します。reshapeは東西南北格子の2次元に整形しています。
実行すると以下のように出力されます。

[[289.67145 289.53082 289.37457 ... 299.09332 299.0777 299.06207]
[289.54645 289.37457 289.18707 ... 299.09332 299.0777 299.04645]
[289.40582 289.21832 288.98395 ... 299.09332 299.0777 299.04645]
...
[243.9683 241.67143 241.82768 ... 271.43707 271.43707 271.43707]
[246.14018 242.39018 241.79643 ... 271.3902 271.3902 271.3902 ]
[249.37456 245.68706 244.20268 ... 271.35895 271.35895 271.34332]]


このままではちょっとわかりにくいのですが、2次元配列で読み込むことができました。

csv形式に出力する

バイナリファイルを読み込めてしまえば、あとは自由に処理することができますね。
とりあえず、視覚的にデータを確認するために、csv形式に出力してみましょう。

バイナリファイルには緯度経度の情報がないので、計算して2次元配列を作ります。

ilon = 120.0 # 最初の経度
ilat = 22.4 # 最初の緯度
dlon = 0.0625 # 経度方向の格子間隔
dlat = 0.05 # 緯度方向の格子間隔

lons = (ilon + dlon * np.indices((jmax,imax))[1])
lats = (ilat + dlat * np.indices((jmax,imax))[0])

print(lons)
print(lats)


indicesは配列番号を返します。array変数と同じように2次元配列にします。
出力：

[[120. 120.0625 120.125 ... 149.875 149.9375 150. ]
[120. 120.0625 120.125 ... 149.875 149.9375 150. ]
[120. 120.0625 120.125 ... 149.875 149.9375 150. ]
...
[120. 120.0625 120.125 ... 149.875 149.9375 150. ]
[120. 120.0625 120.125 ... 149.875 149.9375 150. ]
[120. 120.0625 120.125 ... 149.875 149.9375 150. ]]
[[22.4 22.4 22.4 ... 22.4 22.4 22.4 ]
[22.45 22.45 22.45 ... 22.45 22.45 22.45]
[22.5 22.5 22.5 ... 22.5 22.5 22.5 ]
...
[47.5 47.5 47.5 ... 47.5 47.5 47.5 ]
[47.55 47.55 47.55 ... 47.55 47.55 47.55]
[47.6 47.6 47.6 ... 47.6 47.6 47.6 ]]


ここまでの変数を用いて、csv形式ファイル（経度,緯度,気温）に出力します。

dataline = ""
for j in range(jmax):
   for i in range(imax):
      dataline = dataline + str(lons[j,i]) + ',' + str(lats[j,i]) + ',' + str(array[j,i]) + "\n"

outfile = open('msm.csv', 'w')
outfile.write(dataline)


ただ、この方法だとwriteにすごく時間がかかります。printで出力して、リダイレクトでファイル保存した方がはるかに高速です。例えば、このプログラムの名前がbin2csv.pyとすると、プログラムは、

for j in range(jmax):
   for i in range(imax):
      print(str(lons[j,i]) + ',' + str(lats[j,i]) + ',' + str(array[j,i]))

と書いておいて、

$ python3 bin2csv.py > msm.csv

このように打ちます。writeでもprintでも結果は同じです。私の環境だと、writeでは3分程度、printでリダイレクトでは3秒程度でした。

おまけ：作図して確認する

読み込みが正しいかどうかは分布図を作ってみると良いと思います。
pythonのmatplotlibは様々な可視化ができます。また、Basemapを使って、地図上にデータをプロットしてみました。

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

# プロット領域
fig = plt.figure()
# マップを設定（中程度の解像度の海岸線、等緯度経度投影、東経135度を中心、東経115.5～155.5度、北緯20～50度）
m = Basemap(resolution='i', projection='cyl', lon_0=135, llcrnrlon=115.5, urcrnrlon=155.5, llcrnrlat=20, urcrnrlat=50)
# 緯度経度から投影座標に変換
x,y = m(lons,lats)
# 経線（東経110～160度、5度間隔、グレー色、図の南にラベル表示）
m.drawmeridians(np.arange(110,160,5), color="gray", fontsize='small', labels=[False, False, False, True])
# 緯線（緯度10～50度、5度間隔、グレー色、図の西にラベル表示）
m.drawparallels(np.arange(10,60,5), color="gray", fontsize='small', labels=[True, False, False, False])
# 海岸線
m.drawcoastlines()
# データプロット
levels = np.arange(np.floor(array.min()),np.ceil(array.max()),4) # 最小値～最大値まで4K毎にラベル表示
cs = m.contour(x,y,array,levels=levels) # 等値線描画
cs.clabel() # 等値線にラベルを表示
# 出力
plt.title('MSM-GPV t2 anl')
plt.show()


こんな感じにプロットされます。東西南北の格納順も合っていますので、適切に読み込めていることがわかります。

まとめ

pythonを使って格子点データを読み込み、csv形式出力や作図の方法を解説しました。
データによってはUNDEF値や値変換が必要な場合があるので、2次元配列に読み込んでから処理を入れた方が良いかもしれません。
これで格子点データを使って、機械学習やディープラーニングもできるでしょうか。応用編としてまたお話できたらと思います。