気象データをもとに「天気図っぽい前線」を機械学習で描いてみる(2)

### 気象データをもとに「天気図っぽい前線」を機械学習で描いてみる（2) 入力となる気象データ可視化の話

前回、[気象データをもとに「天気図っぽい前線」を機械学習で描いてみる（1）](https://qiita.com/m-taque/items/80a87052fde00fb96fab "気象データをもとに「天気図っぽい前線」を機械学習で描いてみる (1)")として、全体の流れを投稿しました。
第2回目は、GPV形式の数値データを可視化するまでの話をまとめます。

・2020.6.26 5.4の鉛直流可視化例の注釈の誤記修正

## 0.私の環境

私が利用しているのはMacOSの環境です。
それまでMac book proやiPhone(Pythonista)を使って勉強していましたが、Mac miniがリニューアルされたのを機に購入しました。メモリを増強したのが結果的に良かった。

###ハードウェアとソフトウェア
Mac mini(2018)
　プロセッサ　3.2GHz 6コア Intel Core i7
　メモリ　　　32GB 2667 MHz DDR4

OS macOS Catalina
導入当初は1世代古いmojaveでした。

###開発環境

統合開発環境的なものはatomも入ってはおりますが、あまり使っておりません。
viで何でもやって来た世代ですので・・・


## 1.前線に関わる気象データ
### 1.1 前線を描きこむ対象は？

**前線**は気象庁で気象状況の解析を行い、**アジア地上解析図(ASAS)**や**速報天気図(SPAS)**の中に描き込まれています。

[気象庁ホームページ 天気図について](http://www.jma.go.jp/jma/kishou/know/kurashi/tenkizu.html "気象庁ホームページ 天気図について")

サンプルとして2019年12月2日21時UTCのものを示します。
![SPAS_COLOR_201912021200.v8.512.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/221deeb2-c439-272a-4e00-b4818af94367.png)

これが目指す**天気図**です。

一方で、今回京都大学から入手する気象データは気象庁の数値予報モデルのうちの一つである、**全球予報モデル（GSM）**です。

数値予報などの解説は下記をご覧ください。
[気象庁ホームページ 数値予報とは](https://www.jma.go.jp/jma/kishou/know/whitep/1-3-1.html "数値予報とは")

GSMのデータには、**気温**や**気圧**といった気象要素のデータは含まれていますが、**どこが前線であるか**といったようなことはわかりません。
今回の機械学習は、このGSMのデータに対して、**天気図**のように「前線を描くべき場所がどこかを認識して前線を描く」ということを学習させるものと言えます。

そして、GSMの初期時刻のデータを対象とすることにします。気象庁が解析を行い天気図を作成している対象時刻と同じであり、ほぼ同じ気象状況が表現されていると考えられるからです。

<details><summary>初期時刻データと観測値</summary>
GSMは、初期時刻のデータをもとにして、例えば３時間後、６時間後、・・・というように将来の大気の状態を計算していきます。ここでいう「初期時刻のデータ」とは、ある時間の有限の観測点から規則正しい格子の位置での、”最もたしからしい大気の状態”を統計的な手法を駆使して求めたものです。これを「客観解析」といい、最近様々な分野へ応用され始めている「データ同化」が以前から用いられています。
[気象庁ホームページ 客観解析](https://www.jma.go.jp/jma/kishou/know/whitep/1-3-3.html "客観解析")

</details>

### 1.2 前線を決める気象要素

前線は性質の違う気団同士の境目であって通常は**暖気**と**寒気**の境目になりますが、気象データを可視化するとそこら中にそういった境目は存在しているように見えます。

まずはどうやって前線を描く場所を決めているのか？
気象学の専門書にはあまり具体的な話は載っておりませんが、気象予報士試験の参考書にそれらしい記載があることがわかりました（気象予報士試験では、天気図を見て気象前線を描かせる問題が出るため、前線解析はひとつの大きなジャンルなのです）。

[気象予報士 かんたん合格テキスト](https://www.amazon.co.jp/気象予報士かんたん合格テキスト-%E3%80%88実技編〉-らくらく突破-気象予報士試験受験支援会/dp/4774149950/ref=pd_sbs_14_img_1/355-1239468-1368320?_encoding=UTF8&pd_rd_i=4774149950&pd_rd_r=5f510193-a636-44b6-9ee7-011e7278ba9c&pd_rd_w=wCuCw&pd_rd_wg=ZGywm&pf_rd_p=ca22fd73-0f1e-4b39-9917-c84a20b3f3a8&pf_rd_r=R77VP37VM7V6TX7BMPWE&psc=1&refRID=R77VP37VM7V6TX7BMPWE "気象予報士 かんたん合格テキスト")
![51RlzcuwIBL._SX351_BO1,204,203,200_.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/3d1fa85b-d92f-1ad3-3ae2-b53755d9f330.jpeg)


それによると、

1. **等温線**が集中している部分の暖気側の境界(地上、850hPa面)
2. **等相当温位線**が集中している部分の南端(850hPa面）
3. **風速や風向**が急に変化している部分
4. **湿潤域**や**上昇流域**なども総合的に見て判断

ということをもとに前線を描くようです。
ということで、これらに相当する気象要素を使用するのがよさそうです。

## 2. pygrib

### 2.1 pygribとは？

pygribはGRIB2形式のファイルを読み取ることができるライブラリです。
インストールの方法は例えば以下のような先達の方の解説があります。

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

### 2.2 GRIBファイルをオープンする

`open`メソッドによりファイルを開いてオブジェクトを得ます。

```python:pygrib-open.py
import pygrib
import sys
import numpy as np
import math

(中略）

# GRIB file settings.
_sourcedir = "./hogehoge/"
_GSMfilename = _sourcedir + "Z__C_RJTD_" + _yearstr + _monstr + _daystr + _hrstr + "0000_GSM_GPV_Rgl_FD0000_grib2.bin"

# Open GRIB file
grbs = pygrib.open(_GSMfilename)

```

`_yearstr`、`_monstr`、 `_daystr`、 `_hrstr`には対象とする年、月、日、時を文字列として与えます。
オブジェクト`grbs`から、いくつかのパラメータを指定して`select`メソッドを使用することで目的のデータを抽出します。

### 2.3 気象要素、高度面、予報時刻を指定して抽出する

今回使うのは、`parameterName`、`level`、`forecastTime`の３つです。

データを公開している機関や気象モデルによって具体的なパラメータ名称は異なります。
日本の気象庁が公開しているデータの場合、

[気象庁ホームページ：配信資料に関する技術情報](https://www.data.jma.go.jp/add/suishin/cgi-bin/jyouhou/jyouhou.cgi "配信資料に関する技術情報")

に解説がありますが、pygribの使用は想定されていないのか、parameterNameに何を指定するべきか？については窺い知ることはできません。

表 2-1. parameterNameの指定パラメタ

| 指定パラメタ | 内容                   |
|:----------|:-----------------------|
|Pressure reduced to MSL | 海面更正気圧(Pa単位になっています）  |
| Temperature  | 気温（絶対温度になっているので注意）  |
| Relative humidity  | 相対湿度  |
| Vertical velocity [pressure]  | 鉛直風速度（気圧座標）  |
| u-component of wind  | 水平風速の東西成分  |
| v-component of wind  | 水平風速の南北成分  |
| Geopotential height  | ジオポテンシャル高度  |
| Low cloud cover      | 下層雲雲量         |
| Medium cloud cover   | 中層雲雲量        |
| High cloud cover     | 上層雲雲量　　　　|
| Total cloud cover    | 全雲量　|

表 2-2. level の指定パラメタ

| 指定パラメタ | 内容                   |
|:---------:|:-----------------------|
|850| 850hPa面のデータ |
|700| 700hPa面のデータ |
|500| 500hPa面のデータ |
|300| 300hPa面のデータ |
|0  | 地表面データ（海面更正気圧で使用）|
|2  | 地表面データ（気温および相対湿度で使用）|
|10 | 地表面データ（水平風速で使用）|


表 2-3. forecastTime の指定パラメタ

| 指定パラメタ | 内容                   |
|:---------:|:-----------------------|
|0|予報初期時刻|
|6|初期時刻から６時間後のデータ|
|12|初期時刻から12時間後のデータ|
|| （以降6時間間隔で続く）|





これらを組み合わせて指定することで、GRIB2形式ファイルから使用したい気象データの種類と高度を取り出すことができます。

```python:pygrib-select.py

#- parameterName for select from GRIB file

letter_mslp = "Pressure reduced to MSL"
letter_tmp  = "Temperature"
letter_rh   = "Relative humidity"
letter_vv   = "Vertical velocity [pressure]"
letter_wu   = "u-component of wind"
letter_wv   = "v-component of wind"
letter_gh   = "Geopotential height"
letter_lc   = "Low cloud cover"
letter_mc   = "Medium cloud cover"
letter_hc   = "High cloud cover"
letter_tc   = "Total cloud cover"

#-- Surface

grb_tmp = grbs.select(parameterName=letter_tmp  , level=2 , forecastTime=0)
grb_prs = grbs.select(parameterName=letter_mslp , level=0 , forecastTime=0)
grb_wu  = grbs.select(parameterName=letter_wu   , level=10, forecastTime=0)
grb_wv  = grbs.select(parameterName=letter_wv   , level=10, forecastTime=0)
grb_rh  = grbs.select(parameterName=letter_rh   , level=2 , forecastTime=0)

#-- 850 hPa level

grb_tmp85 = grbs.select(parameterName=letter_tmp  , level=850 , forecastTime=0)
grb_wu85  = grbs.select(parameterName=letter_wu   , level=850 , forecastTime=0)
grb_wv85  = grbs.select(parameterName=letter_wv   , level=850 , forecastTime=0)
grb_rh85  = grbs.select(parameterName=letter_rh   , level=850 , forecastTime=0)
grb_gh85  = grbs.select(parameterName=letter_gh   , level=850 , forecastTime=0) 

#-- 700 hPa level 

grb_tmp70 = grbs.select(parameterName=letter_tmp  , level=700 , forecastTime=0)
grb_rh70  = grbs.select(parameterName=letter_rh   , level=700 , forecastTime=0)
grb_vv70  = grbs.select(parameterName=letter_vv   , level=700 , forecastTime=0)

#-- 500 hPa level

grb_tmp50 = grbs.select(parameterName=letter_tmp  , level=500 , forecastTime=0)
grb_wu50  = grbs.select(parameterName=letter_wu   , level=500 , forecastTime=0)
grb_wv50  = grbs.select(parameterName=letter_wv   , level=500 , forecastTime=0)
grb_rh50  = grbs.select(parameterName=letter_rh   , level=500 , forecastTime=0)
grb_gh50  = grbs.select(parameterName=letter_gh   , level=500 , forecastTime=0)

#-- 300 hPa level

grb_wu30  = grbs.select(parameterName=letter_wu   , level=300 , forecastTime=0)
grb_wv30  = grbs.select(parameterName=letter_wv   , level=300 , forecastTime=0)
grb_rh30  = grbs.select(parameterName=letter_rh   , level=300 , forecastTime=0) 

#-- Cloud coverage

grb_lc  = grbs.select(parameterName=letter_lc   , forecastTime=0)
grb_mc  = grbs.select(parameterName=letter_mc   , forecastTime=0)
grb_hc  = grbs.select(parameterName=letter_hc   , forecastTime=0)
grb_tc  = grbs.select(parameterName=letter_tc   , forecastTime=0)

```

例えば850hPa面の気温データは`grb_tmp85.values`という形でアクセスすることができます。


## 3. MetPy

### 3.1 MetPyについて
Metpyは気象関係の様々なツールを持っているプロジェクトです。

[Metpy ホームページ](https://unidata.github.io/MetPy/latest/index.html “Metpy”)

May, R. M., Arms, S. C., Marsh, P., Bruning, E., Leeman, J. R., Goebbert, K., Thielen, J. E.,
    and Bruick, Z., 2020: MetPy: A Python Package for Meteorological Data.
    Version 1.0.0rc1, Unidata, Accessed 14 January 2020.
    [Available online at https://github.com/Unidata/MetPy.]
    doi:10.5065/D6WW7G29.

### 3.2 相当温位などの計算

**相当温位**は、GPVから直接得られる**気温**、**相対湿度**の値から計算によって導き出さねばなりませんが、MetPyのcalcに含まれる以下の計算ライブラリを用いることで簡単に計算することが出来ます。

`dewpoint_rh` **気温**と**相対湿度**から**露点温度**を計算する
`equivalent_potential_temperature`　**露点温度**、**気圧**、**気温**から**相当温位**を計算する。

## 4.地図のあれこれ
### 4.1 投影法
速報天気図から前線位置を学習するので、可視化画像の地図は、速報天気図と形式・描画範囲を合わせておく必要があります。

<details><summary>合わせておく必要があるか？</summary>
というのは言い過ぎかもしれません。私の手法ではその必要があった、ということです。範囲が違っても沢山のデータから適正な前線を描画する機械学習を作成することができるかもしれません。
</details>

地球が球体なので、平面の地図にするには様々な手法があるというのはご存知と思いますが、速報天気図もよくみると高緯度側ほど狭く表記されています。

速報天気図（SPAS)は、メルカトル図法ではなくポーラーステレオ投影法によっています。
気象庁が提供している、下記の技術情報には、
**「地図の投影方法は、北緯 60 度、東経 140 度を基準としたポーラーステレオ図法を用いている。」**
との記載があります。
[配信資料に関する技術情報(気象編)第358号](http://www.data.jma.go.jp/add/suishin/jyouhou/pdf/358.pdf "配信資料に関する技術情報(気象編)第 358 号")
この後matplotlibのBasemapを用いてデータを地図に重ね書きする際に、この情報を使用します。

### 4.2 描画範囲

これはもう速報天気図から読み取るしかありませんでした。
Basemapのパラメタとして指定するのは左下、右上の緯度経度なのです。
実際のところどのように画像化されるかを**トライアンドエラー**で指定する値を探り出しました。

結果は下記になります。ポーラステレオの基点は、lat_0, lon_0 で指定します。
左下の緯度経度は（115°E, 9°N）、右上は（178°E, 54°N）としています。

```python:
Basemap(projection='stere', llcrnrlat=9, urcrnrlat=54, llcrnrlon=115, urcrnrlon=178, lat_0=60,  lon_0=140, resolution='i' )
```
これをもとに12月2日12時UTCの地上気圧（海面更正気圧）を描画しました。
![gsm_prs_srf_2019120212.v8.512.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/f86f2c64-4ed7-cb58-75b3-31d25e4bd78d.png)
速報天気図は下記です。
![SPAS_COLOR_201912021200.v8.512.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/f43eb766-c3b9-4434-8c0b-10e83fd378df.png)

GSMは本来は全球データなので地球上全てをカバーするだけのデータがあるのですが、使用するのはこの範囲だけです。少しもったいない気もします。
可視化した地図ですが、天気図と比べると端の方が少し気になるところもありますが、まあこれでよしとしました。

## 5.matplotlibによる可視化

### 5.1 可視化について

matplotlibについてはもう解説は不要だと思います。

気象データは地表面だったり、高層のある気圧等圧面だったりの上で（緯度、経度）に並んでいる数値データです。
これを人間が見て分かりやすいというより、CNNが学習しやすいように可視化します。

２次元のデータの可視化としてはコンターマップや等圧線図があります。
天気予報だと気温のカラー範囲を夏と冬でひっそりと変えたりしていますが、今回は学習することを考えて、年間を通して同じ範囲とします。

```python:value_range.py
levels_prs    = np.arange(930.0,1080.0,4.0) #海面更正気圧　  930-1080hPa 4hPa間隔
levels_tmp    = np.arange(210,316,2)        #気温(ケルビン） 210-316K 2K間隔
levels_tmps   = np.arange(-35,45,2)         #気温(摂氏）    -35°Cから45°C 2°C間隔
levels_dp     = np.arange(0,80,10)          #露点温度(摂氏)  可視化には使わない
levels_dp2    = np.arange(0,3,1)            #湿数(摂氏)     0°Cから3°C 1°C間隔
levels_vv     = np.arange(-5,5,0.5)         #鉛直速度(P/s)  -5P/s-5P/s 0.5P/s間隔
levels_ept    = np.arange(200,400,3)        #相当温位(ケルビン) 200-400K 3K間隔
levels_cld    = np.arange(-20,110,5)        #雲量(%)        0-110
levels_rh     = np.arange(0.0,110.0,5)      #相対湿度(%)     0-110 5%間隔
levels_gh_300 = np.arange(7000,10000,40)    #ジオポテンシャル高度 7000-10000m 40m間隔
levels_gh_500 = np.arange(4500,6200,40)     #ジオポテンシャル高度 4500- 6200m 40m間隔
levels_gh_850 = np.arange(1000,2000,40)     #ジオポテンシャル高度 1000- 2000m 40m間隔
levels_gh_700 = np.arange(2000,4000,40)     #ジオポテンシャル高度 2000- 4000m 40m間隔
```

### 5.2 矢羽・コンター・等値線での可視化

今回CNNの入力とするのは、以下の種類です。

1. 地表面　　**風　気温　気圧**
2. 850hPa面    **風　気温  気圧**（**ジオポテンシャル高度**）
3. 500hPa面    **風　気温  気圧**（**ジオポテンシャル高度**） 
風は風速と風向を矢羽図で表現します（matplotlibのbarbsメソッドを利用します）。
気圧はコンター線図、気温はコンター線図とカラーコンターを重ね書きします。
気圧は上空ではジオポテンシャル高度という要素で表されます。気象の世界では、高層の気象要素は、ある気圧となる高さの分布として表します。気圧は高度が上がるに従って指数関数的に下がっていきます。高気圧や低気圧に従って、指定した気圧面の高さは高くなったり低くなったりします。これをジオポテンシャル高度と呼びます。
**前線**となる候補は、気温が急激に変化している領域や、風向が急激に変化している領域です。
これら３種類のデータを重ねて描くソースは下記に示します。

```Python:windmap

def windmap( u , v , values , levels , cmap , _save_filename , values2 , levels2 ):
        # u, v 風速東西成分、風速南北成分
        # values  気圧データ
        # values2 気温データ
        # levels, levels2 カラーリング範囲
        # cmap    カラーマップ
        # _save_filename  保存ファイル名

        fig,ax = plt.subplots(figsize=(6,5.8))　# 図を作成する
        plt.subplots_adjust(left=0.01, right=0.98, top=0.99, bottom=0.01)　#余白を調整

        # Basemapでポーラステレオ投影地図を作成する
        m = Basemap(projection='stere', llcrnrlat=9, urcrnrlat=54, llcrnrlon=115, urcrnrlon=178, lat_0=60,  lon_0=140, resolution='i' )
        m.drawparallels(np.arange(-80.,81.,10.))      # 緯度線を10°刻みで描く
        m.drawmeridians(np.arange(-180.,181.,10.))    # 経度線を10°刻みで描く
        m.drawcoastlines()                            # 海岸線を描く
        
        lons_s = lons[::5, ::5]   # 風矢羽描画ポイントを５要素飛ばしとする
        lats_s = lats[::5, ::5]   # 風矢羽描画ポイントを５要素飛ばしとする
        
        x , y    = m(lons, lats)     # 緯度経度から、Basemapの描画座標へ変換する
        xs , ys  = m(lons_s, lats_s) # 緯度経度から、Basemapの描画座標へ変換する
        us = u[::5, ::5]             # 風データ（東西成分）を５要素飛ばしとする
        vs = v[::5, ::5]             # 風データ（南北成分）を５要素飛ばしとする
        
        m.contour(  x , y , values  , levels=levels  , linewidths=0.7 , colors='k' )
        # valuesの等値線を作成する
        m.contour(  x , y , values2 , levels=levels2 , linewidths=0.3 , colors='k' )
        # values2の等値線を作成する
        m.contourf( x , y , values2 , levels2 , cmap=cmap )
        # カラーのコンターマップ(values2) を上書きで作成する
        m.barbs( xs , ys , us , vs  , length=4.5 )
        # 矢羽図を上書きで作成する
        fig.savefig(_save_filename)
```

### 5.3 コンターマップでの可視化

1. 850hPa面　　**相当温位**
まず**温位**という量があって、これは注目している空気を1000hPaまで持ってきた時の温度です。
空気が下降して断熱圧縮すると温度が上がりますし、上昇して断熱膨張すると温度が下がります。同じ気圧で比べないと、その空気が持っている熱的なエネルギーを比較しづらいためにこのような概念の量が使用されます。
さらに、空気が水蒸気を含んでいると、気温が下がって凝結する際には潜熱を放出して空気の温度が上がります。そこで水蒸気を凝結させることも考慮した量を**相当温位**と呼びます。
天気予報で梅雨の終わりの頃など、「南から暖かく湿った空気が流れ込み前線を刺激するために非常に激しい雨が降るでしょう」などと聞く時の「暖かさと湿り具合」を同時に評価する指標です。
**前線**となる候補ですが、相当温位が急に変わっている領域は、前線候補の有望領域です。

2. 700hPa面　　**湿数**
湿数は、気温と露点温度の差です。露点温度は水蒸気が凝結する温度ですので、この**湿数**が小さいほど空気が凝結しやすい、つまり水蒸気量が飽和水蒸気量に近くなっていて、空気が湿っていることを示します。
気象庁が発表している高層気象図では、湿数が3°C以下の部分に編みかけがされていて、**湿潤域**とされています。
低気圧中心付近や**温暖前線前面**では空気が湿潤で雨が降っていることが多いです。

3. 700hPa面　　**鉛直風**の風速
鉛直方向の風速は、空気の上下方向への動きを表しています。空気が集まってきて行き場がなくなったり、山の斜面にぶつかったりすると、上昇することがあります。上昇すると気温が下がり、露点温度になると空気に含まれている水蒸気が凝結します。これは降水活動の始まりになりますので、鉛直風は降水活動と密接な関係になります。
また、**温帯低気圧の発達には前面で温暖前線面に沿った空気の上昇、後面で寒冷前線に沿った空気の下降が必要となります**ので、前線を含む気象現象には鉛直風が大きく関わっています。

これらは下記のソースにより可視化されます。先ほどに比べると、矢羽図ともう一つのコンターマップが存在しないだけです。

```Python:drawmap
def drawmap( values , levels , cmap , _save_filename ) :
        fig,ax = plt.subplots(figsize=(6,5.8))
        plt.subplots_adjust(left=0.01, right=0.98, top=0.99, bottom=0.01)
        m = Basemap(projection='stere', llcrnrlat=9, urcrnrlat=54, llcrnrlon=115, urcrnrlon=178, lat_0=60,  lon_0=140, resolution='i' )
        
        m.drawparallels(np.arange(-80.,81.,10.))
        m.drawmeridians(np.arange(-180.,181.,10.))
        m.drawcoastlines()
        x , y = m(lons, lats)
        #
        m.contour(  x , y , values , levels=levels , linewidths=0.5 , colors='k' )
        m.contourf( x , y , values , levels , cmap=cmap )
        #
        fig.savefig(_save_filename)
```

### 5.4 GSM可視化データ作成頻度

1日にGSMは4回初期値を作ります。時刻はUTCでいうと、0時、6時、12時および18時です。
一方で速報天気図は、0時を除いて3時間ごとに作成されます。
つまり、同じ時刻で初期値と天気図が揃うのは、6時, 12時, 18時の３回です。
これを3年分、入力と教師データを作ることにしました。

2019年の1月分と8月分を動画化したものを下記に示します。
データは京都大学生存圏研究所が公開しているGPVデータベースからダウンロードしました。

[京都大学生存圏研究所 生存圏データベース](http://database.rish.kyoto-u.ac.jp/arch/jmadata/gpv-original.html "京都大学生存圏研究所 生存圏データベース　気象庁データ")

カラーリングレンジは変えていないので、夏と冬の違いがよくわかるかと思います。
普通の可視化図と違って、学習のノイズになるので、凡例もつけていません。

地上風・気圧・気温 2019年1月（赤いほど暖かい、青いほど寒い）
![ww_wndsrf01.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/37122a41-af46-2fbc-6f67-c4f80628f86b.gif)
地上風・気圧・気温 2019年8月（赤いほど暖かい、青いほど寒い）
![ww_wndsrf08.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/526166d6-45fa-5c30-d97e-588fc3b70afd.gif)
850hPa 風・気圧・気温 2019年1月（赤いほど暖かい、青いほど寒い）
![ww_wnd8501.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/e3dcdedf-d014-a894-dc7f-e66c360cb56e.gif)
850hPa 風・気圧・気温 2019年8月（赤いほど暖かい、青いほど寒い）
![ww_wnd8508.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/9b23179c-f44b-2f77-6f80-fe5d07e25746.gif)
500hPa 風・気圧・気温 2019年1月（赤いほど暖かい、青いほど寒い）
![ww_wnd5001.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/420e84e7-b5d0-3399-1bb5-e2c2bf675f89.gif)
500hPa 風・気圧・気温 2019年8月（赤いほど暖かい、青いほど寒い）
![ww_wnd5008.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/eb7bcaf6-3ddd-76c0-c350-57b95e81a5af.gif)
850hPa 相当温位 2019年1月（赤いほど暖湿）
![ww_ept8501.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/a57e793c-c3a4-24fe-c481-3d9e015452c0.gif)
850hPa 相当温位 2019年8月（赤いほど暖湿）
![ww_ept8508.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/c2f3c0af-4fb9-9036-71c5-6366a6ab20dc.gif)
700hPa 湿数 2019年1月（青いほど湿数小＝湿潤）
![ww_situ7001.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/49343eca-7183-2c53-3398-2f6860e01494.gif)
700hPa 湿数 2019年8月（青いほど湿数小＝湿潤）
![ww_situ7008.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/6d5dbcaf-bb5d-25ac-00cd-5ee3d56d5f2a.gif)
700hPa 鉛直風 2019年1月（青いほど上昇流、赤いほど下降流）
![ww_vvl7001.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/691ebb9d-3ecb-b439-4102-e8aaabc14d6d.gif)
700hPa 鉛直風 2019年8月（青いほど上昇流、赤いほど下降流）
![ww_vvl7008.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/f7c6a29e-1d99-6fbb-5dc6-69373f961631.gif)

速報天気図 2019年1月
![ww201901.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/14b35028-466c-7690-92fc-082bfde45d70.gif)
速報天気図 2019年8月
![ww201908.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/411112/c517aaf1-3816-86fb-8f74-adc953241d8c.gif)


## まとめ

今回はとても長くなってしまいましたが、機械学習の入力となる気象データの可視化までを投稿しました。
次回は、機械学習の教師データとなる、前線を速報天気図からどのように取り出したかについて投稿する予定です。

[次回：気象データをもとに「天気図っぽい前線」を機械学習で描いてみる(3)](https://qiita.com/m-taque/items/4d5bb45e6b5dc42dc833 "気象データをもとに「天気図っぽい前線」を機械学習で描いてみる(3)")


## 補遺

### parameterNameをGRIB2ファイルからどのように特定するか

少々泥臭いやり方ではありますが、GRIB2ファイルをopenして、含まれているデータの一覧を出します。

```python
>>> import pygrib
>>> grbs=pygrib.open("Z__C_RJTD_20181209000000_GSM_GPV_Rgl_FD0000_grib2.bin")
>>> for grb in grbs:
...     grb
... 
1:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 0 hrs:from 201812090000
2:Surface pressure:Pa (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 201812090000
3:10 metre U wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 10 m:fcst time 0 hrs:from 201812090000
4:10 metre V wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 10 m:fcst time 0 hrs:from 201812090000
5:2 metre temperature:K (instant):regular_ll:heightAboveGround:level 2 m:fcst time 0 hrs:from 201812090000
6:2 metre relative humidity:% (instant):regular_ll:heightAboveGround:level 2 m:fcst time 0 hrs:from 201812090000
7:Low cloud cover:% (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 201812090000
8:Medium cloud cover:% (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 201812090000
9:High cloud cover:% (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 201812090000
10:Total cloud cover:% (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 201812090000
11:Geopotential Height:gpm (instant):regular_ll:isobaricInhPa:level 100000.0 Pa:fcst time 0 hrs:from 201812090000
12:U component of wind:m s**-1 (instant):regular_ll:isobaricInhPa:level 100000.0 Pa:fcst time 0 hrs:from 201812090000
13:V component of wind:m s**-1 (instant):regular_ll:isobaricInhPa:level 100000.0 Pa:fcst time 0 hrs:from 201812090000

(中略）

98:Vertical velocity:Pa s**-1 (instant):regular_ll:isobaricInhPa:level 2000.0 Pa:fcst time 0 hrs:from 201812090000
99:Geopotential Height:gpm (instant):regular_ll:isobaricInhPa:level 1000.0 Pa:fcst time 0 hrs:from 201812090000
100:U component of wind:m s**-1 (instant):regular_ll:isobaricInhPa:level 1000.0 Pa:fcst time 0 hrs:from 201812090000
101:V component of wind:m s**-1 (instant):regular_ll:isobaricInhPa:level 1000.0 Pa:fcst time 0 hrs:from 201812090000
102:Temperature:K (instant):regular_ll:isobaricInhPa:level 1000.0 Pa:fcst time 0 hrs:from 201812090000
103:Vertical velocity:Pa s**-1 (instant):regular_ll:isobaricInhPa:level 1000.0 Pa:fcst time 0 hrs:from 201812090000
```

ここで例えば100番レコード(GRIBでmessagesというようです）の``paramaterName``を取得するには、

```python
>>> f=grbs[100]
>>> f['parameterName']
'u-component of wind'
```

のように、keyに``parameterName``を指定することで参照することができます。

[pygrib documentation](https://jswhit.github.io/pygrib/docs/)
を見ると、本来は一覧を出した後の最初のフィールドが``name``というパラメタですので、``grb = grbs.select(name='U component of wind')[0]``のように取り出すのが使い方のようです。

ところがある時、nameがセットされていないmessageに遭遇したことがあり、色々探した結果parameterNameに行き着きました。これがセットされていないmessageは、私は今のところ遭遇したことはありません。

余談ですが、各messageが持つkeyの一覧は、下記のように取り出すことができます。

```python
>>> f.keys()
['globalDomain', 'GRIBEditionNumber', 'tablesVersionLatest', 'grib2divider', 'is_efas', 'angleSubdivisions', 'missingValue', 'ieeeFloats', 'isHindcast', 'section0Length', 'identifier', 'discipline', 'editionNumber', 'totalLength', 'sectionNumber', 'section1Length', 'numberOfSection', 'centre', 'centreDescription', 'subCentre', 'tablesVersion', 'masterDir', 'localTablesVersion', 'significanceOfReferenceTime', 'year', 'month', 'day', 'hour', 'minute', 'second', 'dataDate', 'julianDay', 'dataTime', 'productionStatusOfProcessedData', 'typeOfProcessedData', 'md5Section1', 'selectStepTemplateInterval', 'selectStepTemplateInstant', 'stepType', 'is_chemical', 'is_chemical_distfn', 'is_aerosol', 'is_aerosol_optical', 'setCalendarId', 'deleteCalendarId', 'is_uerra', 'sectionNumber', 'grib2LocalSectionPresent', 'deleteLocalDefinition', 'sectionNumber', 'gridDescriptionSectionPresent', 'section3Length', 'numberOfSection', 'sourceOfGridDefinition', 'numberOfDataPoints', 'numberOfOctectsForNumberOfPoints', 'interpretationOfNumberOfPoints', 'PLPresent', 'gridDefinitionTemplateNumber', 'gridDefinitionDescription', 'shapeOfTheEarth', 'scaleFactorOfRadiusOfSphericalEarth', 'scaledValueOfRadiusOfSphericalEarth', 'scaleFactorOfEarthMajorAxis', 'scaledValueOfEarthMajorAxis', 'scaleFactorOfEarthMinorAxis', 'scaledValueOfEarthMinorAxis', 'radius', 'Ni', 'Nj', 'basicAngleOfTheInitialProductionDomain', 'mBasicAngle', 'angleMultiplier', 'mAngleMultiplier', 'subdivisionsOfBasicAngle', 'angleDivisor', 'latitudeOfFirstGridPoint', 'longitudeOfFirstGridPoint', 'resolutionAndComponentFlags', 'resolutionAndComponentFlags1', 'resolutionAndComponentFlags2', 'iDirectionIncrementGiven', 'jDirectionIncrementGiven', 'uvRelativeToGrid', 'resolutionAndComponentFlags6', 'resolutionAndComponentFlags7', 'resolutionAndComponentFlags8', 'ijDirectionIncrementGiven', 'latitudeOfLastGridPoint', 'longitudeOfLastGridPoint', 'iDirectionIncrement', 'jDirectionIncrement', 'scanningMode', 'iScansNegatively', 'jScansPositively', 'jPointsAreConsecutive', 'alternativeRowScanning', 'iScansPositively', 'scanningMode5', 'scanningMode6', 'scanningMode7', 'scanningMode8', 'g2grid', 'latitudeOfFirstGridPointInDegrees', 'longitudeOfFirstGridPointInDegrees', 'latitudeOfLastGridPointInDegrees', 'longitudeOfLastGridPointInDegrees', 'iDirectionIncrementInDegrees', 'jDirectionIncrementInDegrees', 'latLonValues', 'latitudes', 'longitudes', 'distinctLatitudes', 'distinctLongitudes', 'gridType', 'md5Section3', 'sectionNumber', 'section4Length', 'numberOfSection', 'NV', 'neitherPresent', 'productDefinitionTemplateNumber', 'genVertHeightCoords', 'parameterCategory', 'parameterNumber', 'parameterUnits', 'parameterName', 'typeOfGeneratingProcess', 'backgroundProcess', 'generatingProcessIdentifier', 'hoursAfterDataCutoff', 'minutesAfterDataCutoff', 'indicatorOfUnitOfTimeRange', 'stepUnits', 'forecastTime', 'startStep', 'endStep', 'stepRange', 'stepTypeInternal', 'validityDate', 'validityTime', 'typeOfFirstFixedSurface', 'unitsOfFirstFixedSurface', 'nameOfFirstFixedSurface', 'scaleFactorOfFirstFixedSurface', 'scaledValueOfFirstFixedSurface', 'typeOfSecondFixedSurface', 'unitsOfSecondFixedSurface', 'nameOfSecondFixedSurface', 'scaleFactorOfSecondFixedSurface', 'scaledValueOfSecondFixedSurface', 'pressureUnits', 'typeOfLevel', 'level', 'bottomLevel', 'topLevel', 'tempPressureUnits', 'paramIdECMF', 'paramId', 'shortNameECMF', 'shortName', 'unitsECMF', 'units', 'nameECMF', 'name', 'cfNameECMF', 'cfName', 'cfVarNameECMF', 'cfVarName', 'modelName', 'ifsParam', 'PVPresent', 'deletePV', 'md5Section4', 'lengthOfHeaders', 'md5Headers', 'sectionNumber', 'section5Length', 'numberOfSection', 'numberOfValues', 'dataRepresentationTemplateNumber', 'packingType', 'referenceValue', 'referenceValueError', 'binaryScaleFactor', 'decimalScaleFactor', 'optimizeScaleFactor', 'bitsPerValue', 'typeOfOriginalFieldValues', 'md5Section5', 'sectionNumber', 'section6Length', 'numberOfSection', 'bitMapIndicator', 'bitmapPresent', 'md5Section6', 'sectionNumber', 'section7Length', 'numberOfSection', 'codedValues', 'values', 'packingError', 'unpackedError', 'maximum', 'minimum', 'average', 'numberOfMissing', 'standardDeviation', 'skewness', 'kurtosis', 'isConstant', 'changeDecimalPrecision', 'decimalPrecision', 'setBitsPerValue', 'getNumberOfValues', 'scaleValuesBy', 'offsetValuesBy', 'productType', 'md5Section7', 'section8Length', 'analDate', 'validDate']
```

今回の前線描画では使用していませんが、MSMを可視化したときに使用したparameterNameです。

| 指定パラメタ | 内容                   |
|:----------|:-----------------------|
|Pressure reduced to MSL | 海面更正気圧|
|Temperature | 気温 |
|Relative humidity|相対湿度|
|Vertical velocity [pressure]|鉛直速度（気圧座標）|
|u-component of wind|風速の東西成分|
|v-component of wind|風速の南北成分|
|Geopotential height|ジオポテンシャル高度|
|Total precipitation|降水量|
|Downward short-wave radiation flux|下向き短波放射|
| Low cloud cover      | 下層雲雲量         |
| Medium cloud cover   | 中層雲雲量        |
| High cloud cover     | 上層雲雲量　　　　|
| Total cloud cover    | 全雲量　|