気象予報値（MSM-GPV）をpythonで音にしてみたい。

はじめに

R&Dアドベントカレンダーの6日目を担当します、NTTドコモクロステック開発部3年目社員の角谷昌恭です。
普段の業務では再生可能エネルギーやら電力やら、、通信会社では珍しい分野の研究・開発をしてますので、
本記事では普段の業務の中で目にすることが多い気象関係のデータの中から、特に気象予報値にフォーカスした記事を作成したいと思います。

そこで本記事では、MSM-GPVという気象庁作成の気象予報値の入手方法や使用ルール、データ構造に関する内容をメインに記載しますが、
ただ、それだとあまり面白い記事になりそうにないので、本データの活用方法の例を一つ記載したいと思います。

その活用方法の例というのがタイトルの通り、気象予報値を音にするってものなのですが、
これだけを聞くと何を言ってるんだ君は？って話だと思いますので、何故こんなことをしたいのかというお気持ちを簡単にイメージ図にしました。

イメージ通りのことができる様になった暁には、おもむろに電子音を聞いて、

「この音、、雨が来るな、、」

みたいな一発芸ができそう✌️
飲み会とかでヒーローになれそう✌️✌️

という妄想が膨らみ、本記事ではデータの活用例として、気象予報値（MSM-GPV）の音声化に取り組みます。

気象予報値（MSM-GPV）とは？

本章では気象予報値とは何ぞや？MSM-GPVとは何ぞや？といったことを少し真面目に記載していきます。

気象予報値とはその名の通り、実際に観測装置等を用いて測定したリアルタイムの気象情報ではなく、何らかの手法（衛生画像を用いた予測や、ゴリゴリの物理演算による予測など様々）により将来の気象情報を予測したデータのことを指します。
使用用途としては、例えば各地の時間別の天気予報を作成したり、太陽光パネルの発電量の予測を行う際の説明変数として用いられています。
今回はその気象予報値の中でも、過去データがアーカイブされ、公開されている「MSM-GPV」を取り上げたいと思います。

まずは「MSM-GPV」とは何ぞや？について簡単に説明します。
正式名称は「メソ数値予報モデルGPV（Grid Point Value）」で、気象庁がゴリゴリの物理演算で算出している気象予報値です。笑
基本的にリアルタイムの予報値を提供するサービスであり、本データは「一般財団法人　気象業務支援センター」が有料にて提供してますので、参考までに以下にURLを記載します。

「一般財団法人　気象業務支援センター」　->　http://www.jmbsc.or.jp/jp/online/file/f-online10200.htm

過去データについては後ほど説明致しますので、まずはリアルタイムに取得した場合のデータの概要について説明します。
データの種類には、39時間先まで予測したデータと、51時間先まで予測したデータの2種類あり（過去データに関しては15時間と33時間だったりします）、それぞれ以下のタイミングで計算が開始されおよそ2時間半で計算が終了します。

　39時間先　：　03、06、09、15、18、21UTC　（１日６回）
　51時間先　：　00、12UTC　（１日2回）

ちなみにUTCとは協定世界時のことで日本との時差は＋9時間あります。
なので03UTCで計算されたデータを日本で取得する時刻は以下の様になります。

　3(UTC) + 9(時差) + 2.5(計算時間) = 14.5(日本時間)

提供しているデータは上のURLにも記載されていますが以下の様な内容となっております。

　データ形式　：　国際気象通報式FM92 GRIB 二進形式格子点資料気象通報式(第2版)　＊通称GRIB2、後ほど説明致します。
　配信領域　　：　北緯 22.4度～47.6度、東経 120度～150度の区間を5kmメッシュ
　予報時間　　：　地上は1時間毎、気圧面は3時間毎
　要素　　　　：　海面更正気圧、地上気圧、風（東西成分と南北成分）、気温、相対湿度、時間降水量、雲量（全雲量、上層雲、中層雲、下層雲）、日射量

リアルタイム版は有料となっておりますが、過去データであれば「DIAS」というデータ統合・解析システムにより公開されております。
URLは以下の通りで、実際には東京大学生産技術研究所喜連川研究室が運営しているサイトにてダウンロードすることが可能です。
＊ただし本データの公開は「非商業的な研究及び教育用途のために提供している」ため、その利用については確認の後、自己責任でお願いします。
　また、DIASに掲載されているデータの使用についての注意書きもURL先のサイトにて確認できますので、必ずよく読んでからご使用ください。

「DIAS」　->　https://diasjp.net/
「GPV Data Archive」　->　http://apps.diasjp.net/gpv/

また、「GPV Data Archive」にてダウンロードする際に複数のファイルが存在し、どれを選択して良いか迷う可能性がありますので、その際は以下のファイル名の作成ルールをご参考ください。（個人の解釈によるものなので間違ってたらごめんなさい、、笑）

Z_C_RJTDYYYYMMDD030000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin

記号	意味
YYYY	西暦
MM	月
DD	日
03	UTC
Lsurf	地上：Lsurf、気圧面：L-pall
FH00-15	予測時間先情報、先の場合00-15、過去データは上述の予測時間の間隔とは異なる可能性があるので要注意

簡単なデータ構造の説明

前述した様にMSM-GPVは「国際気象通報式FM92 GRIB 二進形式格子点資料気象通報式(第2版)　＊通称GRIB2」という、WMO（世界気象機関）が定めるバイナリ通報式を採用しております。
何ぞや？と思う方が大半だと思いますが、公式のドキュメントを読む限りですと、GRIB2表を参照し、要素の情報をGRIB2データ中で記述する⾃⼰記述型のデータ形式のことの様です。
平たく言ってしまえば、バイナリデータとしてファイルフォーマット化し伝送する⽅式を採用しており、それを人語に翻訳するためにはGRIB２表を参照してくれってことだと思います。
以下に公式のドキュメントのURLを記載致しますので、興味のある方は確認してみてください。

「国際気象通報式」　->　https://www.jma.go.jp/jma/kishou/books/tsuhoshiki/kokusaibet/kokusaibet_22.pdf

本形式のバイナリデータを可視化するツールとして例えば「Wgrib2」「ecCodes」などもありますが、本記事ではpythonの「pygrib」というライブラリを使用して可視化する方法についてご紹介します。
pygribが何ものかについては本記事では割愛致しますので興味のある方は下記のURLで確認してみてください。

「Module pygrib」　->　https://jswhit.github.io/pygrib/docs/pygrib-module.html

では早速pygribの使用方法ですが、まずはconda等を使用して実行したい環境にライブラリをインストールして下さい。
＊可能であればdockerを使用すると良いです。というのもこのpygrib、インストールしimportするまでに環境によっては色々とエラーがでる可能性がありますので、融通が効きやすいdockerでの環境構築をお勧めします。以下のQiita記事にpygribを含んだdockerfileについて記載してある記事を見つけたので是非参考にしてみて下さい。

「pythonでgrib2フォーマットのファイルを触れる環境を用意する(Docker編)」 ->　https://qiita.com/mhangyo/items/8494a8039973ba220ce5#fnref1

インストールが完了致しましたら早速MSM-GPVのデータ構造をざっくり見ていきましょう。
まずは「open」を使用してバイナリデータを開封します。（例では2018年2月5日の地上、予測時間先00-15時間先を使用しております）

sample.py

import pygrib

gpv = pygrib.open\
('input/msmgpv_raw/Z__C_RJTD_20180205150000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin')
gpv_list = gpv.select(forecastTime=0, name='Surface pressure')
print(type(gpv_list))
gpv_message = gpv.select(forecastTime=0, name='Surface pressure')[0]
print(type(gpv_message))
print(gpv_message)

実行結果

sample.py

<class 'list'>
<class 'pygrib.gribmessage'>
2:Surface pressure:Pa (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 201802051500

ここで行っている作業を簡単に説明しますと、pygrib.openクラスを使ってファイルとのインターフェースを作成し、selectメソッドを使用し条件に合致するpygrib.gribmessageのリストを取得、取得したリストからpygrib.gribmessageを取得しています。
forecastTimeで指定できるのは予測時間先情報、例で使用しているデータでは0-15の範囲、nameで使用できるのはデータの要素になります。参考までにnameで指定できる変数名とその意味について簡単に記載します。

id	変数名	意味
0	Pressure reduced to MSL	海面更正気圧
1	Surface pressure	気圧
2	10 metre U wind component	風の東西成分
3	10 metre V wind component	風の南北成分
4	Temperature	気温
5	Relative humidity	湿度
6	Low cloud cover	低層の雲量
7	Medium cloud cover	中層の雲量
8	High cloud cover	上層の雲量
9	Total cloud cover	全雲量
10	Total precipitation	積算降水量
11	Downward short-wave radiation flux	日射量

＊selectで指定できる検索条件についてはいくつかあるみたいですが、自分はforecastTimeとname以外使用したことがありませんのでもし他の検索条件をご存知な方がいらっしゃったら教えて頂けると助かります。笑

現在の状態は選択した予測時間先の選択したデータ項目の全メッシュのデータが含まれている状態です。
ですので次は、必要な区間（緯度経度で指定）のデータだけを抽出したいと思います。
まずはlatolonsメソッドを使用すると、メッシュの緯度・経度を羅列したデータが確認できますので見てみましょう。

sample.py

lats, lons = gpv_message.latlons()
print(lats)
print(lons)

出力結果

sample.py

[[47.6  47.6  47.6  ... 47.6  47.6  47.6 ]
 [47.55 47.55 47.55 ... 47.55 47.55 47.55]
 [47.5  47.5  47.5  ... 47.5  47.5  47.5 ]
 ...
 [22.5  22.5  22.5  ... 22.5  22.5  22.5 ]
 [22.45 22.45 22.45 ... 22.45 22.45 22.45]
 [22.4  22.4  22.4  ... 22.4  22.4  22.4 ]]
[[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.    ]]

latsとlonsで格納形式が異なりますが、行列として表記した場合と地図上に表示した場合が視覚的に一致するようにこのような順番になっているのでは、、と思いますが真相はよくわかりませんので知っている方いらっしゃったら教えてください。笑
次にdataメソッドを使用して、指定した緯度・経度と最も近いメッシュ上のデータを取得してみます。
引数はlat1(緯度最小)、lat2(緯度最大)、lon1(経度最小)、lon2(経度最大)を指定します。
ただし、dataメソッドはメッシュ上に存在するlatとlonを指定する必要がありますので、少し工夫します。

sample.py

import numpy as np

def getNearestLatLon(lat, lon, message):
    """
    概要: 指定した緯度経度から最も近いメッシュ上の緯度経度を算出する
    @param lat: 対象の緯度
    @param lon: 対象の経度
    @param message: 対象のpygrib.gribmessage
    @return 対象値に最も近いメッシュ状のlat,lon
    """

    # メッシュ上のlatとlonのリストを作成
    lats, lons = message.latlons()
    lat_list = []
    for n in range(0, len(lats)):
        tmp = lats[n][0]
        lat_list.append(float(tmp))
    lon_list = lons[0]

    # リスト要素と対象値の差分を計算し最小値の値を取得
    lat_near = lat_list[np.abs(np.asarray(lat_list) - lat).argmin()]
    lon_near = lon_list[np.abs(np.asarray(lon_list) - lon).argmin()]

    return lat_near, lon_near

# 対象としたいlat、lon
lat = 35.0
lon = 139.0
lat_near, lon_near = getNearestLatLon(lat, lon, gpv_message)
print("lat_near : {} \nlon_near : {}".format(lat_near, lon_near))

# dataメソッドによる抽出
print(gpv_message.data(lat1=lat_near, lat2=lat_near, lon1=lon_near, lon2=lon_near))

出力結果

sample.py

lat_near : 35.00000000000072 
lon_near : 139.0
(array([[98104.17480469]]), array([[35.]]), array([[139.]]))

lat1とlat2、lon1とlon2に同じ値を指定することで指定したlatとlonに近いメッシュ上の点を取得することができました。

音にしてみよう！

MSM-GPVのダウンロードと大まかなデータ構造がわかったところで、早速音で気象予報値を表現してみたいと思います！
やりたいこととしては、MSM-GPVデータで取得可能な気象予報値に対して楽器を割り当て、その値の変動に応じて音階を変化させるというものです。
音の雰囲気から穏やかな1日なのか、慌ただしく天候が変化するかなどがわかるんじゃないかなと思ってます。（小並感）

sample.py

# 基本ライブラリ
import numpy as np
import pandas as pd
# pygrib
import pygrib
# pretty_midi
import pretty_midi

def getNearestLatLon(lat, lon, message):
    """
    概要: 指定した緯度経度から最も近いメッシュ上の緯度経度を算出する
    @param lat: 対象の緯度
    @param lon: 対象の経度
    @param message: 対象のpygrib.gribmessage
    @return 対象値に最も近いメッシュ上のlat,lon
    """

    # メッシュ上のlatとlonのリストを作成
    lats, lons = message.latlons()
    lat_list = []
    for n in range(0, len(lats)):
        tmp = lats[n][0]
        lat_list.append(float(tmp))
    lon_list = lons[0]

    # リスト要素と対象値の差分を計算し最小値の値を取得
    lat_near = lat_list[np.abs(np.asarray(lat_list) - lat).argmin()]
    lon_near = lon_list[np.abs(np.asarray(lon_list) - lon).argmin()]

    return lat_near, lon_near

def getGpvDf(lat, lon, gpv, max_time = 15):
    """
    概要: 指定した緯度経度から最も近いメッシュ上の緯度経度のGPVデータを抽出する
    @param lat: 対象の緯度
    @param lon: 対象の経度
    @param gpv: 対象日時のmsm-gpv
    @param max_time: 予測時間の最大値
    @return 対象値に最も近いメッシュ上のGPVデータをdataframeにて返す
    　　　　（index=変数、columns=予測時間の形式）
    """

    # 予測時間の範囲に応じて抽出するデータの範囲を変更
    if max_time == 15:
        for t in range(0, 15, 1):
            gpv_message = gpv.select(forecastTime=t)
            gpv_var_list = []
            for x in range(0, 11, 1):
                gpv_message_var = gpv_message[x]
                lat_near, lon_near = getNearestLatLon(lat, lon, gpv_message_var)
                gpv_data_var = gpv_message_var.data\
                    (lat1=lat_near, lat2=lat_near, lon1=lon_near, lon2=lon_near)
                gpv_var_list.append(gpv_data_var[0][0])
            if t == 0:
                gpv_df = pd.DataFrame(gpv_var_list)
            else:
                gpv_df = pd.concat([gpv_df, pd.DataFrame(gpv_var_list)], axis=1)
    else:
        for t in range(16, 33, 1):
            gpv_message = gpv.select(forecastTime=t)
            gpv_var_list = []
            for x in range(0, 11, 1):
                gpv_message_var = gpv_message[x]
                lat_near, lon_near = getNearestLatLon(lat, lon, gpv_message_var)
                gpv_data_var = gpv_message_var.data\
                    (lat1=lat_near, lat2=lat_near, lon1=lon_near, lon2=lon_near)
                gpv_var_list.append(gpv_data_var[0][0])
            if t == 16:
                gpv_df = pd.DataFrame(gpv_var_list)
            else:
                gpv_df = pd.concat([gpv_df, pd.DataFrame(gpv_var_list)], axis=1)

    return gpv_df

def getMelody(gpv_df):
    """
    概要: GPVデータのデータフレームから、各変数の変化を表現する音声ファイルを作成する
    @param gpv_df: 対象値に最も近いメッシュ上のGPVデータのdataframe
    @return GPVデータから作成したmidiファイル
    """

    #pretty_midiオブジェクトの作成
    pm = pretty_midi.PrettyMIDI(resolution=960, initial_tempo=120)

    # 各変数に楽器を割り当て、奏でる音階を値の変化から決定する
    for n in range(0, len(gpv_df.index), 1):
        # 楽器の割り当て、listは楽器の種別が重なりにくくするための工夫です
        instrument_list = [x*5 for x in range(0, 11, 1)]
        instrument = pretty_midi.Instrument(instrument_list[n])

        # 最小値0,最大値1にMin-Max Normalization
        gpv_df_norm = gpv_df.apply(lambda x : ((x - x.min())/(x.max()-x.min())),axis=1)
        for t in range(0, len(gpv_df), 1):
            tmp = gpv_df_norm.iloc[n,t]
            if 0<=tmp and tmp<=1/8:
                note_number = pretty_midi.note_name_to_number('C4')
            elif 1/8<tmp and tmp<=2/8:
                note_number = pretty_midi.note_name_to_number('D4')
            elif 2/8<tmp and tmp<=3/8:
                note_number = pretty_midi.note_name_to_number('E4')
            elif 3/8<tmp and tmp<=4/8:
                note_number = pretty_midi.note_name_to_number('F4')
            elif 4/8<tmp and tmp<=5/8:
                note_number = pretty_midi.note_name_to_number('G4')
            elif 5/8<tmp and tmp<=6/8:
                note_number = pretty_midi.note_name_to_number('A4')
            elif 6/8<tmp and tmp<=7/8:
                note_number = pretty_midi.note_name_to_number('B4')
            elif 7/8<tmp and tmp<=1:
                note_number = pretty_midi.note_name_to_number('C5')
            else:
                # 全ての値が0の場合適当な音を代入（正規化でNaNとなるため）
                note_number = pretty_midi.note_name_to_number('C5')

            note = pretty_midi.Note(velocity=100, pitch=note_number, start=t, end=t+1)
            instrument.notes.append(note)
        pm.instruments.append(instrument)
    pm.write('output/midi/test.mid') #midiファイルを書き込みます。

# 対象としたいlat、lon
lat = 35.0
lon = 139.0
# 対象としたい日時のMSM-GPVデータ
gpv = pygrib.open\
('input/msmgpv_raw/Z__C_RJTD_20180205150000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin')

gpv_df = getGpvDf(lat, lon, gpv, max_time = 15)
getMelody(gpv_df)

出力結果
以下のURLのoutputの中に、作成したmidiファイルとmp3ファイルを格納しましたので聞いてみてください。

「msmgpv_getmelody git」　->　https://github.com/sumiya-masayasu/msmgpv_getmelody

作成元となるMSM-GPVデータを変更することで音がしっかりと変わっていることはわかりますね。

おわりに

やりたいことの主軸はできたので筆者的には満足しています。
ただ、これだけでは到底一発芸への昇華は不可能ですので、今後取り組みたい改善点について以下に記載致します。

• 予測期間の中での変動を音で表現しているが、過去数年分データから抽出した最大値、最初値で正規化した方が良さそう
• 楽器をランダムに設定しているが、楽器にも相性があると思うのでちゃんと選んだ方が良さそう。
• 数値によって音階を設定しているが、音の強弱や音階の幅なども細かく表現できた方が表現に幅ができて良さそう
• 気象状況のイメージと音のイメージがリンクすると、もっと聴きやすくなりそう（雨の時に不協和音が流れるなど）

これら改善点を解決した暁には、実際に自身の体で本当に音を聞くだけで気象予報が理解できる人間になれるのかを検証してみたいと思います。笑
進展があれば以上の内容について追記できたらなと思いますが、一旦は本内容にてR&Dアドベントカレンダーの6日目を終了したいと思います。
最後まで読んでいただき、ありがとうございました！