格子点型の気象データのつかいかたは、 各時刻、もし3次元ならば鉛直方向の各レベルについて、水平2次元の空間分布を見ることが多い。 しかし、場所を固定して、数値の時系列として見たいばあいおある。 (数字で書きだすのならば、作図の部分は不要になる。) そのようなプログラム例は、まだ汎用の教材としてつくっていないが、 つくりかけの例をしめすことにする。
つかったデータ例は、JRA-55 の月平均値 (Monthly) の 地表面解析値の 緯度経度 1.25度格子 (surf125) データ (GRIB2形式で提供されている) のうち 2010-2020年の気温を、あらかじめ単純バイナリ (4バイト浮動小数点数) 形式に変換しておいたものである。 データは 実習用計算機の /home/masuda/jra55/monthly/tmp/ の下にある。 このウェブページにはサンプルとして 1個だけデータファイルをつけておく。
また、同じ変数の 1991-2020年の平年値 (clim9120) もつかっている。 実習用計算機の /home/masuda/jra55/clim9120monthly/tmp/ の下に置いた。 【clim9120というディレクトリの下にもあるのだが、ちがうファイル編成になっている。】
プログラム例の datapath のところは (実習用計算機ではなく) 増田のパソコンのフォルダ名のままになっている。 各自、データを置いたところにあわせて書きなおしてほしい。
もし Daily (実際は6時間間隔) のデータについて同様なことをするならば、 year, month のループの内側に day, hour のループ があるような構造にすればよい。
JRA-55の緯度経度 1.25度格子とはちがう格子のデータをあつかうときに書きかえる必要があるところは、 データを読むところ (np.fromfile() の呼び出しをふくむ一連の作業)、 lon, lat の配列に値をあたえる計算式、 それに (a) では losel, lasel に値をあたえる計算式である。
プログラムのはじめのほうで、経度、緯度が lonsel, latsel の値として指定されている。 (単位は「度」、東経・北緯が正である。) 各年各月のデータを読み、 print文で、(年、月、格子点の緯度、経度のあとに) 指定された位置に近い格子点 (端数の処理は 関数 int にまかせた) の値を 標準出力に書きだしている。 (「python timeseries_tmp_onepoint.py > tmp.txt」のように標準出力をファイルに redirect するとよい。)
経度と緯度の範囲が、lonmin, lonmax, latmin, latmax として指定されている。 各年各月のデータを読み、 print文で、(年、月、latmin, latmax, lonmin, lonmax のあとに) その時刻の、指定範囲の空間的な平均値、最小値、最大値を標準出力に書きだしている。
平均値の計算は、ライブラリ関数をつかわず、 格子点ひとつひとつの値をたしこむ方法で書いてみた。 全球の格子点について走査して指定範囲にあるかどうか判断するという、能率面からは最適でないプログラムになっている。 【緯度・経度の範囲から、2次元配列の行番号・列番号の範囲をあらかじめもとめておいて、 その範囲を走査するほうが、よいプログラムだろう。】 欠損の判断は np.isnan() によっている。
緯度経度の格子では、ひとつの格子点が代表する面積が緯度によってちがう。 まず、この面積のちがいを無視して単純に平均するプログラムをつくった (コメントとしてのこしてある)が、 いまの版では、ます目の中央の緯度の cos を重みとするようにしてある。
データセットの気温の単位は K だが、出力は ℃ にした。換算は結果を書きだす直前に入れた。
各年各月の値のデータのほかに、平年値のその月のデータを読んで、偏差を計算する。 偏差についても (c) と同様に、領域平均を計算して書きだす。 (ちがう配列変数について同じ平均操作をおこなうので、関数 a_average() としてくくりだした。)
【[2022-10-22 例文訂正] これまでリンクしていたプログラム例 timeseries_tmp_anomaly_box_average.py に まちがいがあったので、訂正したものにさしかえました。 関数 a_average の最後が「return sumval」となっていましたが、正しくは「return aval」です。】
(b) にふくめた最大値・最小値は (c) にはふくめていないが、追加することはむずかしくない。
JRA-55 の緯度経度1.25度格子でそれぞれの格子点が海としてあつかわれたか陸としてあつかわれたかの情報は、 Const というデータグループの LL125.grib というファイルの LANDsfc という変数にある。
ここではそのデータを単純バイナリファイルに変換してつぎのところに置いた。 南北 145、東西 288 の格子で、値が 0 のところが海、 1のところが陸である。
基本的に (b) と同じだが、あらかじめ海陸分布のデータを読みこんでおき、海の格子点についてだけ集計するようにした。 コメントとしてある文にさしかえると、陸の格子点についてだけ集計するものにもできる。
(b) からはもうひとつ変更した。 (b) のプログラムでは、指定した緯度または経度の境界が格子線に一致しているとき、境界にあたるます目の値も同じ重みで計算する。 (d) のプログラムでは、緯度方向、経度方向それぞれ、境界にあたるます目の重みを半分にした。 その線を境界とする両方の領域の集計をさらに集計するとき、境界上の値が二重になるのはうまくないと考えたからである。
(b) の例文は変更していない。もし海陸の区別はしないが境界のます目の重みを半分にしたいばあいは、 (d)のプログラムでコメントにしてある「if True:」を生かすか、(d)の修正点にあわせて(b)を修正してつかってほしい。