気象×Python 〜衛星降水データ〜

JAXAさんと共同研究していることもあって、衛星データを使う機会が多いんで、今回は数ある衛星データの中でもJAXAから提供されているGSMaPデータをpythonで読み出して、描画していきたいと思います〜
#1.　GSMaPとは

▶複数マイクロ波放射計搭載の衛星群観測によって、時間降水量を推定したプロダクトで、衛星全球降水マップ Global Satellite Mapping of Precipitation の略。
▶日本列島は気象庁AMeDASによる密な時間降水量データが存在しているものの、雨量観測の乏しい山岳域や海上域における降水量は得られないため、衛星により推定された降水量を使うメリットは大きいのでは。

てなわけで、普段はというとGrADSやGMTなどのツールを使って気象データを解析していますが、pythonの描画プログラムを見たところ、比較的容易にできそうだったので試みました。

#2.　データ入手方法
基本的にG-Portalとゆーサイトからユーザー登録を済ませておかないと取得できないと思います。
※登録しなくてもHDF5形式のサンプルデータならあります。今回はユーザー登録後にftpソフトウェア経由で手軽に取得できるbinaryデータを用います。
▶データセット詳細

物理量	空間分解能	時間分解能
降雨強度[mm/hr]	緯度経度0.1°格子	1時間

#3.　サンプルプログラム実行
早速、以下のサンプルプログラムを少々いじり、実行。
Basemapとか足りないモジュールがあれば予めインストール。

# basemapのインストールコマンド
$ conda install -c conda-forge basemap
$ conda install -c conda-forge basemap-data-hires

GSMaP.py

import matplotlib
matplotlib.rcParams['backend'] = 'TkAgg'
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
from mpl_toolkits.basemap import Basemap
import gzip
from struct import unpack

#各々自由にファイルの指定
filename = "gsmap_gauge.20141217.0000.v7.0001.0.dat.gz"
i = 0
rain = [0]*3600*1200
with gzip.open(filename, "rb") as f:
    while True:
        data = f.read(4)
        if data == "":
            break
        rain[i] = unpack("<f", data)[0]
        i = i+1

#読み込んだ降水量データを1200×3600の配列に変換
hprecipRateGC = np.reshape(rain, (1200, 3600))
#欠損値があれば0に置換
hprecipRateGC = np.where(hprecipRateGC < 0, 0, hprecipRateGC)

# binaryデータには緯度経度情報がないので、meshgrid関数で、データに一致する格子状の座標を生成。
lo = [5+10*i for i in range(0, 1800)]
lo.extend([-17995+10*i for i in range(0, 1800)])
Lo = 0.01*np.array(lo)
la = [5995-10*i for i in range(0, 600)]
la.extend([-5-10*i for i in range(0, 600)])
La = 0.01*np.array(la)
Lon, Lat = np.meshgrid(Lo, La)



fig = plt.figure()
#カラーの間隔を設定
interval = list(np.arange(1, 12, 1))
interval.insert(0, 0.1)
cmap = cm.jet
cmap.set_under("w", alpha=0)

m = Basemap(projection="cyl", resolution="i", llcrnrlat=20,
            urcrnrlat=50, llcrnrlon=120, urcrnrlon=150)
m.drawcoastlines(color='black')
m.drawmeridians(np.arange(0, 360, 30), labels=[True, False, False, True])
m.drawparallels(np.arange(-90, 90, 30), labels=[True, False, True, False])
x, y = m(Lon, Lat) 
im = plt.contourf(x, y, hprecipRateGC,
                  interval, cmap=cmap, latlon=True)

cb = m.colorbar(im, "right", size="2.5%")
plt.show()

実行結果がこちら、、、んん？？

なんか横に伸びてます
カラースケール変えたり、欠損値調べたり、配列変えたりしましたが中々改善しません。
致命的なエラーでもなさそうなんやけどなー

諦めて違うデータ形式で読み出そうとしたとき、プログラムガイドに以下の注意書きがありました。

とりあえず、経度方向の配列を反転させれば上手くいきそう。てか、Basemapちゃんと最新のバージョンなんだが、、

#4.　最終的なプログラム

new_GSMaP.py

# -*- coding: utf-8 -*-
import gzip
from struct import unpack
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap


class GSMaP:
    def __init__(self):
        filename = "gsmap_gauge.20141217.0000.v7.0001.0.dat.gz" 
        i = 0
        rain = [0]*3600*1200
        with gzip.open(filename, "rb") as f:
            while True:
                data = f.read(4)
                if len(data) < 4:
                    break
                rain[i] = unpack("<f", data)[0]
                i = i+1
        
        hprecipRateGC = np.reshape(rain, (1200, 3600)) #読み込んだ降水量データを1200×3600の配列に変換
        hprecipRateGC = np.where(hprecipRateGC < 0, 0, hprecipRateGC) #欠損値があれば0に置換
        a = hprecipRateGC[:, 0:1801][:, ::-1]
        b = hprecipRateGC[:, 1801:3601][:, ::-1]
        self.hprecipRateGC = np.concatenate([a, b],1)
    # meshgrid関数で、データに一致する格子状の座標を生成。
    def make_grid(self):
        lo = [5+10*i for i in range(0, 1800)[::-1]] 
        lo.extend([-17995+10*i for i in range(0, 1800)[::-1]])

        Lo = 0.01*np.array(lo)
        la = [5995-10*i for i in range(0, 600)]
        la.extend([-5-10*i for i in range(0, 600)])
        La = 0.01*np.array(la)
        self.Lon, self.Lat = np.meshgrid(Lo, La)

    def draw(self):
        fig = plt.figure(dpi=100)
        levels = [0.1, 0.5, 1, 3, 5, 8, 10]
        cmap = plt.get_cmap('jet')
        cmap.set_under("w")
        cmap.set_over("firebrick")

        m = Basemap(projection="cyl", resolution="i", llcrnrlat=20,urcrnrlat=50, llcrnrlon=120, urcrnrlon=150)
        m.drawcoastlines(color='black')
        m.drawmeridians(np.arange(0, 360, 5), labels=[True, False, False, True],linewidth=0.0)
        m.drawparallels(np.arange(-90, 90, 5), labels=[True, False, True, False],linewidth=0.0)
        im = plt.contourf(self.Lon, self.Lat, self.hprecipRateGC, levels, cmap=cmap, extend='both')
        cb = m.colorbar(im, "right", size="2.5%")
        cb.set_label('mm/hr')
        plt.title("GSMaP  9:00JST 17DEC2014")
        plt.savefig('gsmap_201412170900.jpeg')

if __name__ == '__main__':
    g = GSMaP()
    g.make_grid()
    g.draw()

変更前の配列	変更後の配列
無事描画できました。わっしょい。
#5.　参考
▶GSMaP関連の資料

・JAXAの公式HP
https://sharaku.eorc.jaxa.jp/GSMaP/index_j.htm
・GSMaPの全体的な説明書
https://www.eorc.jaxa.jp/GPM/doc/algorithm/GSMaPforGPM_20140902_J.pdf
・GSMaPのデータフォーマット解説書
https://sharaku.eorc.jaxa.jp/GSMaP/document/DataFormatDescription_MVK&RNL_v7.0000.0.pdf

▶GSMaP関連の論文

・GSMaPについての解説
https://ieeexplore.ieee.org/document/4261062
・GSMaP作成アルゴリズムの解説
https://www.jstage.jst.go.jp/article/jmsj/87A/0/87A_0_119/_article