GRIB2をCartopyで地図上に可視化してみる

気象業界では、数値予報や解析値などの面的なデータにGRIB2というバイナリフォーマットが使われます。@mhangyoさんの[こちらの記事](https://qiita.com/mhangyo/items/f06debce3975a269a658)では[Matplotlib Basemap Toolkit](https://matplotlib.org/basemap/index.html)を用いてGRIB2を可視化していますが、Basemapは開発終了となるため、ここでは[Cartopy](https://scitools.org.uk/cartopy/docs/latest/)を用いた可視化に変えてみます。

（2019-02-28に書いた記事です。）

## はじめに：Matplotlib Basemap ToolkitからCartopyへ

pandas等で読み込んだデータを地図にプロットしたいとき、Matplotlib Basemap Toolkitが広く使われていると思うのですが、最近はCartopyへの移行が強く推奨されています。以下に、[公式サイトのアナウンス](https://matplotlib.org/basemap/users/intro.html#cartopy-new-management-and-eol-announcement)の訳（一部意訳）を載せておきます。

> ### Cartopyとプロジェクト管理の新たなフェーズ、そしてサポート終了のお知らせ
>
> 2016年以降、Basemapのソフトウェアプロジェクト管理は新たなフェーズに入りました。CartopyというプロジェクトがBasemapを置き換えていく予定です。CartopyはまだBasemapのすべての機能を実装しているわけではありません。しかし、新しいソフトウェア開発においては可能な限りCartopyを使ったほうがよいですし、Basemapを使っている既存のソフトウェアはCartopyへの移行プロセスを開始したほうがよいでしょう。すべての保守や開発の取り組みはCartopyに集中させるべきです。
>
> Ben Rootさんが2020年までのBasemapの保守を進んで引き受けてくれましたので、プルリクエストの査読やバグ修正は行われる予定です。また、NumPyやMatplotlibなどのパッケージとの互換性を保証するための保守もなされます。新機能の追加も受け入れます。しかし、強調しておきますが、新たな開発の取り組みはCartopyに集中させるべきです。というのも、Python 2.7の公式サポートが2020年に終了するときに、Ben RootさんはBasemapのリリースをしてサポートを終えるつもりだからです。

## 環境準備

まず、pygribとCartopyの入った環境を用意しましょう。

自分の場合、最近、環境作りではDockerを用いることが多く、特にJupyterを用いたデータサイエンス環境を作りたい場合はたいてい[scipy-notebook](https://hub.docker.com/r/jupyter/scipy-notebook/)を拡張しています。今回は、以下のようなDockerfileで行ないました。ハマったところもあるので、おまけとして記事の最後に書いてあります。

```dockerfile
FROM jupyter/scipy-notebook

USER root

RUN conda install --quiet --yes \
    'pygrib' \
    'cartopy' && \
    conda clean -tipsy && \
    fix-permissions $CONDA_DIR

RUN pip install --upgrade numpy
```

データは、[@mhangyoさんの記事](https://qiita.com/mhangyo/items/f06debce3975a269a658)と同じMSM（Meso-Scale Model, 気象庁メソモデル）のデータファイルを用意しておきます。

## GRIB2の読み込みから可視化まで

Pythonでの読み込みと可視化に入っていきます。

以下のコードは、先程のDockerfileからビルドしたDockerイメージのJupyterで動作を確認していますが、matplotlibの出力先などを適切に設定してあげれば他の環境でも動くはずです。

### import

まずはimportです。pygribとCartopy以外はよく見るimport文です。

```python

%matplotlib inline

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import pygrib
import seaborn as sns
```

### 読み込み

GRIB2を読み込んで、緯度、経度、格子点値を変数に格納します。読み込みなので、当然ながらBasemap使用時からの変化は特にありません。

```python
grbs = pygrib.open('Z__C_RJTD_20170102000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin')
grb = grbs.select(forecastTime=0)[0]
lats, lons = grb.latlons()
mslp = grb.values
```

### 可視化の準備

続いて、可視化のためにデータを調整していきます。

今回は、等圧線にラベルも入れることにします。気圧の値は、GRIB2内にはPa単位で入っていますが、気象業界ではhPa単位で表すのが通例です。今回もラベルはhPa単位で表示したいので、Pa単位の格子点値をすべてhPa単位に変換しておきましょう。

```python
mslp_hpa = mslp / 100
```

また、等高線はhPaが2の倍数となるところに引きたいので、2の倍数のみをレベルとするように調整します。

```python
levels = np.arange(mslp_hpa.min()//2*2, mslp_hpa.max()//2*2+2, 2)
```

### 正距円筒図法でのプロット

まずは、Basemapのデフォルト（`projection='cyl'`）と同じ正距円筒図法でプロットしてみましょう。正距円筒図法は、緯度と経度が直角かつ等間隔に交差する最もシンプルな図法で、緯度と経度の格子データとして表される気象データとは相性がよい図法でもあります。

```python
lat_s = lats.min()
lat_n = lats.max()
lon_w = lons.min()
lon_e = lons.max()

proj = ccrs.PlateCarree() # 正距円筒図法を指定

ax = plt.axes(projection=proj)
ax.gridlines()
ax.coastlines(resolution='10m')
ax.set_extent((lon_w, lon_e, lat_s, lat_n), proj) # 緯度経度の範囲を指定

CS = ax.contour(lons, lats, mslp_hpa, levels, transform=proj)
ax.clabel(CS, fmt='%.f') # 気圧を整数で表示するよう、fmtを指定
```

![PlateCarree.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/39513/345150fd-e540-4952-4a15-0b82d57d40f3.png)


正距円筒図法は気象データとの相性はよいのですが、高緯度の地域で潰れた感じになるのが少し気になりますね。

### 正射図法でのプロット

続いて、正射図法で表示してみます。正射図法は、地球の表面をそのまま真横の平面へと投影した図法です。遥か彼方から地球を眺めたときの半球面を描くことになるため、自然な感じになります。

正距円筒図法は、すべての緯線と経線を平等に扱うため、特にパラメータの指定は必要ありませんでした。しかし正射図法では、「どこを中心に投影するか」、すなわち「どこを中心に眺めるか」で見たい地域の形が変わってくるため、投影の中心の緯度経度を指定してあげる必要があります。例えば、日本の中心を投影の中心とすれば日本はきれいに描けますが、太平洋の中心を投影の中心とすると日本はかなり歪んでしまいます。

ここでは、描画範囲の中心を投影の中心にしてプロットしてみましょう。

```python
lat_s = lats.min()
lat_n = lats.max()
lon_w = lons.min()
lon_e = lons.max()

proj = ccrs.Orthographic((lon_w + lon_e) / 2, (lat_s + lat_n) / 2) # 対象緯度範囲、経度範囲の中心を投影の中心とした正射図法を指定
latlon_proj = ccrs.PlateCarree() # 緯度経度の処理用に正距円筒図法も使う

ax = plt.axes(projection=proj) # 地図は正射図法で表示
ax.gridlines()
ax.coastlines(resolution='10m')
ax.set_extent((lon_w, lon_e, lat_s, lat_n), latlon_proj) # 緯度経度で指定した描画範囲は正距円筒図法で処理

CS = ax.contour(lons, lats, mslp_hpa, levels, transform=latlon_proj) # 緯度経度の座標値は正距円筒図法で処理
ax.clabel(CS, fmt='%.f')
```

![Orthographic.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/39513/1252f970-f3fa-1dc1-201c-fefcf3f1d026.png)

一つの描画用のソースコードで図法が2つ登場するので、最初は混乱しやすいかもしれません。
**「緯度経度を指定する箇所では、正距円筒図法（`cartopy.crs.PlateCarree`）を必ず使う」**と考えるとよいかと思います。

## おまけ：Dockerでのpygrib環境作りでハマったところの共有

今回、pygribを入れるため、最初は以下のようなDockerfileを用意しました。

```dockerfile
FROM jupyter/scipy-notebook

USER root

RUN conda install --quiet --yes \
    'pygrib' \
    'cartopy' && \
    conda clean -tipsy && \
    fix-permissions $CONDA_DIR
```

ところが、イメージのビルドは成功したように見えるものの、立ち上げると、以下のようにpygribのimportができません。

```python
import pygrib

print(pygrib.__version__)
```

出力：

```
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
RuntimeError: module compiled against API version 0xc but this version of numpy is 0xb
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-1-7b90afe9041e> in <module>
----> 1 import pygrib
      2 
      3 print(pygrib.__version__)

ImportError: numpy.core.multiarray failed to import
```

Anacondaでこのようなことが発生した場合、NumPyを入れ直してあげるとよいようです。本記事冒頭のDockerfileを再掲します。

```dockerfile
FROM jupyter/scipy-notebook

USER root

RUN conda install --quiet --yes \
    'pygrib' \
    'cartopy' && \
    conda clean -tipsy && \
    fix-permissions $CONDA_DIR

RUN pip install --upgrade numpy
```

このイメージだと、きちんとimportが通ります。

```python
import pygrib

print(pygrib.__version__)
```

出力：

```
2.0.3
```

一つの記事にするほど深く調べたわけでもないので、ついでとしてこちらに書いておきました。