バイナリダンプされた実数値の二次元配列をPythonから読む方法を紹介します。
具体的にはここにある地形データを読み込みます。
今回、対象とするのはMSMに使われている地形データです。グリッド数や緯度経度を見る限り、このディレクトリのデータが現状のMSMに対応していそうです。
READMEによると、
4. ファイル形式
(1) データ形式
4バイト実数(IEEE754)の2次元配列です。
バイトの並びはbig endianです。
(2) 座標系
等緯度経度格子座標系です。
(a) メソ数値予報モデルGPV
格子数: 東西方向 481 南北方向 505
格子間隔:東西方向 0.0625度 南北方向 0.05度
先頭の格子点:北緯47.6度 東経120度
(中略)
(3) データ格納順序
先頭の格子点から緯度の同じ格子点を経度方向東向きに格納し,
そのすぐ南の緯度で同様に繰り返し格納しています。
とのことです。
以下の例では、
地表ジオポテンシャル高度データ
(a) "TOPO.MSM_5K_20071121" メソ数値予報モデルGPV用
このデータを例に説明します。
まずはデータのダウンロード。
$ wget http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/etc/200709/TOPO.MSM_5K_20071121
あとはpythonで作業します。numpy
がインストールされているものとします。
import numpy as np
f = open("TOPO.MSM_5K_20071121",mode='rb')
topo = np.fromfile(f, dtype='>f',sep='').reshape(505,481)
これで、505x481のMSMと同じサイズで地形の高度データを読み込むことができます。
fromfileでdtype='>f'
としてbigendianの4バイト浮動小数としてデータを読み込んで、1次元配列に保存して、それをreshape
しています。
fromfileはdtypeを指定することで、1次元配列以外もとれるのですが、dt = np.dtype(('>f',(505,481)))
こういう感じで指定してもtopo.shape = (1,505,481)
のように3次元配列になってしまいます。
dtypeのドキュメントも確認したのですが、いい感じに一発で二次元配列にはできないみたいです。
正しく読み込まれているかを確認するために可視化してみましょう。
詳細はこちらの記事を参考にしてください。
%matplotlib inline
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import pygrib
grbs = pygrib.open('/temp/Z__C_RJTD_20170102000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin')
grb = grbs.select(forecastTime=0)[0]
lats, lons = grb.latlons()
flat_lats= np.ravel(lats)
flat_lons= np.ravel(lons)
m = Basemap(llcrnrlat=lats.min(),urcrnrlat=lats.max(), llcrnrlon=lons.min(),urcrnrlon=lons.max())
m.drawcoastlines()
m.contourf(flat_lons, flat_lats, b, latlon=True, tri=True)
plt.colorbar()
plt.show()
ちなみに海陸分布データだと以下みたいになります。