2023年度 WXBC 人材育成WG テクノロジー研修
気象庁GPVデータ分析チャレンジ！入門¶
目次：¶
- はじめに
- 1 GPV データ
- 2 WXBCオリジナルPythonライブラリwxbcgribX
- 3 GPVデータの読み込み
- 4 GPVデータの処理
- 4.1 オブジェクト xarray.DataArray
- 4.2 GPV構成要素の取り出し
- 4.3 気象データへのアクセス
- (1) Fancy indexing-Fancy-indexing)
- (2) メソッド isel-メソッド-isel)
- (3) メソッド sel-メソッド-sel)
- 4.4 気象データの演算
- (1) 絶対温度からセルシウス度への変更-絶対温度からセルシウス度への変更)
- (2) 水蒸気圧の計算-水蒸気圧の計算)
- 5 降水量と日射量に気を付けましょう
- 6 その他の気象庁GPVプロダクト
- 付録 分布図の地図投影
- おわりに
- 著作権について
はじめに¶
「GPV」とは聞きなれない言葉ですね。GPVは、緯度、経度、高度、時刻について整然と並ぶ格子点毎の値として提供されるデータを意味します。気象データといえば気象庁のアメダスがとても身近ですが、気象庁はGPV形式の気象データも多種配信しています。
例えば、メソ数値予報モデルGPVは、最長78時間先までの時別気象予測データで、海上を含む日本周辺域を約5km間隔で並ぶ格子点毎に予報値が収録されています。そして、収録される気象要素は、海面更正気圧、地上気圧、風ベクトル（U、V）、気温、相対湿度、降水量、雲量(全雲量、上層雲量、中層雲量、下層雲量)、日射量と多彩です。
気象予報値であること、特定地点での時間変化や、面的な広がりを把握できることから、様々なビジネスで活用できそうです。ただ、残念なことに 気象庁のGPVデータは、GRIB2 と呼ばれる特殊な形式で電子ファイル化されているなど、利用するには少々のスキルが必要とされます。
そこで、この研修では、気所庁のいくつかのGPVプロダクトを取り上げ、実際にデータをファイルから取り出して加工することを通し、気象庁GPVデータへの理解を深めます。
1 GPVデータ¶
1.1 GPVとは¶
GPVは、Grid Poin Data の略語とされ、緯度、経度、高度、時刻について整然と並ぶ格子点毎の値を意味します。そして、この形式で提供されるデータをGPVデータと呼びます。
時間的、空間的な広がりと連続性を持つ現象を近似的に表現する方法で、コンピュータ処理との親和性が高いので、気象予報モデルによる大気のシミュレーションや、気象レーダーによる降水観測などでこのデータ形式が用いられます。そして、気象庁は、これら予測や観測の結果を外部に提供する際にも、目的によりGPV形式を使用しています（ただし、シミュレーションやデータ処理に使用される格子配置と結果の配信で使われている格子配置とは同じとは限りません）。
なお、Grid Point Value データと類似の概念として メッシュ データがあります。英語で表現すると、Grid Square Value ないしは Grid Cell Value となります。これは、連続的に分布する事物を、縦の格子と横の格子で区切られる矩形の範囲で計量して一つの値とし、その集合として表現しようとする考え方のデータです。どちらで考えても大差ないことがほとんどです。例えば、気象庁の 推計気象分布（気温） などはGPVデータですが、格子点を中心とする矩形を想定しその代表的な気温として取り扱っても実用上全く問題は起きません。しかし、例えば、人口の分布など、メッシュデータとしては表現できるがGPVでは表現できない量もあることは、知っておいてください。
1.2 気象庁のGPVプロダクト¶
気象庁が対外的に作成する気象データは、気象庁ホームページ、気象庁情報カタログにおける気象に関する情報一覧で知ることができます。一覧から特定のカテゴリーをクリックするとそれに分類される気象情報の諸元が表の形で示されます。この表の、フォーマット等の解説 の列の、< 形式 > の下に「GRIB2」と記されていたら、それはGPV形式のデータです。
これから、大変多くの気象情報がGPVデータとして対外的に作成されていることがわかります。
「GRIB2」(グリブツー)は、GPVデータを格納しているファイルフォーマットの略号です。GRIBやその後継のGRIB2フォーマットは世界気象機構(WMO)が策定したもので、気象関係機関では広く使われていますがそこから一歩出ると殆ど使われていません。
気象庁は、気象情報プロダクトごとに、冊子「配信資料に関する仕様」（例）を作成し、それぞれのプロダクトの諸量が、GRIB2フォーマットにどのように載せられているかについて記載しています。ファイルの何ビット目から何ビット目にはには何がどのような数値で記録されているかといった情報が細かく書かれており、これに基づいてプログラムを書き下ろせばデータを取り出すことができます。しかし、このような行為は現代のコストパフォーマンスに全く適合しないので、この研修では、アメリカ大気海洋局(NOAA)の気候予測センター(CPC)が開発し公開するソフトウエア wgrib2 を使用してデータをデコードします。講師が検証した範囲では、これで実用上問題なくデコードできます。
1.3 気象庁のGPVプロダクトの入手¶
「アメダス気象データ分析チャレンジ！」で実習した通り、アメダスデータは、任意の地点、任意の期間のデータを気象庁ホームページからダウンロードできますが、GPVデータはそれができません。気象庁は、GPVプロダクトを気象業務支援センターを通じて配信しており、利用を希望する者は受信設備を用意したうえでセンターと契約を結び配信を受けます(有償です)。配信は最新のものを即時に受け取る方式のみで、過去のGPVデータを選択的に得ることはできません。
民間の気象会社のいくつかは、センターから受信したデータを再配布しており同様な配信を受けることが可能です。中には、アーカイブしている会社もあり、個別の交渉により、過去のGPVデータを有償で得ることができる場合もあります。
このほかとして、京都大学生存圏研究所では、生存圏データベースを運営しており、いくつかの気象庁GPVプロダクトをアーカイブしてダウンロードも可能にしていますが、利用できるのは、教育研究機関に限定されています。
すなわち、ビジネスにおける利用可能性を検討するために過去の一定期間のGPVデータを分析してみたい、というような用途に適したGPVデータは、きわめて入手しずらいのが現状です。
2 WXBCオリジナルPythonライブラリwxbcgribX¶
この講習では、読み込みを含めたGPVデータの処理に、Pythonライブラリ wxbcgribX を使用します。このライブラリは、WXBCの有志が開発したもので、GPVデータの処理工数を減らすことのできる関数を提供するものです。この講習の事前準備において、いくつかのPythonライブラリをインストールしましたが、 wxbcgribX は小規模なものなのでインストール作業(ファイルの配置とその場所をPythonに認識させること)は不要です。ライブラリの実体である、ファイル wxbcgribx.py をプログラムと同じフォルダ内に置き、プログラムにおいてインポートすれば利用できます。
インポートは、import 文を用いて以下のようにします。ライブラリをこのようにインポートすると、ライブラリが提供する関数を呼び出す際、「wx」を「wxbcgribx」の省略名として使えるようになります。
以下を実行してライブラリ wxbcgribX を利用する準備をしてください。
import wxbcgribx as wx
- ご注意
ライブラリ wxbcgribX は、内部でソフトウエア wgrib2 を使用します。 このため、このライブラリを使用するには、wgrib2 をインストールし、ファイル wxbcgribx.py の35行目に wgrib2 がインストールされているパスを指定する必要があります。
なお、wxbcgribX のライセンスは下記のとおりです。
Copyright (c) 2022 気象データ×IT勉強会
以下に定める条件に従い、本ソフトウェアおよび関連文書のファイル（以下「ソフトウェア」）の複製を取得するすべての人に対し、ソフトウェアを無制限に扱うことを無償で許可します。これには、ソフトウェアの複製を使用、複写、変更、結合、掲載、頒布、サブライセンス、および/または販売する権利、およびソフトウェアを提供する相手に同じことを許可する権利も無制限に含まれます。
上記の著作権表示および本許諾表示を、ソフトウェアのすべての複製または重要な部分に記載するものとします。
ソフトウェアは「現状のまま」で、明示であるか暗黙であるかを問わず、何らの保証もなく提供されます。ここでいう保証とは、商品性、特定の目的への適合性、および権利非侵害についての保証も含みますが、それに限定されるものではありません。 作者または著作権者は、契約行為、不法行為、またはそれ以外であろうと、ソフトウェアに起因または関連し、あるいはソフトウェアの使用またはその他の扱いによって生じる一切の請求、損害、その他の義務について何らの責任も負わないものとします。
3 GPVデータの読み込み¶
前置きはこれくらいにして、GRIBファイルからGPVデータを実際に取り出してみましょう。ここでは、気象庁 局地数値予報モデルGPV(通称LFM-GPV)のデータを例にとります。
3.1 局地数値予報モデルGPVについて¶
これから読むデータは、気象庁の正式名では「 LFM格子点データ（ファイル形式）」です。気象庁の局地数値予報モデルから作られるプロダクトの1つで、通常は、「局地数値予報モデルGPV」「LFM-GPV」などと呼ばれています。 このプロダクトに関する完全な説明は、「配信資料に関する仕様 No.12701」で得ることができます。
このプロダクトのデータは地上面の気象と、指定等圧面の気象とに大きく分けられます。ここでは、地上面気象データを読みます。地上気象の気象要素は、海面更正気圧、地上気圧、風ベクトル（U、V）、気温、相対湿度、降水量、雲量(全雲量、上層雲量、中層雲量、下層雲量)、日射量です。
GPVがカバーする緯度経度範囲は、北緯22.4～47.6度、東経120～150Eの範囲ですが、範囲の端に若干の無情報領域があります。地上気象の格子点間隔は南北方向0.02度、東西方向が0.025度です。
30分間隔、10時間先までの予報で、各時刻の予報値が1つのGRIBファイルに収められています(全21ファイル)。これが1時間毎に作成されています。
今回読むのは、2023年6月1日00Zを初期値とするものです。以下に、これを構成するGRIBファイルのファイル名を示します。このファイルは、フォルダChallenge3/jmadata/lfm に、年と月で整理されて保存されています。
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0000_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0030_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0100_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0130_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0200_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0230_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0300_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0330_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0400_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0430_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0500_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0530_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0600_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0630_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0700_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0730_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0800_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0830_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0900_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0930_grib2.bin
Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH1000_grib2.bin
ファイル名は、いくつかの文字列がアンダースコア「 _ 」で区切られた形となっています。中ほど、数字が14桁並んでいるのが目を引きます。これは予報の初期時刻を示し 中央協定時が 年(4桁)月(2桁)日(2桁)時(2桁)分(2桁)秒(2桁) で表現されています。その次、LFM は、局地数値予報モデルのプロダクトであることを示します。その２つ後ろ、Lsurf は、地上気象であることを示します。その次、FH と4桁の数字は、予報時間（初期時刻からどれくらい先の予報か）を示しています。F の後ろの文字が H なので、この場合は、数字4桁は 時(2桁)分(2桁) を示しています。
3.2 GRIBファイルからの読み込み¶
それではLFM-GPVのデータを読み込んでみましょう。読み込みにはwxbcgribXが提供する関数 getgpv を使用します。この関数には、オプションが様々ありますが、以下を基本として使用してください。
ds = wx.getgpv(path_to_file, element, timezone='Asia/Tokyo')
ここで、path_to_file は、読み込むGRIBファイルへのパス(path:フォルダ名とファイル名を繋げたもの、コンピュータにおけるデータの住所)、element は取り出す気象要素の記号です。結果は ds に入れられます。
LFM-GPVのGRIBデータは、21ファイルで一組です。このように、複数のGRIBファイルをまとめて読み込むには、ファイルへのパスをリストにして与えます。
気象要素の記号は、気象庁の数値予報モデルGPVプロダクトで共通しており下表のとおりです。複数の気象要素記号を、リストにして与えることも可能です。
| 記号 | 気象要素 |
|---|---|
| PRMSL | 海面更正気圧 |
| PRES | 地上気圧 |
| UGRD | 風ベクトル(U:南北方向) |
| VGRD | 風ベクトル(V:東西方向) |
| TMP | 気温 |
| RH | 相対湿度 |
| APCP | 降水量 |
| TCDC | 全雲量 |
| HCDC | 上層雲量 |
| MCDC | 中層雲量 |
| LCDC | 下層雲量 |
| DSWRF | 日射量 |
なお、関数 getgpv のオプションの全てとその用法について知りたいときには、Code Cellを一つ開いて、以下を実行すると表示させることができます。
help( wx.getgpv )
さて、日付や時刻の数字が入っているパスのリストを作るのは単純ながら結構面倒です。先々自動処理も考えると、初期時刻の文字列を入れるだけでファイルパスのリストが作れるようにしておいたほうがよいでしょう。そのようなとき、wxbcgribXが提供する関数 strft_range が便利です。この関数は、日時を示す様々な形式の文字列のリストを作り出すことができます。
このケースでのキモは、30分づつ増加する書式hhmmの文字列のリストです。この場合、関数は以下のように使います。実際にスクリプトを動かしてみましょう。
hhmm = wx.strft_range(start='0000', iter=21, format='%H%M', minutes=30)
hhmm が期待通りのリストになっているか、結果を確認してみましょう。以下を実行してください。
print(hhmm)
['0000', '0030', '0100', '0130', '0200', '0230', '0300', '0330', '0400', '0430', '0500', '0530', '0600', '0630', '0700', '0730', '0800', '0830', '0900', '0930', '1000']
期待通りの結果となったので、他の文字列も付け加えて、GRIBファイルへパスを作りましょう。
yyyymmdd = '20230601'
prdct = 'lfm'
grb_path = f'./jmadata/{prdct}/{yyyymmdd[0:4]}/{yyyymmdd[0:6]}'
grb_files = [f'{grb_path}/Z__C_RJTD_{yyyymmdd}000000_LFM_GPV_Rjp_Lsurf_FH{xxxx}_grib2.bin'
for xxxx in hhmm]
grb_files が期待通りのリストになっているか、結果を確認してみましょう。以下を実行してください。
grb_files
['./jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0000_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0030_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0100_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0130_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0200_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0230_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0300_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0330_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0400_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0430_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0500_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0530_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0600_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0630_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0700_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0730_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0800_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0830_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0900_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0930_grib2.bin', './jmadata/lfm/2023/202306/Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH1000_grib2.bin']
これでGRIBファイルへのリストができました。
ここで、文字列に関する、すごく基本的な事項を二つ復習します。
1つ目は、文字列からその一部を取り出す方法です。文字列からその一部を取り出すには、リストからその一部を取り出すのと同じ スライス という方法を用います。「20230601」という文字列が保存されている変数 yyyymmdd の、左の４文字だけを指定するには、yyyymmdd[0:4] とします。Python がリストや文字列を"スライス"する場所は以下のように決められているからです。
以下のCellに追記して、yyyymmddから、「日」にあたる「01」を取り出すスライスを作ってみてください。
yyyymmdd[6:8]
'01'
2つ目は、変数を埋め込んだ文字列の作り方です。Pythonでは、文字列はクォーテーションマーク「 " または ' 」で挟んで定義しますが、始まりを「 f" または f' 」とすると、文字列の中に変数値を埋め込むことができます。この時、埋め込む変数を中括弧で包みます。
例えば、変数 yyyymmdd に「20230601」という文字列が、xxxx に「0930」という文字列が記憶されている時、
「 Z__C_RJTD_20230601000000_LFM_GPV_Rjp_Lsurf_FH0930_grib2.bin 」
という文字列をつくるには、以下のようにします。
f'Z__C_RJTD_{yyyymmdd}000000_LFM_GPV_Rjp_Lsurf_FH{xxxx}_grib2.bin' は、文字列
この例では、変数を記憶する変数の値を埋め込みましたが、数値が記憶されている変数でもエラーなく値を埋め込むことができます。
以下のCellに手を加え、変数 message の内容が、「LFM-GPVのGRIBデータは、21ファイルで一組です。」となるようにしてください。
prdct = 'LFM-GPV'
num = 21
message = f'{prdct}のGRIBデータは、{num}ファイルで一組です。'
message
'LFM-GPVのGRIBデータは、21ファイルで一組です。'
このような文字列の作り方を、f-string と呼びます。
寄り道をしてしまいましたが、LFM-GPVのデータの読み込みに進みましょう。ここでは、気温と相対湿度を読むことにします。以下を実行してください。
elements = ['TMP', 'RH']
ds = wx.getgpv(grb_files, elements, timezone='Asia/Tokyo')
ファイルの読み出しには少し時間がかかります。Jupyterでは、コードが実行中の時、Code Cell の左側にある大括弧 [ ] の中の表示がアスタリスク「＊」となり、終了すると連番となります。表示が数字になるまで少し待ってください。
関数 getgpv が返すのは、変数 ds たった一つです。読み込んだのは、格子間隔2kmで日本を覆う10時間分の気象予報、しかも気温と湿度の2要素ですから、 ds はただの数値や配列ではないことは想像がつきます。では実際、どのようなものでしょうか。
Python では、変数名だけを"実行"すると、その変数の概要が表示されます。以下のセルを実行してください。
ds
<xarray.Dataset>
Dimensions: (latitude: 1261, longitude: 1201, time: 21)
Coordinates:
* latitude (latitude) float64 22.4 22.42 22.44 ... 47.58 47.6
* longitude (longitude) float64 120.0 120.0 120.1 ... 150.0 150.0
* time (time) datetime64[ns] 2023-06-01 ... 2023-06-01T10:0...
awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-...
Data variables:
TMP_1D5maboveground (time, latitude, longitude) float32 301.1 301.1 ... nan
RH_1D5maboveground (time, latitude, longitude) float32 82.76 82.73 ... nan
Attributes:
Conventions: COARDS
History: created by wgrib2
GRIB2_grid_template: 0ds の概要を示すフォームが表示されました。では、記載事項について確認してゆきましょう。
１行目、Dimensions の行には、GPVの次元が示されています。このGPVは、３次元であり、それぞれに、latitude、longitude、timeという名前が付けられていて、配列の要素数は、この順に、1261、1201、21であることがわかります。
その下、Coordinates(座標)のセクションには、各次元における座標(格子点の位置)がまとめられています。latitude の行の右端にあるお皿が重なったようなアイコンをクリックしてください。座標の緯度値が配列で示されます。これを見ると、GPVの南北範囲は北緯22.4度～47.6度となっており、3.1章で示した、LFM-GPVの概要ににおける値と一致していることが確認できます。
次元 time の座標についてもみておきましょう。"お皿アイコン"をクリックしてください。2023年06月01日00時から30分間隔で同日10時までのデータであることが分かります。アイコンをもう一度クリックして畳んでから、その左隣の"書類アイコン"をクリックしてください。reference_date の項目に「2023.06.01 00:00:00 UTC」とあり、この時刻がUTCであることが確認できます。
座標 time の下にawaretime というもう一つの座標があります。横の列に (time) とあるのは、この座標が次元 time に結びついていることを示します。この座標は、指定されたタイムゾーンにおける時刻です。GPVを読み込むとき、関数 getgpv のキーワード引数に「timezone='Asia/Tokyo'」を指定したので、データは日本標準時です。
その下、Data variables (データ値)のセクションには、今般読み込んだ気象データ本体の概要が示されています。左側のアイコンでデータやその正式名、単位などが確認できるのは座標と同じです。ここで、気象要素の記号が、GRIBファイルからGPVデータを取り出すときに使用したものと異なることに注意してください。
例えば、GRIBファイルから気温を取り出すときに使用した記号は「TMP」です。関数 getgpv は、データを読み込む際、この記号に、高度の情報を連結して「TMP_1D5maboveground」に変更します(Dは小数点)。これは、様々な高度におけるGPVが同じ記号となってしまうのを防ぐためです。
このように、GRIBファイルには、気象データの本体だけでなく、名称は単位、それが準拠する時空間格子など、気象データを分析する上で必要となる情報が収められており、関数 getgpv で作られる変数(この場合は ds )には、それらも一緒に格納されています。
4 GPVデータの処理¶
4.1 オブジェクト xarray.DataArray¶
上で見た通り、ds は、一つの変数の中に様々な数値や情報を抱えています。Pythonでは、このようなものを オブジェクト と呼びます。"データ形式"とでも捉えておけばよいでしょう。関数 getgpv が返したこのオブジェクトには正式な名称があり、xarray.Dataset といいます(以降、Datasetオブジェクトと略記)。Datasetオブジェクトは、座標が同じ複数のGPVデータを保持するコンテナのようなもの です。今回作成された ds の場合、気温と相対湿度のGPVデータを保持しています。
そして、当然ながら、Datasetオブジェクト からは、単一のGPVデータを取り出すことができます。Datasetオブジェクトから取り出した単一のGPVは、xarray.DataArray という別なオブジェクトになります (以降、DataArrayオブジェクトと略記)。
Dataset、DataArray の完全な説明は xarray Developers のホームページ にあります。
それでは、ds から気温のデータセットを取り出して、DataArray オブジェクトについての理解を深めましょう。
以下のようにすると、ds から気温のGPVデータを取り出すことができます。ここで、「TMP_1D5maboveground」は、 ds で使用されている気温の記号です。
da = ds['TMP_1D5maboveground']
ds のときと同様に、以下を実行して da の概要フォームをみてみましょう。DataArray の概要フォームでは、Dimensions(次元)のすぐ下にデータ本体の数値が表示されます。長くなることも多いので、邪魔であれば、左上にある「お皿アイコン」をクリックして畳んでください。
da
<xarray.DataArray 'TMP_1D5maboveground' (time: 21, latitude: 1261,
longitude: 1201)>
array([[[301.13507, 301.14288, 301.1507 , ..., nan, nan,
nan],
[301.11945, 301.12726, 301.12726, ..., nan, nan,
nan],
[301.096 , 301.10382, 301.10382, ..., nan, nan,
nan],
...,
[ nan, nan, nan, ..., nan, nan,
nan],
[ nan, nan, nan, ..., nan, nan,
nan],
[ nan, nan, nan, ..., nan, nan,
nan]],
[[301.1196 , 301.1274 , 301.13522, ..., nan, nan,
nan],
[301.09616, 301.1118 , 301.1118 , ..., nan, nan,
nan],
[301.08054, 301.08835, 301.09616, ..., nan, nan,
nan],
...
[ nan, nan, nan, ..., nan, nan,
nan],
[ nan, nan, nan, ..., nan, nan,
nan],
[ nan, nan, nan, ..., nan, nan,
nan]],
[[301.31738, 301.30957, 301.30176, ..., nan, nan,
nan],
[301.28613, 301.27832, 301.27832, ..., nan, nan,
nan],
[301.2627 , 301.25488, 301.24707, ..., nan, nan,
nan],
...,
[ nan, nan, nan, ..., nan, nan,
nan],
[ nan, nan, nan, ..., nan, nan,
nan],
[ nan, nan, nan, ..., nan, nan,
nan]]], dtype=float32)
Coordinates:
* latitude (latitude) float64 22.4 22.42 22.44 22.46 ... 47.56 47.58 47.6
* longitude (longitude) float64 120.0 120.0 120.1 120.1 ... 150.0 150.0 150.0
* time (time) datetime64[ns] 2023-06-01 ... 2023-06-01T10:00:00
awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-01T19:00:0...
Attributes:
short_name: TMP_1D5maboveground
long_name: Temperature
level: 1.5 m above ground
units: K4.2 GPV構成要素の取り出し¶
概要フォームからわかるように、 DataArray には、単一の気温データ本体だけでなく、名前や単位、準拠する緯度経度など、様々なデータが一纏まりになって保持されています。まずは、それら構成要素を個別に取り出す方法をおさえておきましょう。以下にその方法を列挙するので、実行して確認してみてください。
# 気象データの座標名を取り出す(タプル)
x = da.dims
print(x)
('time', 'latitude', 'longitude')
# 気象データの配列サイズを取り出す(タプル)
x = da.shape
print(x)
(21, 1261, 1201)
# 気象データ本体を取り出す(裸の配列(numpy.ndarrayオブジェクト)で)
x = da.data
print(x)
[[[301.13507 301.14288 301.1507 ... nan nan nan] [301.11945 301.12726 301.12726 ... nan nan nan] [301.096 301.10382 301.10382 ... nan nan nan] ... [ nan nan nan ... nan nan nan] [ nan nan nan ... nan nan nan] [ nan nan nan ... nan nan nan]] [[301.1196 301.1274 301.13522 ... nan nan nan] [301.09616 301.1118 301.1118 ... nan nan nan] [301.08054 301.08835 301.09616 ... nan nan nan] ... [ nan nan nan ... nan nan nan] [ nan nan nan ... nan nan nan] [ nan nan nan ... nan nan nan]] [[301.11456 301.1302 301.138 ... nan nan nan] [301.09894 301.10675 301.12238 ... nan nan nan] [301.0755 301.09113 301.09894 ... nan nan nan] ... [ nan nan nan ... nan nan nan] [ nan nan nan ... nan nan nan] [ nan nan nan ... nan nan nan]] ... [[301.14267 301.14267 301.13486 ... nan nan nan] [301.11923 301.11923 301.11923 ... nan nan nan] [301.0958 301.0958 301.0958 ... nan nan nan] ... [ nan nan nan ... nan nan nan] [ nan nan nan ... nan nan nan] [ nan nan nan ... nan nan nan]] [[301.22662 301.22662 301.2188 ... nan nan nan] [301.2032 301.2032 301.19537 ... nan nan nan] [301.17975 301.17194 301.17194 ... nan nan nan] ... [ nan nan nan ... nan nan nan] [ nan nan nan ... nan nan nan] [ nan nan nan ... nan nan nan]] [[301.31738 301.30957 301.30176 ... nan nan nan] [301.28613 301.27832 301.27832 ... nan nan nan] [301.2627 301.25488 301.24707 ... nan nan nan] ... [ nan nan nan ... nan nan nan] [ nan nan nan ... nan nan nan] [ nan nan nan ... nan nan nan]]]
# 緯度座標を取り出す(座標だけからなる Xarray.DataArray オブジェクトとして)
x = da['latitude']
x #概要フォームを表示させるため print文 を使わなかった
<xarray.DataArray 'latitude' (latitude: 1261)>
array([22.4 , 22.42, 22.44, ..., 47.56, 47.58, 47.6 ])
Coordinates:
* latitude (latitude) float64 22.4 22.42 22.44 22.46 ... 47.56 47.58 47.6
Attributes:
units: degrees_north
long_name: latitude# 緯度座標を取り出す(裸の配列(numpy.ndarrayオブジェクト)で)
x = da['latitude'].data
print(x)
[22.4 22.42 22.44 ... 47.56 47.58 47.6 ]
# 気象要素の正式名を取り出す(文字列)
x = da.attrs['long_name']
print(x)
Temperature
では、このGPVの単位を取り出すにはどのようにしたらよいでしょうか。下のセルを完成させて実行してみてください。
ヒント：単位の文字列が、概要フォームのどこに示されているかに注意してください
# 気象要素の単位を取り出す(文字列)
x = da.attrs['units']
print(x)
K
4.3 気象データへのアクセス¶
GPVデータの分析では、特定の時刻の分布を見たり、特定の場所での時間変化を調べたりなど、特定の部分を取り出すことが多くあります。また、その部分のデータを書き換えることもあります。xarray.DataArrayオブジェクトには、特定の場所のデータにアクセスする様々な方法が用意されています。
(1) Fancy indexing¶
最もシンプルに、各次元の要素番号を数字の並びで指定してデータを取り出すことができます。DataArrayオブジェクトは、多次元配列計算モジュール NumPy の ndarrayで配列要素を指定するのと同じ方法(Fancy indexing と呼ばれています)で要素を指定することができます。
Fancy indexing の方法で、緯度も0番目、経度も0番目、時刻は初期時刻を含めて10ステップ分のデータを取り出すことを考えてみましょう。DataArray オブジェクト da は、[時刻、経度、緯度] の形で配列化されています。そして、Fancy indexing では、0から9番目までの場合、「 0:10 」と記します。したがって、データにアクセスするには、da[0:10, 0, 0] と書きます。
下のCellを実行し、気温データが 10 要素の 1 次元のデータとなっていること、1 次元になってしまったあとでも、その格子の緯度経度や属性は付随していることを確認してください。
da[0:10,0,0]
<xarray.DataArray 'TMP_1D5maboveground' (time: 10)>
array([301.13507, 301.1196 , 301.11456, 301.09723, 301.08594, 301.08533,
301.06885, 301.0407 , 300.9884 , 300.9485 ], dtype=float32)
Coordinates:
latitude float64 22.4
longitude float64 120.0
* time (time) datetime64[ns] 2023-06-01 ... 2023-06-01T04:30:00
awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-01T13:30:0...
Attributes:
short_name: TMP_1D5maboveground
long_name: Temperature
level: 1.5 m above ground
units: K上の指定で取り出した気温変化を可視化してみましょう。オブジェクト DataArray は、GPVデータを簡単に可視化することができます。下のCellを実行してください。
da[0:10,0,0].plot(figsize=(4,3))
[<matplotlib.lines.Line2D at 0x2b9c37c1750>]
可視化に使用したのは、plot というメソッドです。メソッドとはオブジェクト専用の関数のことです。普通の関数は、関数を先に書いてその後ろ括弧を付け、対象をその中に入れますが、メソッドは、対象のオブジェクトの後ろにピリオドを置きメソッド名をその後ろ書きます。メソッドにオプションがあれば、その後ろに括弧を付けてオプションを与えます。
普通の関数：
関数名(対象のオブジェクト, 必要に応じオプション )
メソッド：
対象のオブジェクト.メソッド名( 必要に応じオプション )
対象となるオブジェクトは、オリジナルのものである必要はなく、同じオブジェクト形式を維持していればそれにも使用できます。上の例は、daから一部を切り出したもの(正確には一部しか見えないようにしたもの)ですが、これもDataArrayオブジェクトなので、メソッド plot が使用できます。
メソッド plot は、対象とするDataArrayオブジェクトが1次元の場合、自動的に折れ線グラフを返します。今回はグラフの大きさを指定するために、キーワード引数、「figsize=(4,3)」を与えています。
次の例として、最後の時刻、全ての緯度範囲、全ての経度範囲のデータを取り出して見ましょう。この場合はどのように指定したらよいでしょうか。
Fancy indexing では、「-1」は最後、「：」は全部を意味するので、[-1, :, : ] で指定できます。
では、この範囲で気温データを取り出して分布を描いてみましょう。メソッド plot は、2次元のオブジェクトに対し、自動的に分布図を作成します。分布図と一口に言っても表現方法は様々です。今回は、色合い、スケール、大きさをオプションで指定してみました。オプションは他にもたくさんあるので、興味のある方は、xarryのホームページで確認してください。
da[-1,:,:].plot(cmap='RdYlGn_r',vmin=280, vmax=305, figsize=(7,7))
<matplotlib.collections.QuadMesh at 0x2b9c39a3d60>
データアクセスからちょっと脱線しましょう。
メソッド plot は図のラベル類をうまく表示してくれますが、タイトルなど、カスタマイズしたい場合もあります。そのようなときは、描画ライブラリ MatPlotlib の力を借ります。
まずライブリをインポートしておきます。
import matplotlib.pyplot as plt
変更は、以下のスクリプトのように、まずメソッド plot で図を描き、そのあとで変更したい部分を上書きします。
実行して図の違いを確認してください。
da[-1,:,:].plot(cmap='RdYlGn_r',vmin=280, vmax=305, figsize=(7,7))
plt.title('New Title')
plt.ylabel('New y-label')
plt.draw()
(2) メソッド isel¶
Fancy indexing でデータを指定する場合、目的の要素に正しくアクセスするには、次元とその順序を予め知っている必要がありますが、メソッド isel を使うと、次元の名前で要素を指定できます。次元の並び順を気にする必要もありません。これはスクリプトの可読性の向上に役立ちます。
isel を使う場合、範囲は関数 slice で指定します。また、記載しない次元については、その次元の要素すべてを選択したことになります。
なお、iselはメソッド(=関数)なので、引数は大括弧「[ ]」ではなく小括弧「( )」で括ります。
下のCellを実行し、その結果が Fancy Indexig での結果と同じであることを確認してください。
da.isel(time=slice(0,10), latitude=0, longitude=0)
<xarray.DataArray 'TMP_1D5maboveground' (time: 10)>
array([301.13507, 301.1196 , 301.11456, 301.09723, 301.08594, 301.08533,
301.06885, 301.0407 , 300.9884 , 300.9485 ], dtype=float32)
Coordinates:
latitude float64 22.4
longitude float64 120.0
* time (time) datetime64[ns] 2023-06-01 ... 2023-06-01T04:30:00
awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-01T13:30:0...
Attributes:
short_name: TMP_1D5maboveground
long_name: Temperature
level: 1.5 m above ground
units: K(3) メソッド sel¶
メソッド sel を用いると、要素を、要素番号ではなく座標値で選び出すことができます。メソッド sel に続く括弧の中に、時刻、緯度、経度を指定するキーワード引数を記入します。緯度、経度の指定では、指定する格子の緯度経度値を与えます。範囲を指定する場合は関数 slice を使用します。
da には日本標準時の座標 awaretime が付いているので、この値で時刻を指定したいところですが、メソッド sel に使用できる座標は、次元として使われているもの(time, latitude, longitude)だけで、残念ながら、awaretime は使用できません。
これを利用して、2023年6月1日03時UTCにおける中部地方の一部(１次メッシュ5135～5337の範囲)の予測気温分布図を描いてみましょう。下に示すのが、図を描かせるスクリプトです。
緯度経度範囲は、最大値と最小値を数値で一つ一つ書きこんでもよいのですが、応用が利くように、別途リストで用意し、それを与えるようにしています。また、Pythonは、括弧の中ではどのように改行してもよいので、見やすいように、時刻の指定、緯度範囲の指定、経度範囲の指定を1行ずつにしました。 これを実行してください。
area = [34, 36, 135, 138] #１次メッシュ5135～5337
Ta = da.sel(time='2023-06-01T03:00',
latitude=slice(area[0],area[1]),
longitude=slice(area[2],area[3]) )
Ta.plot(cmap='RdYlGn_r')
plt.title(Ta.awaretime.data) # 図のタイトルを日本標準時に
plt.draw()
次に、この図の中心近くに位置する名古屋地方気象台付近の気温の推移を折れ線グラフにしてみましょう。気象台の位置は北緯35.166667度、東経136.965度です。この緯度経度の値を取り出せばよいのですが、残念ながら(というか普通は)GPVの格子点にこれと一致するものはありません。そのようなとき、メソッド sel にキーワード引数「 method='nearest' 」を加えることで、指定した緯度経度に最も近い格子点を選ぶことができます。
lat, lon = 35.166667, 136.965 # 名古屋地方気象台
Ta = da.sel(latitude=lat, longitude=lon, method='nearest')
Ta.plot(x='awaretime',figsize=(10,3))
plt.title(f'(N{Ta.latitude.data:.4f}, E{Ta.longitude.data:.4f})')
plt.draw()
3.2 の'横道'で説明した f-string を上のスクリプトの5行目でも使用しています。浮動小数の値を文字列に埋め込むとき、この例のようにすると、表示桁数を指定することができます。
4.4 気象データの演算¶
DataArrayオブジェクトのデータに対する演算は、Numpy の配列計算と同じ作法で行うことができます。オブジェクト同士やオブジェクトとスカラー（単一の数値）の演算が可能です。
(1) 絶対温度からセルシウス度への変更¶
DataArrayオブジェクトとスカラーの演算の一例として、絶対温度(K)で記録されている気温データから、273.15を一律に引いて摂氏(℃)に変換し、上と同じグラフを描いてみましょう。
# 絶対温度をセルシウス度への変更
Tc = Ta - 273.15
Tc.plot(x='awaretime',figsize=(10,3))
plt.title(f'(N{Tc.latitude:.4f}, E{Tc.longitude:.4f})')
plt.draw()
値自体は、正しく摂氏に変更されていますが、縦軸のタイトルが変わってしまいました。概要フォームを表示させて何が起きているのか確認してみましょう。
Tc
<xarray.DataArray 'TMP_1D5maboveground' (time: 21)>
array([22.703827, 23.563354, 23.941132, 24.111298, 24.412506, 24.42752 ,
24.395416, 24.359467, 24.354034, 24.173492, 23.778748, 23.418915,
22.938202, 22.602905, 22.321442, 21.64264 , 21.30368 , 21.232788,
21.187988, 21.178192, 19.972076], dtype=float32)
Coordinates:
latitude float64 35.16
longitude float64 137.0
* time (time) datetime64[ns] 2023-06-01 ... 2023-06-01T10:00:00
awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-01T19:00:0... Attributes(付随情報)が削除されています。
DataArray オブジェクトは、データ本体に変更が加えられると付随情報を引き継ぎません。演算して異なるものを作ったのですから属性が保存されることは逆に問題であり、これを避けるために引き継がれないのです。演算結果の付随情報は自ら書き込む必要があります。
付随情報は、DataArray オブジェクトの内部にあらかじめ用意されている変数 attrs に記録する決まりになっているので、辞書の形で代入します。
Tc.attrs = {"long_name":"Air temperature near Nagoya Loc. Met. Obs. ", "units":"Celsius"} # 付随情報を２つ設定
今度は意図通りのグラフが描けるはずです。 下のCellを実行し確認してみましょう。
Tc.plot(x='awaretime', figsize=(10,3))
plt.title(f'(N{Tc.latitude:.4f}, E{Tc.longitude:.4f})')
plt.draw()
(2) 水蒸気圧の計算¶
次は、２つ気象要素を使用する少し複雑な演算をしてみましょう。気象庁の数値予報モデルGPVは、湿度を相対湿度で提供しますが、湿度を表現する指標はこれ以外にもいくつかあります。その一つである水蒸気圧は、着目する気塊から結露や降水で水蒸気が除かれてしまわない限り一定に保たれるという特徴をもっています。相対湿度 H [%]、気温 T [℃] の空気の水蒸気圧 e [hPa] は、近似的に下の式で表現できます。
$$ e = 6.1078 \times 10^{\left( \frac{7.5T}{T+237.3} \right)} \times \frac{H}{100} $$
気温と湿度のGPVデータを演算して、水蒸気圧を計算してみましょう。この式は下のようにコーディングすることができます。このCellを実行してください。
#気象要素の取り出し
T = ds['TMP_1D5maboveground'] - 273.15
H = ds['RH_1D5maboveground']
#水蒸気圧の計算
e = 6.1078 * 10**((7.5*T)/(T+237.3)) * H/100
#付随情報の設定
e.attrs = {'long_name':'water vapor pressure', 'units':'hPa'}
準拠する座標が同じDataArrayオブジェクトは、ndarray 配列と同じ記法で計算でき、計算結果はDataArrayオブジェクトになります。もちろん、Fancy indexing で一部だけにしてもDataArrayオブジェクトです。従って、メソッド plot を適用すれば、分布図を描くことができます。
#分布の可視化
e[-1,:,:].plot(cmap='RdYlGn_r')
plt.title(e[-1].awaretime.data) # 図のタイトルを日本標準時に
plt.draw
<function matplotlib.pyplot.draw()>
・注意してください¶
Pythonでは、b = a のようなスクリプトがあった場合、通常は、変数 b が新規に作られ、a の値が代入されますが、DataAttayオブジェクトの場合、新規には作られません。「 b という変数名でも a が保持するデータにアクセスできるようになる」という処理がされます。このため、別物と思ってb の内容を変更してしまうと、おおもとの a の内容まで変更されてしまいます。
例えば、先に学習した、温度の単位の kelvin から Celsius度 への変更において、以下のようにしてしまうと、Ta と Tc の両方が Celsius度になってしまいます。
Tc = Ta
Tc = Ta - 273.15
あるDataAttayオブジェクトの内容に手を加えた別なオブジェクトを作りたいときは、以下のようにメソッド copy を使ってオブジェクトの複製をまず作り、それに対し変更を加えます。
b = a.copy()
5 降水量と日射量に気を付けましょう¶
気温や湿度などの気象量は、特定の格子点(時間で言えば特定の時刻)における値を定義できます。降水や日射についても、それぞれ、降水強度、日射放射強度により、格子点値を定義できますが、これらの量は、時間・空間的変化が大きいことから、瞬時値よりもそれを時間積算した値のほうがより実用価値があります。このため、気象庁のGPVプロダクトでは、 降水量と日射量は積算値または平均値で提供されています。これは、第1章で説明したメッシュ の考え方に他なりません。
すなわち、注意すべき第1点目として、降水量と日射は、空間については格子点値だが、時間においては"メッシュ値"であるということがあげられます。
注意すべき第2点目は、気象庁GPVでは、降水量と日射量に関し、「どのような積算か？」がプロダクトによって異なる ことです。
LFM-GPVでは、降水量については初期時刻からの積算値 が記録され、日射量については 前予報時刻(30分前)からの平均値 が記録されていますが、GPVプロダクトによってそれはまちまちです。
注意すべき第3点目は、これらが積算や平均であることから、初期時刻の値が存在しないことです。これを反映して、各予報時間ステップが独立したGRIBファイルになっていLFM-GPV では、初期時刻のGRIBファイルには、律儀にも降水量と日射にかかわる情報が全く入っていません。
これが何を引き起こすかというと、他の要素と同じつもりで降水量や日射量を読み込もうとすると、以下に示すように、関数 getgpv は最初のGRIBファイルで躓き、エラーを起こします。したがって、降水量と日射量をGRIBファイルから取り出すときは、初期値のファイルはリストに含めないでください。
GRIBファイルのリストから初期値のファイルだけを除くには、以下２行目にあるようにスライスを使います。
elements = ['APCP', 'DSWRF']
ds5 = wx.getgpv( grb_files[1:] , elements, timezone='Asia/Tokyo')
#スライスにより、初期時刻のGRIBファイル名を除外する
それでは、2023年6月1日に雨が降ったアメダス「那覇」に最も近い格子点における降水量の時間変化をグラフにしてみましょう。以下を実行してください。表示されたグラフの、縦軸の名称と単位にも注意してください。
lat, lon = 26.20666667, 127.6866667 #「那覇」の緯度経度
Pr = ds5['APCP_surface'].sel(latitude=lat,longitude=lon, method='nearest')
Pr.plot(x='awaretime', figsize=(10,3))
plt.title(f'(N{Pr.latitude.data:.4f}, E{Pr.longitude.data:.4f})')
plt.draw()
確かに、LFM-GPVの降水量は初期時刻からの積算値です。時間雨量を出すためには、差分を計算する必要があります。
次は、日射量の時間変化をグラフにしてみましょう。那覇とは対照的によく晴れた、アメダス「札幌」付近の格子点について可視化します。
lat, lon = 43.06, 141.3283333 #「札幌」の緯度経度
Sr = ds5['DSWRF_surface'].sel(latitude=lat, longitude=lon, method='nearest')
Sr.plot(x='awaretime', figsize=(10,3))
plt.title(f'(N{Sr.latitude.data:.4f}, E{Sr.longitude.data:.4f})')
plt.draw
<function matplotlib.pyplot.draw()>
6 その他の気象庁GPVプロダクト¶
これまで、GPVデータの基本的な取り扱い方を、LFM-GPV を素材として学びましたが、1.2 章で見た通り、気象庁は多種多様なGPVデータを対外的に作成しています。そこで、これ以外のいくつかのGPVプロダクトについて簡単な処理をしながら概観してみましょう。
この章では、これまでの実習で使用したPythonカーネルをクリアして、改めて開始することにします。
ツールバーにある輪のようなアイコンをクリックし、さらに、確認の [Restart] ボタンをクリックして今実行中のPythonカーネルをりセットし、下記のインポート文を実行してください。
# ライブラリ のインポート
import wxbcgribx as wx
import matplotlib.pyplot as plt
6.1 MSM格子点データ（ファイル形式）¶
このプロダクトは、メソ数値予報モデルから作られるGPVの一つです。通常は、「メソ数値予報モデルGPV」「MSM-GPV」などと呼ばれています。主に明日、明後日の天気予報を行うのに用いられます。 このプロダクトに関する完全な説明は、 配信資料に関する仕様No.12601 にあります。
このプロダクトのデータも地上面の気象と、指定等圧面の気象とに大きく分けられます。地上気象の気象要素は、海面更正気圧、地上気圧、風ベクトル（U、V）、気温、相対湿度、降水量、雲量(全雲量、上層雲量、中層雲量、下層雲量)、日射量です。なお、MSM-GPVの降水量は直前1時間の積算値 で記録されています。また、日射量は前予報時間からの平均値です。
GPVがカバーする緯度経度範囲は、北緯22.4～47.6度、東経120～150Eの範囲です。地上気象の格子点間隔は南北方向0.05度、東西方向が0.0625度です。
毎日03,06,09,15,18,21UTCを初期時刻として39時間先までの予報が1時間ステップで作られています。また、00,12UTCを初期時刻として78時間先までの予報が作られています。予報は3個または5個のファイルに分割されて配信されます。
ここで読むのは、2023年6月1日00時UTCを初期値とする地上気象です。この時刻の予報は、78時間先までで、下に示すような5個のファイルで構成されます。サンプルデータのファイルは、フォルダ Challenge3/jmadata/msm に、年と月で整理されて保存されています。
Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin
Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH16-33_grib2.bin
Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH34-39_grib2.bin
Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH40-51_grib2.bin
Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH52-78_grib2.bin
ファイル名は、いくつかの文字列がアンダースコア「 _ 」で区切られた形となっています。中ほど、14桁の数字の並びは予報の初期時刻で 年(4桁)月(2桁)日(2桁)時(2桁)分(2桁)秒(2桁) です。中央協定時で表現されています。その次、MSM は、メソ数値予報モデルのプロダクトであることを示します。その２つ後ろ、Lsurf は、地上気象であることを示します。
その次、FH に続くハイフンを挟む2桁の数字は、それぞれファイルが収録する予報の初めの時刻と最後の時刻を示しています。F の後ろの文字が H なので、2桁の数字は、時(2桁)と解釈します。つまり、一番上のファイルの場合、初期時刻から15時間先予報までのデータが収録されていることを示しています。
上記GRIBファイルへのパスのリストは以下のようにすると作ることができます。
yyyymmdd = '20230601'
prdct = 'msm'
fh = ['00-15', '16-33', '34-39', '40-51', '52-78']
grb_path = f'./jmadata/{prdct}/{yyyymmdd[0:4]}/{yyyymmdd[0:6]}'
grb_files = [f'{grb_path}/Z__C_RJTD_{yyyymmdd}000000_MSM_GPV_Rjp_Lsurf_FH{xx}_grib2.bin'
for xx in fh]
# 結果を確認したい人はこのセルを実行して内容を表示させてください
grb_files
['./jmadata/msm/2023/202306/Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin', './jmadata/msm/2023/202306/Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH16-33_grib2.bin', './jmadata/msm/2023/202306/Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH34-39_grib2.bin', './jmadata/msm/2023/202306/Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH40-51_grib2.bin', './jmadata/msm/2023/202306/Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH52-78_grib2.bin']
GRIBファイルに保持されている気象要素は、ライブラリ wxbcgribx が提供する関数 getvarname で確認することができます。それには以下のようにします。
print( wx.getvarname(grb_files[0]) )
['APCP', 'DSWRF', 'HCDC', 'LCDC', 'MCDC', 'PRES', 'PRMSL', 'RH', 'TCDC', 'TMP', 'UGRD', 'VGRD']
MSM-GPVのGRIBファイルは、複数の時刻の予報が一つのファイルにまとめられているので、LFM-GPV のように、GRIBファイルによって、収録要素が異なるということはありません。すなわち、関数 getvarname の引数に grb_files[0] を与えても grb_files[1] を与えても結果は同じです。
それでは、気温、相対湿度、降水量を読んでみましょう。以下を実行してください。概要フォームが表示されるので、読み込まれたGPVの座標や気象要素を確認してみてください。
elements = ['TMP','RH','APCP']
ds = wx.getgpv(grb_files,elements, timezone='Asia/Tokyo')
ds
<xarray.Dataset>
Dimensions: (latitude: 505, longitude: 481, time: 79)
Coordinates:
* latitude (latitude) float64 22.4 22.45 22.5 ... 47.5 47.55 47.6
* longitude (longitude) float64 120.0 120.1 120.1 ... 149.9 150.0
* time (time) datetime64[ns] 2023-06-01 ... 2023-06-04T06:0...
awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-...
Data variables:
TMP_1D5maboveground (time, latitude, longitude) float32 301.2 ... 277.9
RH_1D5maboveground (time, latitude, longitude) float32 83.08 ... 90.82
APCP_surface (time, latitude, longitude) float32 nan nan ... 0.0 0.0
Attributes:
Conventions: COARDS
History: created by wgrib2
GRIB2_grid_template: 0セクション Data variables、APCP_surfac の項目の"お皿"アイコンをクリックして、降水量のデータを表示させてください。データの初めの部分が、値ではなく nan と表示されていることが分かります。これは初期時刻の値がないことを反映しています。
それでは、GPVデータを可視化しましょう。
まず、空間範囲全体の分布図を見ることにします。以下を実行して、2023年6月1日3:00UTCにおける相対湿度の分布図を表示させてください。台湾の東方に、台風3号が予測されているのがわかります。
H_msm = ds['RH_1D5maboveground']
H_msm_23060103 = H_msm.sel(time='2023-06-01T03:00')
H_msm_23060103.plot(cmap='RdYlGn_r', vmin=30,vmax=100)
plt.title(H_msm_23060103.awaretime.data)
plt.draw
<function matplotlib.pyplot.draw()>
次に、4.3 (3) において、LFM-GPVで作図したのと同じ条件で気温データを可視化し、両者を比較してみましょう。
まず、2023年06月01日03時(UTC)における中部地方の気温分布図を作成します。
以下を実行して下さい。
area = [34, 36, 135, 138] #３次メッシュ5135～5337
T_msm = ds['TMP_1D5maboveground']
T_msm_23060103 = T_msm.sel(time='2023-06-01T03:00',
latitude=slice(area[0],area[1]),
longitude=slice(area[2],area[3]) )
T_msm_23060103.plot(cmap='RdYlGn_r')
plt.title(T_msm_23060103.awaretime.data)
plt.draw
<function matplotlib.pyplot.draw()>
名古屋地方気象台付近の気温推移についても可視化しましょう。
lat, lon = 35.166667, 136.965 #名古屋地方気象台
T_msm_nag = T_msm.sel(latitude=lat, longitude=lon, method='nearest')
T_msm_nag.plot(x='awaretime',figsize=(10,3))
plt.title(f'(N{T_msm_nag.latitude:.4f}, E{T_msm_nag.longitude:.4f})')
plt.draw
<function matplotlib.pyplot.draw()>
数値予報モデルの格子点間隔と予報期間はトレードオフの関係にあります。LFM-GPVは、格子点間隔がおよそ2kmと空間表現性には優れていますが、10時間先までの予報しかありません。一方、MSM-GPVは、格子点間隔は5kmと劣りますが、予報期間は78時間と10倍も長いです。
以下に、同じ地点での相対湿度と降水量との折れ線グラフを作成するスクリプトを示すので、実行してください。
#Hr_msm = ds_msm['RH_1D5maboveground'] # 実行済
H_msm_nag = H_msm.sel(latitude=lat,longitude=lon, method='nearest')
H_msm_nag.plot(x='awaretime',figsize=(10,3))
plt.title(f'(N{H_msm_nag.latitude:.4f}, E{H_msm_nag.longitude:.4f})')
plt.draw
<function matplotlib.pyplot.draw()>
P_msm = ds['APCP_surface']
P_msm_nag = P_msm.sel(latitude=lat,longitude=lon, method='nearest')
P_msm_nag.plot(x='awaretime',figsize=(10,3))
plt.title(f'(N{P_msm_nag.latitude:.4f}, E{P_msm_nag.longitude:.4f})')
plt.draw
<function matplotlib.pyplot.draw()>
MSMの降水量が、LFMのような初期値からの積算にではないことに、改めて注意してください。
6.2 全球数値予報モデルGPV（日本域）¶
このプロダクトは、全球数値予報モデルから作られるGPVの一つです。このモデルは、文字通り全球を計算対象としていて、主に週間天気予報に用いられます。 このプロダクトとは別に、「全球数値予報モデルGPV（全球域）」というものも存在します。もとは同じモデルですが、GPVプロダクトの格子点間隔は同じではありません。このプロダクトに関する説明は、 配信資料に関する仕様 No.12501 です。
このプロダクトのデータも地上面の気象と、指定等圧面の気象とに大きく分けられます。地上気象の気象要素は、海面更正気圧、風ベクトル（U、V）、気温、相対湿度、積算降水量、雲量（全雲量、上層雲量、中層雲量、下層雲量）、日射量です。なお、GSM-GPVの降水量は初期時刻からの積算値 で記録されています。また、日射量は前予報時間からの平均値です。
GPVがカバーする緯度経度範囲は、北緯20～50度、東経120～150度の範囲です。地上気象の格子点間隔は南北方向0.1度、東西方向が0.125度です。
毎日00,06,12,18UTCを初期時刻として132時間先までの予報が1時間ステップで作られています。00,12UTCを初期時刻とする予報につては、135～264時間先までのものが作られていますが、予報の時間間隔が3時間です。
ここで読むのは、2023年6月1日00Zを初期値とする地上気象です。以下に、これを構成するGRIBファイルのファイル名を示します。サンプルファイルは、フォルダ Challenge3/jmadata/gsm に、年と月で整理されて保存されています。
Z__C_RJTD_20230601000000_GSM_GPV_Rjp_Gll0p1deg_Lsurf_FD0000-0100_grib2.bin
Z__C_RJTD_20230601000000_GSM_GPV_Rjp_Gll0p1deg_Lsurf_FD0101-0200_grib2.bin
Z__C_RJTD_20230601000000_GSM_GPV_Rjp_Gll0p1deg_Lsurf_FD0201-0300_grib2.bin
Z__C_RJTD_20230601000000_GSM_GPV_Rjp_Gll0p1deg_Lsurf_FD0301-0400_grib2.bin
Z__C_RJTD_20230601000000_GSM_GPV_Rjp_Gll0p1deg_Lsurf_FD0401-0500_grib2.bin
Z__C_RJTD_20230601000000_GSM_GPV_Rjp_Gll0p1deg_Lsurf_FD0501-0512_grib2.bin
Z__C_RJTD_20230601000000_GSM_GPV_Rjp_Gll0p1deg_Lsurf_FD0515-0700_grib2.bin
Z__C_RJTD_20230601000000_GSM_GPV_Rjp_Gll0p1deg_Lsurf_FD0703-0900_grib2.bin
Z__C_RJTD_20230601000000_GSM_GPV_Rjp_Gll0p1deg_Lsurf_FD0903-1100_grib2.bin
ファイル名は、いくつかの文字列がアンダースコア「 _ 」で区切られた形となっています。中ほど、14桁の数字の並びは予報の初期時刻で 年(4桁)月(2桁)日(2桁)時(2桁)分(2桁)秒(2桁) です。中央協定時で表現されています。その次、GSM は、全球数値予報モデルのプロダクトであることを示します。その4つ後ろ、Lsurf は、地上気象であることを示します。その次、FD に続くハイフンを挟む4桁の数字は、それぞれファイルが収録する予報の最初と最後の時刻を示しています。F の後ろの文字が D なので、4桁の数字は、日(2桁)時(2桁)と解釈します。つまり、一番上のファイルの場合、初期時刻から1日先予報までのデータが収録されていることを示しています。
上記GRIBファイルへのパスのリストは以下のようにすると作ることができます。
yyyymmdd = '20230601'
prdct = 'gsm'
fd = ['0000-0100','0101-0200','0201-0300','0301-0400','0401-0500','0501-0512',
'0515-0700','0703-0900','0903-1100']
grb_path = f'./jmadata/{prdct}/{yyyymmdd[0:4]}/{yyyymmdd[0:6]}'
grb_files = [f'{grb_path}/Z__C_RJTD_{yyyymmdd}000000_GSM_GPV_Rjp_Gll0p1deg_Lsurf_FD{xx}_grib2.bin'
for xx in fd]
では、気温、相対湿度、降水量を読みましょう。以下を実行してください。概要フォームが表示されるので、読み込まれたGPVの座標や気象要素を確認してみてください。
elements = ['TMP','RH','APCP']
ds = wx.getgpv(grb_files, elements, timezone='Asia/Tokyo')
ds
<xarray.Dataset>
Dimensions: (latitude: 301, longitude: 241, time: 177)
Coordinates:
* latitude (latitude) float64 20.0 20.1 20.2 20.3 ... 49.8 49.9 50.0
* longitude (longitude) float64 120.0 120.1 120.2 ... 149.9 150.0
* time (time) datetime64[ns] 2023-06-01 ... 2023-06-12
awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-12...
Data variables:
TMP_2maboveground (time, latitude, longitude) float32 301.5 301.5 ... 278.5
RH_2maboveground (time, latitude, longitude) float32 83.15 83.59 ... 96.72
APCP_surface (time, latitude, longitude) float32 nan nan ... 33.56
Attributes:
Conventions: COARDS
History: created by wgrib2
GRIB2_grid_template: 0気象要素の補足情報から、GMSでは、地表面に最も近い層の高さが 2 m に設定されていることがわかります。
それでは、MSM-GPVと同じ設定の図を作ってゆきましょう。最初は6月1日03時(UTC)における全範囲の湿度場です。
H_gsm = ds['RH_2maboveground']
H_gsm_23010103 = H_gsm.sel(time='2023-06-01T03:00')
H_gsm_23010103.plot(cmap='RdYlGn_r')
<matplotlib.collections.QuadMesh at 0x2b9c6268160>
次は、中部域の気温分布図です。
area = [34, 36, 135, 138] #１次メッシュ5135～5337
T_gsm = ds['TMP_2maboveground']
T_gsm_23060103 = T_gsm.sel(time='2023-06-01T03:00',
latitude=slice(area[0],area[1]),
longitude=slice(area[2],area[3]) )
T_gsm_23060103.plot(cmap='RdYlGn_r')
plt.title(T_gsm_23060103.awaretime.data)
plt.draw()
格子点間隔が 10 km もなると、心の目を使わないと伊勢湾や若狭湾を見つけ出すのは難しいですね。
それはともかく、名古屋地方気象台付近の気象の推移を可視化しましょう。
まず気温です。
lat, lon = 35.166667, 136.965 #名古屋地方気象台
T_gsm_nag = T_gsm.sel(latitude=lat, longitude=lon, method='nearest')
T_gsm_nag.plot(x='awaretime',figsize=(10,3))
plt.title(f'(N{T_gsm_nag.latitude:.4f}, E{T_gsm_nag.longitude:.4f})')
plt.draw
<function matplotlib.pyplot.draw()>
相対湿度です。
H_gsm = ds['RH_2maboveground']
H_gsm_nag = H_gsm.sel(latitude=lat, longitude=lon, method='nearest')
H_gsm_nag.plot(x='awaretime',figsize=(10,3))
plt.title(f'(N{H_gsm_nag.latitude:.4f}, E{H_gsm_nag.longitude:.4f})')
plt.draw
<function matplotlib.pyplot.draw()>
降水量です。
P_gsm = ds['APCP_surface']
P_gsm_nag = P_gsm.sel(latitude=lat, longitude=lon, method='nearest')
P_gsm_nag.plot(x='awaretime',figsize=(10,3))
plt.title(f'(N{P_gsm_nag.latitude:.4f}, E{P_gsm_nag.longitude:.4f})')
plt.draw
<function matplotlib.pyplot.draw()>
GSM-GPVのデータは、132時間先までは1時間ステップですがその先は3時間ステップになるので、今回読み込んだGPVの場合、6月6日12UTC (6月6日21JST) までは1時間ステップで、その次のデータは同日15UTCで3時間ステップです。気温や湿度の折れ線を注意深く見てください。これに対応し、折れ線の細かさがこの前後で変わっていることが分かると思います。
6.3 推計気象分布¶
このプロダクトは、地上観測や地球観測衛星のデータから作られるGPVです。観測結果をGPVにしたものなので、未来のデータはありません。気象要素は気温、日照時間、天気の3種類です。このプロダクトに関する説明資料は、気温については 配信資料に関する仕様 No.13801 、天気については 配信資料に関する仕様 No.13802 、日照時間については 配信資料に関する仕様 No.13803 です。
GPVがカバーする緯度経度範囲は北緯20度～48度、東経118度～150度で、格子点間隔は南北方向0.008333度、東西方向が0.0125度です。この格子は、国土数値情報における標準地域メッシュ(3次メッシュ)の中心点と一致しています。そして、値は3次メッシュの代表的な値とされています。つまり、厳密な意味ではGPVではなくメッシュデータです。
毎正時に、その時刻のGPVが作られていて、約20分後に気象業務支援センターから配信されています。GRIBファイル名は、以下のようなものです。サンプルデータはフォルダ Challenge3/jmadata/obs_gpv に、年と月で整理されて保存されています。
Z__C_RJTD_20230601000000_OBS_GPV_Rjp_Ggis1km_Ptt_A202306010000_grib2.bin
Z__C_RJTD_20230601000000_OBS_GPV_Rjp_Ggis1km_Pds60_A202306010000_grib2.bin
ファイル名は、いくつかの文字列がアンダースコア「 _ 」で区切られた形となっています。14桁の数字の並びが2か所現れます。これらは、GPVの時刻で 年(4桁)月(2桁)日(2桁)時(2桁)分(2桁)秒(2桁) です。中央協定時で表現されています。中ほどの OBS_GPV はプロダクトが推計気象分布であることを示します。後ろの14桁の数字の手前、Ptt や Pds60 は気象要素を示し、それぞれ、気温と日照時間を表しています。
では上記のうち、日照時間のGPVデータひとつを読んでみましょう。まずGRIBファイルで使用されている気象要素の記号を把握します。
# GRIBファイルへのパスを作る
yyyymmdd = '20230601'
prdct = 'obs_gpv'
param = 'Pds60'
grb_path = f'./jmadata/{prdct}/{yyyymmdd[0:4]}/{yyyymmdd[0:6]}'
grb_file = f'{grb_path}/Z__C_RJTD_{yyyymmdd}000000_OBS_GPV_Rjp_Ggis1km_{param}_A{yyyymmdd}0000_grib2.bin'
# 収録される気象要素の記号を取得する
print(wx.getvarname(grb_file))
['SUNSD', 'var discipline=0 center=34 local_table=0 parmcat=6 parm=194']
このプロダクトには、２種類のGPVデータが保持されていることが分かります(二つ目は大変奇妙な名称になっていますが)。これらを読み出してみましょう。
elements = ['SUNSD', 'var discipline=0 center=34 local_table=0 parmcat=6 parm=194']
ds = wx.getgpv(grb_file,elements)
ds
<xarray.Dataset>
Dimensions: (latitude: 3360, longitude: 2560, time: 1)
Coordinates:
* latitude (latitude) float64 20.0 20.01 20.02 ... 47.98 47.99 48.0
* longitude (longitude) float64 118.0 118.0 118.0 ... 150.0 150.0
* time (time) datetime64[ns] 2023-06-01
Data variables:
SUNSD_surface (time, latitude, longitude) float32 nan nan ... nan nan
var0_6_194_surface (time, latitude, longitude) float32 nan nan ... nan nan
Attributes:
Conventions: COARDS
History: created by wgrib2
GRIB2_grid_template: 0 概要フォームから、時間方向の座標は1つだけながら、データ本体の次元は3であることがわかります。そして、格子点数は、実に860万点超と膨大であることも分かります。また、DataSetオブジェクトにおける記号「SUNSD_surface」が日照時間であり、単位は秒であることが分かります。なお、「var0_6_194_surface」は、仕様書によればデータ品質情報です。
まず、日照時間の分布を可視化しましょう。
ds['SUNSD_surface'][0,:,:].plot(cmap='RdYlGn_r')
<matplotlib.collections.QuadMesh at 0x2b9815e6ad0>
推計気象分布（日照時間）のデータは、数値予報モデルのそれとは異なり、海上や小さな島に位置する格子点にはデータがありません。
分布をみると、2023年5月31日23時～6月1日00時(UTC)の間、東日本では広い範囲で晴天である一方、西日本ではほとんど日照がなかった地域も多かったことがわかります。
次に、このデータの品質値を可視化しましょう。
ds['var0_6_194_surface'][0,:,:].plot()
<matplotlib.collections.QuadMesh at 0x2b9816ba4a0>
品質値は、全域にわたって１です。仕様書の表１によれば、値１は、「やや疑わしい」と定義されています。「やや疑わしい」は実用上問題ありません。
さて、今度は、気温データを見てみましょう。1日と15時間分取得してみましょう。今回は日本時間の2023年6月1日00時(UTC)から2日15時(UTC)について読み出します。
まずは、GRIBファイルのリストを作りましょう。数値予報モデルGPVのときと少しプログラムを変えています。これは、観測系のGPVデータの場合、取り出す期間によってはサブフォルダ年や年月をまたぐことがあり、これに対応するためです。
# GRIBファイルへのパスを作る
prdct = 'obs_gpv'
param = 'Ptt'
yyyymmddhhs = wx.strft_range(start='2023060100', end='2023060215', format='%Y%m%d%H', hours=1)
grb_files = []
for yyyymmddhh in yyyymmddhhs:
grb_path = f'./jmadata/{prdct}/{yyyymmddhh[0:4]}/{yyyymmddhh[0:6]}'
grb_files.append(f'{grb_path}/Z__C_RJTD_{yyyymmddhh}0000_OBS_GPV_Rjp_Ggis1km_{param}_A{yyyymmddhh}00_grib2.bin')
リストの最初のGRIBファイルを覗いて、気象要素に使用されている記号を調べてみましょう。
wx.getvarname(grb_files[0])
['TMP']
変な名前でなくて安心しました。
これから、GPVデータを読みこむわけですが、日照時間の概要フォームで確認した通り、推計気象分布のデータはGPVに展開すると非常に大きな配列となるため、GRIBファイルからGPVを読み込む段階で領域を中部地域に限定することにします。
以下を実行してください。
# 全域をオブジェクト化するとサイズが巨大になるので領域を限定して読みこみます
area = [34, 36, 135, 138] #３次メッシュ5135～5337
ds = wx.getgpv(grb_files,'TMP', timezone='Asia/Tokyo', lalomima=area)
T_obg = ds['TMP_surface']
T_obg
<xarray.DataArray 'TMP_surface' (time: 40, latitude: 240, longitude: 240)>
array([[[ nan, nan, nan, ..., nan, nan, nan],
[ nan, nan, nan, ..., nan, nan, nan],
[ nan, nan, nan, ..., nan, nan, nan],
...,
[ nan, nan, nan, ..., 288. , 290. , 290.5],
[ nan, nan, nan, ..., 289. , 290. , 290. ],
[ nan, nan, nan, ..., 289.5, 290. , 288.5]],
[[ nan, nan, nan, ..., nan, nan, nan],
[ nan, nan, nan, ..., nan, nan, nan],
[ nan, nan, nan, ..., nan, nan, nan],
...,
[ nan, nan, nan, ..., 290. , 292. , 292.5],
[ nan, nan, nan, ..., 291. , 292.5, 292. ],
[ nan, nan, nan, ..., 291.5, 292. , 290.5]],
[[ nan, nan, nan, ..., nan, nan, nan],
[ nan, nan, nan, ..., nan, nan, nan],
[ nan, nan, nan, ..., nan, nan, nan],
...,
...
...,
[ nan, nan, nan, ..., 290. , 291. , 291.5],
[ nan, nan, nan, ..., 290.5, 291.5, 291. ],
[ nan, nan, nan, ..., 291. , 291. , 290.5]],
[[ nan, nan, nan, ..., nan, nan, nan],
[ nan, nan, nan, ..., nan, nan, nan],
[ nan, nan, nan, ..., nan, nan, nan],
...,
[ nan, nan, nan, ..., 290. , 291. , 291.5],
[ nan, nan, nan, ..., 290.5, 291.5, 291. ],
[ nan, nan, nan, ..., 291. , 291. , 290. ]],
[[ nan, nan, nan, ..., nan, nan, nan],
[ nan, nan, nan, ..., nan, nan, nan],
[ nan, nan, nan, ..., nan, nan, nan],
...,
[ nan, nan, nan, ..., 289.5, 291. , 291.5],
[ nan, nan, nan, ..., 290. , 291. , 291. ],
[ nan, nan, nan, ..., 290.5, 291. , 290. ]]], dtype=float32)
Coordinates:
* latitude (latitude) float64 34.0 34.01 34.02 34.03 ... 35.98 35.99 36.0
* longitude (longitude) float64 135.0 135.0 135.0 135.0 ... 138.0 138.0 138.0
* time (time) datetime64[ns] 2023-06-01 ... 2023-06-02T15:00:00
awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-03T00:00:0...
Attributes:
short_name: TMP_surface
long_name: Temperature
level: surface
units: K2023年6月1日03時UTCにおける気温分布を可視化しましょう。
T_obg_23060103 = T_obg.sel(time='2023-06-01T03:00')
T_obg_23060103.plot(cmap='RdYlGn_r')
plt.title(T_obg_23060103.awaretime.data)
plt.draw
<function matplotlib.pyplot.draw()>
推計気象分布（気温）のデータは、海上や小さな島に加え、内水面に位置する格子点にもデータがありません。
名古屋地方気象台が位置する3次メッシュの気温の推移を折れ線グラフにしてみましょう。
lat, lon = 35.166667, 136.965 #名古屋地方気象台
T_obg_nag = T_obg.sel(latitude=lat, longitude=lon, method='nearest')
T_obg_nag.plot(x='awaretime',figsize=(10,3))
plt.title(f'(N{T_obg_nag.latitude:.4f}, E{T_obg_nag.longitude:.4f})')
plt.draw
<function matplotlib.pyplot.draw()>
推計気象分布の気温は、0.5度刻みに丸められたGPVです。
6.4 解析雨量¶
このプロダクトは、気象庁のアメダス降水量計、気象レータ、国土交通省河川局・道路局のレーダ・雨量計、外の観測結果を合成して作られる降水量のGPVです。観測結果をGPVにしたものなので、未来のデータはありません。このプロダクトの作成方法の概要は、気象庁ホームページに、データ構造に関する説明は、配信資料に関する技術情報（気象編）第238号 で示されています。
GPVがカバーする緯度経度範囲は北緯20度～48度、東経118度～150度で、格子点間隔は南北方向0.008333度、東西方向が0.0125度です。この格子は、国土数値情報における標準地域メッシュ(3次メッシュ)の中心点にと一致しています。そして、値は3次メッシュの代表的な値とされています。つまり、厳密な意味ではGPVではなくメッシュデータです。
30分毎に、過去60分間の積算降水量のGPVが作られます。30分毎に作成される一方、値は過去60分積算のため、全てのファイルを連結してしまうと降水量は2倍になってしまうので、雨の強弱を知りたいのか、降水量を知りたいのかで、読み込むファイルや処理方法を変える必要があります。
なお、降水は生活や産業活動、気象災害に深くかかわることから、速報版解析雨量 や 高解像度降水ナウキャスト など、様々なGPVプロダクトが別途作られています。
解析GRIBファイル名は以下のようなものです。サンプルデータはフォルダ Challenge3/jmadata/ra に、年と月で整理されて保存されています。
Z__C_RJTD_20230601010000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin
ファイル名は、いくつかの文字列がアンダースコア「 _ 」で区切られた形となっています。14桁の数字の並びがGPVの時刻で 年(4桁)月(2桁)日(2桁)時(2桁)分(2桁)秒(2桁) です。中央協定時で表現されています。Prr60lv_ANAL が降水量解析値であることを示します。
では、推計気象分布と同じ要領で降水量データを読みこんでみましょう。まず、GRIBファイルのリストを作ります。
# GRIBファイルへのパスを作る
prdct = 'ra'
yyyymmddhhs = wx.strft_range(start='2023060100', end='2023060215', format='%Y%m%d%H', hours=1)
grb_files = []
for yyyymmddhh in yyyymmddhhs:
grb_path = f'./jmadata/{prdct}/{yyyymmddhh[0:4]}/{yyyymmddhh[0:6]}'
grb_files.append(f'{grb_path}/Z__C_RJTD_{yyyymmddhh}0000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin')
リストの最初ファイルから、保持されているGPVの記号を取り出します。
print(wx.getvarname(grb_files[0]))
['var discipline=0 center=34 local_table=1 parmcat=1 parm=200']
大変奇妙な名称ですが、GPVデータが一つ見出されました。このデータを読み込み、概要フォームを表示させてみます。
element = ['var discipline=0 center=34 local_table=1 parmcat=1 parm=200']
ds = wx.getgpv(grb_files[0], element)
ds
<xarray.Dataset>
Dimensions: (latitude: 3360, longitude: 2560, time: 1)
Coordinates:
* latitude (latitude) float64 20.0 20.01 20.02 ... 47.98 47.99 48.0
* longitude (longitude) float64 118.0 118.0 118.0 ... 150.0 150.0
* time (time) datetime64[ns] 2023-06-01
Data variables:
var0_1_200_surface (time, latitude, longitude) float32 nan nan ... nan nan
Attributes:
Conventions: COARDS
History: created by wgrib2
GRIB2_grid_template: 0正式名称(long_name)が「desc」、単位(units)が「unit」と読み出されました。付随情報については正しく読み取れていません。ただ、Datasetオブジェクトにおける記号が var0_1_200_surface と分かったので最初このGPVを読み込んで分布図を作ってみましょう。
da = ds['var0_1_200_surface']
da[0,:,:].plot(cmap='Blues') #青のグラデーションにするには 'Blues'
<matplotlib.collections.QuadMesh at 0x2b981a1b430>
解析雨量においても、データがない領域が存在します。データが存在する領域の形は、もととなるレーダーの観測域を合成した形になっています。データの領域には、台風3号によるとみられる同心円状の降水域のような模様が見られます。
しかし、緑～赤のカラーチャートは降水量分布の表現にはちょっとふさわしくありませんね。水っぽく、青のグラデーションで書き直してみましょう。上のスクリプトの2行目、'RdYlGn_r' を 'Blues' に書き直して、Cellを再実行してみてください。
Matplotlibで使用できるカラーチャートは、これらのほかにもさまざまなものがあります。他の配色や自作の仕方などについては、下記を参照してください。
https://matplotlib.org/stable/users/explain/colors/colormaps.html
降水量の話に戻ります。「 unit 」なる単位のデータですが、アメダス観測値の比較の結果、単位を kg/m2 (mm) とする降水量で間違いないことが確認されています。
それでは、推計気象分布と同じように、GRIBファイルからGPVを読み込む段階で領域を中部地域に限定してファイルからGPVを読み込みます。
また、付随情報を付け直します。
以下を実行してください。
# 全域をオブジェクト化するとサイズが巨大になるので領域を限定して読みこみます
area = [34, 36, 135, 138] #３次メッシュ5135～5337
ds = wx.getgpv(grb_files,'var discipline=0 center=34 local_table=1 parmcat=1 parm=200',
timezone='Asia/Tokyo', to_netcdf=False, lalomima=area)
P_ra = ds['var0_1_200_surface']
P_ra.attrs = {'long_name':'Total precipitation', 'units':'kg/m^2','level':'surface'}
それでは、中部地域の分布図を描いてみましょう。他のプロダクトでは、2023年6月1日03時(UTC)の分布図を作りましたが、その時刻、この地域に降水はなかったので、一日後、6月2日03時(UTC)の分布図を作ります。
P_ra_23060203 = P_ra.sel(time='2023-06-02T03:00')
P_ra_23060203.plot(cmap='Blues')
plt.title(P_ra_23060203.awaretime.data)
plt.draw
<function matplotlib.pyplot.draw()>
次に、名古屋地方気象台が位置する3次メッシュの降水量の推移を折れ線グラフにしてみましょう。
lat, lon = 35.166667, 136.965 #名古屋地方気象台
P_ra_nag = P_ra.sel(latitude=lat, longitude=lon, method='nearest')
P_ra_nag.plot(x='awaretime',figsize=(10,3))
plt.title(f'(N{P_ra_nag.latitude:.4f}, E{P_ra_nag.longitude:.4f})')
plt.draw
<function matplotlib.pyplot.draw()>
付録 分布図の地図投影¶
ご注意
地図投影や海岸線の追加を行うには、ライブラリ cartopy が必要です。皆さんにお配りした「事前設定ガイド」には、インストールすべきライブラリにこれを指定していなかったので、まだ皆さんの Python 環境では利用できません。以下を実行し、 cartopy をインストールしてください。
- スタートボタンから、Anconda 3 を開き Anaconda Prompt を開始する
- 「conda install cartopy」と打ち込み、インストールを開始する
- 確認を求められたら、[y] キーを押す
- 「conda list cartopy」と打ち込み、インストールを確認する
- 「exit」と打ち込み、 Anaconda Prompt を終了する
4.3章(1) Fancy indexing においては、plt.title('New Title') 等によって、タイトルやラベルは簡単に付け替えられることを学びましたが、この方法はいわば簡便法であり図の複雑な加工には適しません。
分布図を地図投影するような複雑な加工は、「オブジェクト指向のプロッティング」と呼ばれる違う作法で行います。この作法では、図を、キャンバス、図柄、軸などの構成要素に分けて考え、「まずキャンバスを定義する、図柄(散布図など)を一つ持つ、図柄は軸を持つ、軸はタイトルを持つ、・・。」などのように、図を構成する要素とそれが持つ付随物を一つ一つ設定して、全体を設定します。
では、6.4 全球数値予報モデルGPV（日本域） で描画した中部地域の気温分布図を地図投影し海岸線を付加することを通し、その方法を学びましょう。
まず、Pythonカーネルをクリアして、改めて開始することにしましょう。ツールバーの「輪矢印」アイコンをクリックし、[Restart] ボタンをクリックして今実行中のPythonカーネルをリセットし、その後に下記のインポート文を実行してください。
# ライブラリ のインポート
import numpy as np
import wxbcgribx as wx
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
まず、GSM-GPVのGRIBファイルから気温分布データを読み込みます（変数名は必ずしも 6.2 と同じではありません）。
# プロダクトの読み込み、GPVの切り出し
grb_path = f'./jmadata/gsm/2023/202306'
grb_file = f'{grb_path}/Z__C_RJTD_20230601000000_GSM_GPV_Rjp_Gll0p1deg_Lsurf_FD0000-0100_grib2.bin'
element = ['TMP']
ds = wx.getgpv(grb_file, element, timezone='Asia/Tokyo')
area = [34, 36, 135, 138] #１次メッシュ5135～5337
da = ds['TMP_2maboveground']
var = da.sel(time='2023-06-01T03:00',
latitude=slice(area[0],area[1]),
longitude=slice(area[2],area[3]) )
参照用に xarry.plot メソッドだけで可視化しておきましょう。
var.plot(cmap='RdYlGn_r')
<matplotlib.collections.QuadMesh at 0x2b9c4721db0>
この図に対し、オブジェクト指向の作法で、地図投影、海岸線描画、タイトル変更 を行います。
fig = plt.figure(figsize=(9,5)) # 横8×縦5のキャンバスを定義
ax = fig.add_subplot( # キャンバスの中に図柄を追加
1,1,1, # 1行1列1枚の図柄
projection=ccrs.PlateCarree() # この図柄を正距円筒図法で地図投影する
)
var.plot( # plotメソッドで分布図作成
ax=ax, # 分布図を図柄として引き渡す
transform=ccrs.PlateCarree(), # この図柄は正距円筒図法で緯度経度を結びつくと宣言
cmap='RdYlGn_r'
)
ax.coastlines(linewidth=0.8) #図柄に海岸線を書き込む
gl = ax.gridlines(crs=ccrs.PlateCarree(), linewidth=0.5, linestyle='--', draw_labels=True) # 図柄に緯経線を描くと宣言
gl.xlocator = mticker.FixedLocator(np.arange(132, 140, 1)) # 経線の設定
gl.ylocator = mticker.FixedLocator(np.arange(33, 38, 1)) # 緯線設定
ax.set_title(var.awaretime.data) #図柄のタイトルを設定しなおす
plt.draw()
海岸線を描くためには、投影法を定めておく必要があります。気象庁のGPVプロダクトは緯度経度に対して等間隔に格子を割り当てているので、正距円筒図法(plate carree) で無理なく地図投影が可能なことから、これを使っています。
さて、この図を描くスクリプトには、正距円筒図法の指定が3回出てきました。5, 10, 16行目です。これらのうち、最初の指定と、2番目以降の指定では役割が異なることに注意してください。
最初の指定(5行目)は、図柄をどの図法でキャンバスに投影するかを指定しています。これに対し、２番目(10行目)と３番目(16行目)の指定は、分布図、あるいは、グリッド線が、それぞれどの投影法に基づいて作られているのかを宣言するものです。
したがって、例えば、気象庁のGPVデータをランベルト正角円錐図法で図化する場合には、最初の指定をランベルト正角円錐図法としますが、緯度経度でGPVの格子点が並ぶ気象庁のGPVを可視化する限りは、２番目では常に正距円筒図法を宣言します。３番目は、正距円筒図法における直線を描く指定なので、これについても正距円筒図法を宣言することになります。
上の図を、ランベルト正角円錐図法で描き直す例を下に示します。図法の違いを示すため、図幅をやや広くとっています。
fig = plt.figure(figsize=(10,4)) # 横10×縦4のキャンバスを定義
ax = fig.add_subplot( # キャンバスの中に図柄を追加
1,1,1, # 1行1列1枚の図柄
projection=ccrs.LambertConformal(central_longitude=135, central_latitude=35) # この図柄をランベルト正角円錐図法に投影する
)
ax.set_extent([130, 140, 33, 38]) # 描画範囲を設定
var.plot(
ax=ax, # plotメソッドの図を図柄として引き渡す
transform=ccrs.PlateCarree(), # この図柄は正距円筒図法で緯度経度と結びつくと宣言
cmap='RdYlGn_r'
)
ax.coastlines(linewidth=0.8) #図柄に海岸線を書き込む
gl = ax.gridlines(crs=ccrs.PlateCarree(), linewidth=0.5, linestyle='--', draw_labels=True,color='blue') # 正距円筒図法上の直線を引くと宣言
gl.xlocator = mticker.FixedLocator(np.arange(130, 142, 2)) # 経度線の設定
gl.ylocator = mticker.FixedLocator(np.arange(33, 38, 2)) # 緯度線の設定
ax.set_title(var.awaretime.data) #図柄のタイトルを設定しなおす
plt.draw()
分布図を気象学的に 読みたい 場合には、等値線図が有効な場合があります。等値線図は、DataArrayオブジェクトのメソッド plot.contour で簡単に描くことができます。以下に、GSM-GPVの海面気圧の等値線図を描くスクリプトの例を示します。
# プロダクトの読み込み、GPVの切り出し
grb_path = f'./jmadata/gsm/2023/202306'
grb_file = f'{grb_path}/Z__C_RJTD_20230601000000_GSM_GPV_Rjp_Gll0p1deg_Lsurf_FD0000-0100_grib2.bin'
element = ['PRMSL'] # 平均海面高度気圧
ds = wx.getgpv(grb_file, element, timezone='Asia/Tokyo')
da = ds['PRMSL_meansealevel']
var = da.sel(time='2023-06-01T03:00') # 空間的な切り出しは行わない
# 作図
fig = plt.figure(figsize=(8,8)) # 横8×縦8のキャンバスを定義
ax = fig.add_subplot( # キャンバスの中に図柄を追加
1,1,1, # 1行1列1枚の図柄
projection=ccrs.LambertConformal(central_longitude=135, central_latitude=35) # この図柄をランベルト正角円錐図法に投影する
)
ct = var.plot.contour(
ax=ax, # メソッドの図を図柄として引き渡す
transform=ccrs.PlateCarree(), # この図柄は正距円筒図法で緯度経度と結びつくと宣言
cmap='RdYlGn_r',
levels = 21 # とりあえずレベル分けは21
)
ax.clabel(ct, inline=True, fontsize=10) # メソッドの図に等値線を書き込む
ax.coastlines(linewidth=0.8) #図柄に海岸線を書き込む
gl = ax.gridlines(crs=ccrs.PlateCarree(), linewidth=0.5, linestyle='--', draw_labels=True,color='blue') # 正距円筒図法上の直線を引くと宣言
gl.xlocator = mticker.FixedLocator(np.arange(120, 150+1, 10)) # 経度線の設定
gl.ylocator = mticker.FixedLocator(np.arange(20, 50+1, 10)) # 緯度線の設定
ax.set_title(f'{var.long_name} [{var.units}]\n{var.awaretime.data}') #図柄のタイトルを設定しなおす
plt.draw()
この例では、GRIBファイルに記録されている海面気圧をそのままプロットしたので単位が Pa です。気象学では hPa を使用するのが一般的なので、DattaArray オブジェクト var の値を 1/100 したうえで図化した方がよりよいでしょう。その方法は今回の研修で習得していますので、興味のある方はチャレンジしてみてください。なお、等値線の間隔や個々の値、色や太さなど、等値線図には様々な設定項目があります。それらについては、
- xarray.plot.contorの説明文書
https://xarray.pydata.org/en/v2023.09.0/generated/xarray.DataArray.plot.contour.html や、 - matplotlib.pyplot.clabelの説明文書
https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.clabel.html などを参照してください。
おわりに¶
この研修では、気象庁が様々なGPVデータを対外的に作成していること、それらGPVプロダクトは、GRIB2というフォーマットのファイルに入れられていること、産業利用を検討する目的で過去のGPVを入手する方法は、現時点では限られていることを学習しました。
また、WXBCのオリジナルライブラリ wxbcgribX を利用すると、GRIBファイルからGPVを取り出すことができ、それはxarray.DataArrayオブジェクトとなること、このオブジェクトを操作することにより、気象値の演算や可視化ができることを学習しました。
そして、この研修を通じて、気象庁が保有する3つの数値予報モデル、LFM、MSM、GSM のGPV、ならびに、観測結果ベースのプロダクト、推計気象分布、解析雨量のGPVを実際に読み込んで可視化しました。
この研修では、複数の気象庁GPVプロダクトを実際触っていただくことを大きな目的としたため、処理テクニックの説明については、最低限にとどめています。実際に利用する上では、GPVデータの集計(平均など)、補間(3時間間隔のデータを1時間間隔に)、リグリッド(プロダクトで異なる格子座標を揃える)、結合(複数のGPVを一つのGPVにまとめる)などのテクニックが必要となりますが、これについては、「気象庁GPVデータ分析チャレンジ！基礎編」で取り扱う予定です。さらにスキルを深めたい方はこちらにもご参加ください。
著作権について¶
(C) 2024 WXBC
＜利用条件＞
本書は、本書に記載した要件・技術・方式に関する内容が変更されないこと、および出典を明示いただくことを前提に、無償でその全部または一部を複製、翻案、翻訳、転記、引用、公衆送信等して利用できます。なお、全体または一部を複製、翻案、翻訳された場合は、本書にある著作権表示および利用条件を明示してください。
＜免責事項＞
本書の著作権者は、本書の記載内容に関して、その正確性、商品性、利用目的への適合性等に関して保証するものではなく、特許権、著作権、その他の権利を侵害していないことを保証するものでもありません。本書の利用により生じた損害について、本書の著作権者は、法律上のいかなる責任も負いません。