気象要素の地図上の表示としてよく知られている 地上天気図は、等圧線、つまり、海面気圧の等値線をふくんでいる。 ここでは、それを Python で作図してみよう。 (天気図には、等圧線のほかに、低気圧・高気圧の中心、前線、各観測点での観測値も表示されている。 今回の実習では、そこまではやらない。したがって「天気図をつくる実習」にはなっていない。)
「【...】」のように角かっこをつけた部分は読みとばしてもよい。
Pythonプログラムの例文は、luce の /home/masuda/shared/program/ に置いてある。
各自、自分のホームディレクトリの下の作業用ディレクトリ (前回作成した python_map_practice をまた使ってもよい) に cd してから、 つぎのコマンドによって、例文をそこにコピーしてほしい (スペースのあとのピリオドが重要)。
$ cp /home/masuda/shared/program/*.py .
読まれる気象データは、/home/masuda/shared/jra55/surf125/ に置いてある。
今回のプログラムは、そこに置かれたデータをそのまま読むので、 データを自分のところにコピーする必要はない。 (もし、別のところに置かれたデータを読むばあいは、 プログラムのなかの datapath という変数にあたえる文字列を変更する。)
作図を実行するには、あらかじめ Xming を起動してから Teraterm で luce に login し、 作業用ディレクトリに cd して、 プログラムを必要に応じて編集したあと、つぎのような形で実行する。
$ python3 なになに.py
図を画面に表示したいときは「plt.show()」を、 画像ファイルに保存したいときは「plt.savefig(...)」をふくむ4行を、それぞれ生かし、 つかわないほうは行の頭に「#」をつけてコメントにしておく。 例文のその部分が自分の意図とあっているか確認し、ちがっていたら編集で変更してから実行してほしい。 画像ファイルは、作業用ディレクトリの下の「image」というディレクトリの中につくられる。 画像ファイルに保存したいときは、あらかじめ「mkdir image」で そのディレクトリをつくっておいてから実行する。
格子点型の気象データの例として、日本の気象庁がつくった再解析データ「JRA-55」を利用する。 これは、観測値の統計的内挿と数値予報モデルをくみあわせた「データ同化」技術の産物であり、 1958年から最近まで、6時間間隔で、全球の気象要素の数値を収録している。
再解析データは、3次元の空間分布をもつ量がおもになっているが、 ここでとりあげる海面気圧は、地表面付近の水平2次元の空間分布をもつ量を 緯度・経度 1.25度の格子で収録した「anl_surf125」というサブセットにふくまれている。
【配布されている JRA-55 データの形式は、 WMO (世界気象機関) で決められた GRIB という形式の 第1版 である。 (現行の標準である第2版ではない。) Python では、pygrib というパッケージを入れれば、GRIB 形式のデータを直接読むことができるが、 わたしはまだ pygrib をつかった経験がなく、 また、ここでは一般に格子点型データを読む作業の例題をしめしたいので、 あらかじめ単純バイナリデータに変換しておいたものを読むことにする。 変換には GrADS というソフトウェアを利用した。 変換の手順は、別ぺージ [JRA-55のGRIBデータをGrADSで読み、単純バイナリデータに変換する]に書いた。】
この実習で読むデータは、 ひとつの気象要素の1時刻の値がひとつのファイルになっていて、 数値が4バイトの浮動小数点型であらわされ、 ファイルの大きさは、{格子点の個数 × 4} バイトである。
どんな気象要素のいつの値であるかの情報は、ファイルのなかにはなく、ファイル名にある。 ファイル名は、たとえば「psl2007081606.r4l」のようにした。 ここで「psl」は「海面気圧」という気象要素をあらわす。 「2007081606」は日時で、 この例は、2007年 8月 16日の 世界時 (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 からである。 複数のプログラム言語をつかう人は注意する必要がある。】
データの一例として、あとで図をしめす、 2007年1月7日 世界時0時 (日本時間9時) の海面気圧のデータを置く。
luce には、2007年の毎日4回 (世界時 0, 6, 12, 18時) のデータを置いた。 【2007年をとりあげたのは、前の教材の21節 (プログラム例 16) のアメダスの気温データの日をふくめたかったからだが、 等圧線の事例としては、8月16日のものはわかりにくかったので、温帯低気圧が発達した日をえらんだ。】
海面気圧は、実際の地表面気圧から、「海面更正」という補正手続きによって得られるものである。
英語では "mean sea level pressure"、"sea level pressure" の両方の表現が よくみられる。ここでは mean がないほうを採用する。 ファイル名などにつかう略号は (「気圧」であることをさきにだして) 「psl」とする。
【地表面の標高が高いところでは、海面更正の不確かさが大きいので、 地上天気図よりも 850 hPa や 700 hPa などの高層天気図を見るのがよいとされている。 高層天気図に出てくる量は、JRA-55 では3次元の気象場を収録した 「anl_p125」というサブセットにある。 残念ながら、この教材ぺージではあつかわない。】
全球のデータセットを読みはじめるときには、 格子点の配置が期待どおりであることを確認する必要もあるので、 まず世界地図上に分布を表示してみる。
地球を平面の画面に表示するにあたって、どの経度で切断するかという問題がある。 データにすなおにしたがえば、経度 0度 (グリニッジ子午線) が端になる。 ヨーロッパやアメリカの人がつくった例文では、グリニッジ子午線を中央にもってきて、 経度 180度 (日付変更線) を端にすることが多い。 わたしは、太平洋を分割したくないし、アフリカ大陸を分割したくないので、 西経 30度 (ここでのプログラムでは東経を正とするので、経度 -30度 または 330度)を端、 東経 150度を中央にすることにした。
等値線 (気圧なので「等圧線」でもある) を線でかくこともやってみたが、 気圧の変動幅が、温帯で大きく、熱帯で小さいので、 等値線間隔をどのように設定しても、世界の気圧分布があまり見やすくならない。
色がつかえるので、contourf という機能を利用して、 気圧の値を色でぬりわけた図をつくることにする。
まず、1つの時刻の図をかくプログラム例 18a をしめす。
【[2021-07-20 授業後補足] 授業時に用意したプログラム例 18a にまちがい(修正もれ)がありました。 「plt.show()」の行のはじめに空白がありますが、このままだとエラーになります。 空白をつめて p の字が行の左端にくるのが正しいです。このページからリンクされたプログラム例は修正しました。】
↓ 結果の例 psl2007010700_global.png (2007年 1月 7日 世界時0時の海面気圧)
複数の時刻の同様な図がほしいとき、 1時刻ごとにプログラムを修正して起動するのはてまがかかりすぎるので、 プログラムを、時刻についてのループをふくむ構造にして、図をつづけてつくる。 年にまたがることはないとして、月、日、時の3重ループにした。 Python では、ループの内容を右に (標準的には4文字ぶん) 段下げ (indent) する必要があるので、 まえの例文と見かけがすこし変わってくるが、 プログラム例 18b は、それぞれのループを1回だけ実行するようにしたので、 実行内容は 18a とかわらない。
データは6時間ごとにあるので、1年ぶんだと1460枚の図になるが、 1回の作業で実際に画面で見るのに適した図の数は 数十枚である。 プログラム例 18c は、月単位で期間をえらび、その月の毎日1回 (世界時 0時) の気圧の分布図をつづけてつくる。 なお、月の日数は calendar パッケージの monthrange 関数で得ている。 【[2021-07-23 補足] monthrange を呼び出したところにある「[1]」は、 0が月のはじめ、1が月のおわりという意味であり、月の番号ではない。】 「for month in [8]:」のカッコの中の数値を見たい月に変えて実行してみてほしい。
luce には 「ImageMagick」がインストールされているので、 「display」や「animate」などのコマンドがつかえる。 プログラムの結果を画像ファイルに保存するようにし、 image ディレクトリに cd して、つぎのようにして、動画の形で見てみよう。 (2007につづく文字は作図対象の月によって変わる。)
$ animate psl200708*_global.png
(この節は、くわしいことはわからなくてもかまわないので、概略をつかんでほしい。)
地図の領域を、北緯 20～50度、東経 120～150度にかぎり、 「contour」 (「f」がつかない) という関数を利用して、 気圧の等値線 (つまり「等圧線」) をかいてみる。 間隔は、日本の地上天気図にあわせて、4 hPa にした。 【色ぬりはしないので、contourf や colorbar はつかわない。】
プログラム例 19 の構成は、18b と同様、月、日、時の3重ループにしたうえで、 それぞれのループを1回だけ実行することによって、図を1つだけかくようにしてある。 「for month in [1]:」「for day in [7]:」「for hour in [0]:」の in のあとにくるものを 変えることによって、作図対象の日時を変更することができる。
地図の領域を限定するのは、ax.set_extent でおこなっている。 緯度・経度めもりの位置をあたえる xticks, yticks の値も、それにあわせて変えている。
【[2021-07-20 授業後補足] 授業時に用意した プログラム例 19 に修正もれがありました。 現状はコメントになっているのでエラーではありませんが、有効にするとエラーになります。 「for day ...」のループのところで、変更案のコメントとして 「for day in range(daystart, dayend+1)」と書いてありますが、 ループのはじめの文ですから、かっことじのあとに、コロン(「:」) があるのが正しいです。 このページからリンクされたプログラム例は修正しました。】
↓ 結果の例 psl2007010700_contour_japan.png (2007年 1月 7日 世界時0時 (日本時間 9時) の海面気圧)
格子点の位置がどこにあり、そこにどのような数値がはいっているかを見たいこともある。 それを図示するプログラム例 20 をつくった。 地点の位置に、記号「+」を plot で、気圧の数値を文字列にしたものを text で表示するようにした。
地図の領域は、プログラム例 19 よりもさらに限定して、北緯 34.5～38度、東経 138～142度とした。
【緯度値の目盛りは整数値のほうがわかりやすいので、yticks には latmin (= 34.5) のかわりに 35.0 をあたえることにした。】日時の指定のしかたは、プログラム例 18b, 19 と同様である。
【[2021-07-20 授業後補足] プログラム例 20 も、プログラム例 19 と同様に、 「for day ...」のループのところで、変更案のコメントとして書いてある 「for day in range(daystart, dayend+1)」のかっことじのあとに、コロン(「:」) があるのが正しいです。 このページからリンクされたプログラム例は修正しました。】
↓ 結果の例 psl2007010700_contour_kanto.png (2007年 1月 7日 世界時0時 (日本時間 9時) の海面気圧)
2007年のうちで見たい期間 (1週間から1か月程度) をえらび、 プログラム例 19 のループのまわしかたをプログラム例 18c にならって変更して、 毎日1回の等圧線図を作図して、観察してみよう。
(たとえば、ある月の5日から12日ならば、monthの[]の中にその月の番号を入れ、 dayのループを「for day in range(5,12+1):」とする。 複数の月にまたがるが月の全部ではない期間について作図したいときは、 プログラム例をそのような期間に対応させるのが簡単でないので、 月ごとにわけて実行してほしい。)
(時間間隔が1日ではあらすぎると感じたばあいは、6時間間隔で見てみよう。 hourのループを「for hour in [0, 6, 12, 18]:」とする。)
地図の領域を変えてもよい。 そのばあい、等値線をひく値の設定 (levels にあたえる値) も変える必要があるかもしれない。 たとえば、熱帯の気圧分布を見るばあいは、等圧線間隔をこまかくする必要があるかもしれない。 温帯低気圧が発達する海域では、等値線をひく気圧の下限をもっと低くする必要があるだろう。
2007年の各月の天候の特徴についての参考として、気象庁のウェブサイトのつぎのページをあげておく。