観測点ごとの気象要素の表示と、格子点データによる等値線などの表示を、同じ地図上にかさねてみる。 ここでは、海面気圧について、つぎの 3つのデータによって、それぞれ等圧線をひいてみる。
それと、観測点の気圧を数字で表示したものをかさねてみたり、 風ベクトルを矢印で表示したものをかさねてみたりする。
気圧以外の気象要素についても、同様にあつかうことができる。 海面気圧をあつかうために JRA-55 では psl、MSM では psea としてあるところを、 気温ならば、JRA-55 では tmp、MSM では temp とすればよい。 (JRA-55 のファイル名を決める際に MSM のことを考慮しなかったので不統一になってしまった。) ただし、物理量の単位の問題がある。 JRA-55 や MSM の 気圧は Pa であたえられて hPa で表示したいので、読んだらすぐ 100 でわっている。 気温は K であたえられているので、℃ で表示したいならば、273.15 をひく。 また、とりうる数値の範囲を考えて、等値線レベルをあたえている levels0, levels1 のリストの内容を変更する。
また、MSM の地上気温の分布は、地面の標高の影響ばかりが見えてしまう。 地方スケールの気温の水平分布を知りたいならば、 MSM の気圧面の 3次元の場のデータ (MSM-P) で (たとえば) 850 hPa 面の気圧を見たほうがよさそうだ。 MSM-P のつかいかたの説明は今回の授業ではしめさないが、夏のうちには教材ウェブページに追加したい。
水蒸気については、JRA-55 では比湿 (q)、MSM-S では相対湿度 (rh) がある。
今回のPythonプログラムの例文は、luce の /home/masuda/mapstationgrid に置いてある。
各自、つぎのコマンドによって、プログラム例をコピーしてほしい (最後の空白のあとのピリオドが重要)。
cp -r /home/masuda/mapstationgrid .
そして、mapstationgrid に cd して作業してほしい。
今回つかうデータのうち、気象庁の地上気象観測データは、前回の [観測点ごとの風のデータを図示する]でつかったものと同じなので、 各自の mapstationwind のディレクトリにもあるはずだが、 今回は /home/masuda/mapstationwind にあるものを読むことにする。 (別のところにあるデータを読みたい人は、プログラムの変数 datapath にあたえる文字列を変更してほしい。)
また、JRA-55 のデータ (2023-06-13 の [格子点型の気象データ (JRA-55 海面気圧) の地図上表示]でつかったもの) と、 MSM のデータ (下で説明する) も、/home/masuda/ の下にあるものを読む。
観測点ごとの数値から、等値線 (気圧ならば等圧線、温度ならば等温線) をひきたいことがある。 手がきならば、 近くにある 2つの観測点をむすぶ線分を考えて、そのうちどこに等値線をひくべき値がくるかを判断する ことをくりかえしていくのがふつうの方法だ。 計算機で等値線をひく方法のうちにも、そのような方法がある。 (x,y)平面、あるいは地球の球面を、観測点によって三角形網に分割して、それぞれの三角形の辺上で 等値線がとおる位置をもとめていく。 【わたしは、そのような作業をする Fortran言語によるプログラムをもっているが、まだ Python になおすことができていない。】
それよりもよくつかわれるのは、観測点の位置からの空間内挿によって格子点データをつくり、 格子点データにもとづいて等値線をひく技術を応用することである。 Python の SciPy という大きなパッケージのうちに interpolate.griddata という関数があるので、 ここではそれをつかう。 この interpolate では内挿方法として、nearest, linear, cubic を選択できるが、 ここでは linear (空間座標の1次関数による内挿) を採用する。 【nearest は、内挿というよりも入力値のうちで位置がいちばん近い値をえらぶことで、 データが数量ではなく番号のようなものである場合に適する。 cubic は 3次関数で、入力データがなめらかな関数値である場合に適する。】
scipy.interpolate.griddataをつかう準備として、「from scipy import interpolate」を指定している。
前回と同様に、観測地点表のデータフレーム dfstn をつくり、その列 col00 に、つぎの dfobs での列番号を書きこんでおく。 そして、2023年6月1-6日の毎時の「気象官署等」の海面気圧のデータを、データフレーム dfobs に読みこむ。
時刻は、6月1日 1時 (日本時間) を 0 として、krecstart から krecend まで 1時間ずつ進む。 (例文では、krecstart = 0, krecend = 71 つまり 6月4日 0時までとしてある。)
内挿の準備として、配列を確保する。
各時刻について、観測点の緯度を xp, 経度を yp, 気圧の値を zp につめかえ、 (xp, yp) から zp の値を (xg, yg) に内挿した結果を zg とする。
地図をかき、前回の [31a_map_hourly.py] と同様に、 観測点の位置に 〇印をかき、その上側に観測値を数字でかく。
contour によって zg の等値線をひく。 ここでは、1000 hPa をふくむ 4 hPa おきのレベルに実線、 それから 2 hPa ずれた 4 hPa おきのレベルに破線で等値線をひくことを、 contour を 2回よぶことで実現した。 (くふうすれば 1回ですませることも可能だが、このほうがわかりやすいと判断した。)
結果は、作業ディレクトリの下のディレクトリ image の下に png 画像ファイルができる。 「display *.png」または 「animate -delay 100 *.png」などで見てみよう。
対象地域を変えるには、つぎのところを変更する。
前回の [30c_map_winduv.py] と同様に、 各観測点の風ベクトルを矢印で表示し (数字による表示は省略)、 それに、今回の [33a_map_interp_hourly.py] と 同様に、観測点の海面気圧からの空間内挿でつくった格子点データによって等圧線をひく。
【風のデータと海面気圧のデータの形式は、だいたい同じなのだが、すこしちがう。 いまのプログラムでは、それぞれを読むところを別々に書いており、プログラムが長くなっている。 ひとつの関数にまとめてそれを引き数の値を変えて2回呼び出す形になおしたほうがよいと思うが、 わたしが Python の変数ひきわたしをまだよく理解していないので、できていない。】
再解析データとは、気象の数値予報モデルに観測データをとりこむデータ同化でつくられたデータのうち、 リアルタイムではなく、一定の方式で長期間にわたって作成されたものである。 ここでは、2023-06-13 の [格子点型の気象データ (JRA-55 海面気圧) の地図上表示] と、 2023-06-20 の [格子点型の気象データによる風の矢印 (および、矢羽根) 表示] でつかった、 JRA-55 再解析の、daily (正確には6時間ごと)、surf125 (地表面付近の2次元データで、緯度・経度1.25度格子にしたもの) をつかう。 2021年1月1日から2023年6月10日までの 増田が単純バイナリに変換したデータが /home/masuda/jra55/daily/surf125/ の下に置かれている。
【JRA-55 の数値予報モデルの格子は 1.25度よりもこまかく、その格子でのデータも公開されているのだが、 東西ひとまわりの格子数が緯度によってちがうので、わたしはまだつかいかたがわかっていない。】
JRA-55 のデータは、毎日、世界時 (UTC) の 0, 6, 12, 18 時 (日本標準時 (JST) では 9, 15, 21, 翌日 3時) に値があたえられる。 プログラム例では、地上観測CSVデータにあわせて JST で日時をかぞえ、 datetime パッケージをつかって UTC に換算して JRA-55 のファイル名を得るようにした。 プログラムの krecstart は、このデータのある時刻に対応させておく必要がある。 そこから krecend まで、6時間おきに作図をつづけるようにしてある。
データ同化のプロダクトには、再解析のほかに、リアルタイムの数値予報のためにつくられるものがある。 気象庁では、全球モデル (GSM) と、領域をかぎって格子間隔をこまかくした メソモデル (MSM) をうごかしている。
気象庁の数値予報の出力データを、リアルタイムに得るには、 気象業務支援センターと有料の契約をする必要があるが、過去のデータは再配布可能になっている。 気象庁でつくられた格子点データの形式は WMO (世界気象機関) できめられた GRIB version 2 だが、 京都大学 生存圏研究所では、受信した MSM のデータを、 アメリカの Unidata できめられた NetCDF という形式に変換して、 「生存圏データベース」http://database.rish.kyoto-u.ac.jp のつぎのページから リンクされたところに置いている。
MSM-S と MSM-P の2つのグループのデータがある。
MSM-Sデータは、(おもに地表面の) 水平2次元データで、つぎの変数がふくまれる。 psea (海面気圧)、sp (モデル地表面気圧)、u (10 m 風速の西風成分)、v (10 m 風速の南風成分)、 temp (2 m 気温)、rh (2 m 相対湿度)、r1h (1時間降水量)、 ncld (雲量)、ncld_upper (上層雲量)、ncld_mid(中層雲量)、ncld_low (下層雲量)。 水平格子間隔は実際は緯度経度だが約 5 km である。 時間間隔は 1 時間だが、そのうち hour が (UTC でも JST でも) 3の倍数の時刻は解析値、 他の時刻は予報値がはいっている。 降水量はすべて予報値である。 観測値の情報がとりこまれるのは解析値がつくられるときである。
MSM-Pデータは気圧面の3次元データで、つぎの変数がふくまれる。 z (ジオポテンシャル高度)、w (鉛直p速度)、u (風速の西風成分)、v (風速の南風成分)、temp (気温)、rh (相対湿度)。 【ただし、w は、増田が 2023-06-30 までに見た 2022年のいくつかの日時の netCDF ファイルに関するかぎり、すべて 0 がはいっている。】 水平格子間隔は実際は緯度経度だが約 10 km (MSM-Sの格子のひとつおき)である。 時間間隔は 3 時間で、(2006年3月以降は) すべて解析値である。
ファイルは1日ごとに複数時刻がまとめられており、「0601.nc」が 6月1日をあらわすような形で名づけられている。 年はファイルが置かれたディレクトリ名で区別する。 (netCDFなのでファイルの内容にも日時の情報はある。)
ここでは、2023年5月31日から6月9日までの MSM-S データを、/home/masuda/msm/s/ の下に置いた。
Python で netCDFを読み書きするためのパッケージとしては、netcdf4 というものがある。 (この 4 は netCDF の version 4 をさしている。) そのほかに、xarray というものがつくられ、これでも netCDF を読むことができる。 ここでは、京大防災研の 榎本 剛 [たけし] 先生によるつぎのページ (MSM-P を読むもの) を参考に、 xarray によって numpy の配列に読みこむことにした。
作図対象の日時は、例文では、 krecstart = 0 つまり 6月1日 日本時間 1時から、 krecend = 71 つまり 6月4日 0時まで、1時間ごととしてある。 datetime パッケージによって 世界時 (UTC) に変換し、読むべき MSMデータを判断している。
MSMの格子間隔がこまかいので、intv = 5 として、東西方向・南北方向それぞれ 5つおきに矢印をかくことにした。