気象要素の地図上の表示としてよく知られている 地上天気図は、等圧線、つまり、海面気圧の等値線をふくんでいる。 ここでは、それを Python で作図してみよう。 (天気図には、等圧線のほかに、低気圧・高気圧の中心、前線、各観測点での観測値も表示されている。 今回の実習では、そこまではやらない。したがって「天気図をつくる実習」にはなっていない。)
「【...】」のように角かっこをつけた部分は読みとばしてもよい。
Pythonプログラムの例文は、luce の /home/masuda/contour/ に置いてある。 各自、つぎのコマンドによって、例文を 自分のホームディレクトリの下にコピーしてほしい。
$ cp -r /home/masuda/contour .
読まれる気象データは、/home/masuda/jra55/daily/surf125/ の下に置いてある。 今回のプログラムは、そこに置かれたデータをそのまま読むので、 今学期中に luce で実行するかぎりは、データを自分のところにコピーする必要はない。 【もし、別のところに置かれたデータを読むばあいは、 プログラムのなかの datapath という変数にあたえる文字列を変更する。 また、データのディレクトリ構成にあわせて、 filepathname という変数の値をあたえる実行文も変更する必要があるかもしれない。】
作図を実行するには、あらかじめ Xming を起動してから Teraterm で luce に login し、 作業用ディレクトリ (今回は、名まえをかえていなければ contour) に cd して、 プログラムを必要に応じて編集したあと、つぎのような形で実行する。
$ python3 なになに.py
図を画面に表示したいときは「plt.show()」を、 画像ファイルに保存したいときは「plt.savefig(...)」をふくむ5行 (頭に「##」をいれてある)を、それぞれ生かし、 つかわないほうは行の頭に「#」をつけてコメントにしておく。 例文のその部分が自分の意図とあっているか確認し、ちがっていたら編集で変更してから実行してほしい。 画像ファイルは、作業用ディレクトリの下の「image」というディレクトリの中につくられる (このディレクトリ名はプログラム中の「imagefilepath」という変数で指定してあり、変更可能)。 画像ファイルに保存したいときは、あらかじめ「mkdir image」で そのディレクトリをつくっておいてから実行する。
格子点型の気象データの例として、日本の気象庁がつくった再解析データ「JRA-55」を利用する。 これは、観測値の統計的内挿と数値予報モデルをくみあわせた「データ同化」技術の産物であり、 1958年から最近まで、6時間間隔で、全球の気象要素の数値を収録している。 このデータは (プロジェクトのウェブサイトではなく他の) いくつかのサイトから配布されているが、 増田は DIAS から取得した。
【その後、「JRA-3Q」という再解析プロジェクトが進行しており、 1947年9月から現在までを対象としているが、これまでに公開されたデータの対象期間は 1991年から 2013年4月までなので、 今学期のセミナーの実習教材とするのは見送った。】
再解析データは、3次元の空間分布をもつ量がおもになっているが、 ここでとりあげる海面気圧は、地表面付近の水平2次元の空間分布をもつ量を 緯度・経度 1.25度の格子で収録した「anl_surf125」というサブセットにふくまれている。
【配布されている JRA-55 データの形式は、 WMO (世界気象機関) で決められた GRIB という形式の 第1版 である。 (なお、JRA-3Q の配布形式は GRIB の 第2版である。) Python では、pygrib というパッケージを入れれば、GRIB 形式のデータを直接読むことができるが、 わたしはまだ pygrib をつかった経験がなく、 また、ここでは 一般的な格子点型データを読む作業の例題をしめしたいので、 あらかじめ単純バイナリデータに変換しておいたものを読むことにする。 変換には GrADS というソフトウェアを利用した。 変換の手順は、別ぺージ [JRA-55のGRIBデータをGrADSで読み、単純バイナリデータに変換する]に書いた。】
この実習で読むデータは、 ひとつの気象要素の1時刻の値がひとつのファイルになっていて、 数値が4バイトの浮動小数点型であらわされ、 ファイルの大きさは、{格子点の個数 × 4} バイトである。
どんな気象要素のいつの値であるかの情報は、ファイルのなかにはなく、ファイル名にある。 ファイル名は、たとえば「psl2021082406.r4l」のようにした。 ここで「psl」は「海面気圧」という気象要素をあらわす。 「2021082406」は日時で、 この例は、2012年 8月 24日の 世界時 (UTC) の6時 (日本時間 (JST) では15時) である。 「r4l」(最後の文字はエルの小文字)は、たまたまそうなっていると思ってほしい。 【これは、増田が個人的にデータの形式を思いだすためにつけた符丁で、 4バイト浮動小数点数 (floating-point number なのだが、Fortran言語でいう 「実数型」(REAL) という表現のなごりで r となっている)、 下位バイト先行 (least signifcant byte first, little-endian) であることをしめす。】
数量の単位は、ファイルの内容にもファイル名にもふくまれていない。 気圧の単位は、JRA-55 では (hPa ではなく) Pa である。 単純バイナリデータに変換する際にも変えていない。
格子点のならびかたは、各緯度で東西方向にならんだ行が、南北方向にならんだ形である。 東西方向には、東経 0度 (グリニッジ子午線) からはじまり、東にすすみ、288個の格子点がある。 南北方向には、南緯 90度 (南極) からはじまり、北にすすみ、145の行がある。 【この格子点の順序は、ちかごろの全球の格子点データではふつうのものだが、 かならずしもすべてのデータセットがこの順序になっているわけではない。 南北は北からはじまるものも多い。東西は経度「-180度」 (日付変更線の経度) からはじまるものも多い。】
これを Python の NumPy の2次元配列に読みこむと、 配列の1つめの引き数が 南北方向の位置、2つめの引き数が 東西方向の位置をしめす番号になる。 【Fortran で読むと、1つめが東西方向、2つめが南北方向である。 また、番号のかぞえはじめが、Python では 0 から、Fortran では 1 からである。 複数のプログラム言語をつかう人は注意する必要がある。】
luce には、2021年1月1日から2023年6月10日までの毎日4回 (世界時 0, 6, 12, 18時) のデータを置いた。 【例としては、関東に雪がふった 2023年2月10日 世界時 0時 (日本時間 9時) をあげている。】
海面気圧は、実際の地表面気圧から、「海面更正」という補正手続き (もし地表面と平均海面の高さのあいだに空気があったら平均海面の高さでの気圧はどうなるかの推定) によって得られるものである。
英語では "mean sea level pressure"、"sea level pressure" の両方の表現が よくみられる。ここでは mean がないほうを採用する。 ファイル名などにつかう略号は (「slp」とされることも多いが、ここでは「圧力」であることをさきにだして) 「psl」とする。
【地表面の標高が高いところでは、海面更正の不確かさが大きいので、 地上天気図よりも 850 hPa や 700 hPa などの高層天気図を見るのがよいとされている。 高層天気図に出てくる量は、JRA-55 では3次元の気象場を収録した 「anl_p125」というサブセットにある。残念ながら、この教材ぺージではあつかわない。】
全球のデータセットを読みはじめるときには、 格子点の配置が期待どおりであることを確認する必要もあるので、 まず世界地図上に分布を表示してみる。
地球を平面の画面に表示するにあたって、どの経度で切断するかという問題がある。 データにすなおにしたがえば、経度 0度 (グリニッジ子午線) が端になる。 ヨーロッパやアメリカの人がつくった例文では、グリニッジ子午線を中央にもってきて、 経度 180度 (日付変更線) を端にすることが多い。 わたしは、太平洋を分割したくないし、アフリカ大陸を分割したくないので、 西経 30度 (ここでのプログラムでは東経を正とするので、経度 -30度 または 330度)を端、 東経 150度を中央にすることにした。
等値線 (気圧なので「等圧線」でもある) を線でかくこともやってみたが、 気圧の変動幅が、温帯で大きく、熱帯で小さいので、 等値線間隔をどのように設定しても、世界の気圧分布があまり見やすくならない。
色がつかえるので、contourf という機能を利用して、 気圧の値を色でぬりわけた図をつくることにする。
まず、1つの時刻の図をかくプログラム例 23a をしめす。
複数の時刻の同様な図がほしいとき、 1時刻ごとにプログラムを修正して起動するのはてまがかかりすぎるので、 プログラムを、時刻についてのループをふくむ構造にして、図をつづけてつくる。 年にまたがることはないと仮定して、月、日、時の3重ループにした。 Python では、ループの内容を右に (標準的には4文字ぶん) 段下げ (indent) する必要があるので、 まえの例文と見かけがすこし変わってくるが、 プログラム例 23b は、それぞれのループを1回だけ実行するようにしたので、 実行内容は 23a とかわらない。
データは6時間ごとにあるので、1年ぶんだと1460枚の図になるが、 1回の作業で実際に画面で見るのに適した図の数は 数十枚である。
プログラム例では、時 (hour) のループを「for hour in [0]:」として、 毎日1回 世界時 00時 (日本時間 09時) についてだけ作業するようにした。 コメントにしてある「for hour in [0, 6, 12, 18]:」のほうを生かせば 1日4回のデータについて作業するようになる。 [ ] の中を変更して時刻をえらぶこともできる。
月、日については、あらかじめあたえた monthstart の月の daystart の日から monthend の月の dayend の日まで毎日作図するようにした。 月にまたがるばあいは、前の月の終わりの日から、つぎの月の1日につづく。 月の終わりの日付は、calendar パッケージを import して 「calendar.monthrange(year, month)[1]」で得ている。([0]がはじめ、 [1]が終わりである。)
月の全部の日付について作業したいばあいは、 「if month == monthstart:」から「for day in range(daystartofloop, dayendofloop+1):」 までを消して (あるいはコメントにして)、 コメントになっている「for day in range(1, lastdayofmonth+1):」を生かせばよい。
luce には 「ImageMagick」がインストールされているので、 「display」や「animate」などのコマンドがつかえる。
プログラムの結果を画像ファイルに保存するようにし、 image ディレクトリに cd して、つぎのようにして見てみよう。
動画の形で見るばあい
$ animate -delay 100 *.png 【授業後に例文を変更】
静止画の形で見るばあい (2023-06-14 追記)
$ display *.png 【授業後に例文を追加】
(この節は、くわしいことはわからなくてもかまわないので、概略をつかんでほしい。)
図の範囲を、東経 110～170度、北緯 10～60度に変更してみる。 気圧と色との関係も、日時の指定のしかたも、23b と同様。
地図の領域を限定するのは、ax.set_extent でおこなっている。 緯度・経度めもりの位置をあたえる xticks, yticks の値も、それにあわせて変えている。
気圧の等値線 (つまり「等圧線」) をかく。 「contour」 (fがつかない) という関数をつかう。 間隔は、日本の地上天気図にあわせて、4 hPa にした。 【色ぬりはしないので、contourf や colorbar はつかわない。】
格子点の位置がどこにあり、そこにどのような数値がはいっているかを見たいこともある。 それを図示するプログラム例 23e をつくった。 地点の位置に、記号「+」を plot で、気圧の数値を文字列にしたものを text で表示するようにした。
地図の領域は、プログラム例 23d よりもさらに限定して、 北緯 30～40度、東経 132～142度とした。 【プログラム名を cjapan としたのは便宜上のもので、ここが日本の中央部であると主張するつもりはない。】
2021年1月から2023年6月10日までのうちで見たい期間 (1週間から1か月程度) をえらび、 プログラム例 23d のループのまわしかたを変更して、 東アジア領域の毎日1回の等圧線図を作図して、観察してみよう。
たとえば、ある月の5日から12日ならば、 monthstart, monthend にはいずれもその月の番号を入れ、 「daystart = 5」 「dayend = 12」とする。
(時間間隔が1日ではあらすぎると感じたばあいは、6時間間隔で見てみよう。 hourのループを「for hour in [0, 6, 12, 18]:」とする。)
地図の領域を変えてもよい。 そのばあい、等値線をひく値の設定 (levels にあたえる値) も変える必要があるかもしれない。 たとえば、熱帯の気圧分布を見るばあいは、等圧線間隔をこまかくする必要があるかもしれない。 温帯低気圧が発達する海域では、等値線をひく気圧の下限をもっと低くする必要があるだろう。
天候の特徴についての参考として、気象庁のウェブサイトのつぎのページをあげておく。