Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
Help us understand the problem. What is going on with this article?

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

More than 1 year has passed since last update.

気象業界では、数値予報や解析値などの面的なデータにGRIB2というバイナリフォーマットが使われます。@mhangyoさんのこちらの記事ではMatplotlib Basemap Toolkitを用いてGRIB2を可視化していますが、Basemapは開発終了となるため、ここではCartopyを用いた可視化に変えてみます。

(2019-02-28に書いた記事です。)

はじめに:Matplotlib Basemap ToolkitからCartopyへ

pandas等で読み込んだデータを地図にプロットしたいとき、Matplotlib Basemap Toolkitが広く使われていると思うのですが、最近はCartopyへの移行が強く推奨されています。以下に、公式サイトのアナウンスの訳(一部意訳)を載せておきます。

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を拡張しています。今回は、以下のような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さんの記事と同じMSM(Meso-Scale Model, 気象庁メソモデル)のデータファイルを用意しておきます。

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

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

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

import

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

%matplotlib inline

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

読み込み

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

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単位に変換しておきましょう。

mslp_hpa = mslp / 100

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

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

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

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

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

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

正射図法でのプロット

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

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

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

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

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

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

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

FROM jupyter/scipy-notebook

USER root

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

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

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を再掲します。

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が通ります。

import pygrib

print(pygrib.__version__)

出力:

2.0.3

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

noritada
共訳書として『Pythonによるデータ分析入門』(オライリー・ジャパン)、監訳書として『Subversion実践入門 達人プログラマに学ぶバージョン管理 第2版』(オーム社)があります。
weathernews
ウェザーニューズは、気象リスクを軽減するための気象情報サービスを行なっています。そのための気象予測やその情報活用に関するサービス開発や研究を行っています。
https://weathernews.com
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away