環境気象学アドバンスト 2023年度 増田担当分
増田担当分では、計算機をつかって、気象データを図化してみます。
前期のセミナーで学ぶ Python プログラムを前提としますが、
今年度セミナーをとらない受講生には必要に応じて個別指導します。
2023年5月10日のうちあわせの結果、
気象観測をふくむ日帰り巡検 (10月14日) のあと、
その観測で得られたデータやその周辺の情報を作図してみることになりました。
周辺の情報については、つぎのもののうちから、観測のねらいと、当日の天気とを考えて、
いくつかえらびます。
- 地上気象観測値 (気象官署 および アメダス) 気温、湿度、日照時間、風など
- 高層気象観測値 (館野のラジオゾンデ)
- 全国合成レーダー ... 空間分解能 約 1 km の雨雲の情報
- 気象衛星ひまわり赤外画像 ... 空間分解能 約 2 km の 雲頂温度 または 地表面温度に関連する情報
- 気象衛星ひまわり可視画像 ... 空間分解能 約 0.5 km の 雲の情報
増田担当分 各回の内容
- 10月20日 アメダスの風の時別値を地図に表示するプログラムをうごかしてみる。
- xming, teraterm を起動して、hamana にログインする。
- cp -r /home/masudako/ex/mapamedas . (空白のあとピリオド)
- cd mapamedas
- emacs -nw map_amedas_wind_hourly.py (control-x control-s で保存、control-x control-c で終了)
- python3 map_amedas_wind_hourly.py
- cd image
- display {pngファイル名} または animate -delay 100 *.png
- cd ..
- 10月27日 アメダスデータの作図 (風)
- つぎのプログラムの機能をたしかめ、8月13・14日のデータを作図してみよう。
- [気象庁アメダスの風のデータ (時別値) を地図表示する]
- map_amedas_wind_hourly.py
- graph_amedas_hourly_wind_ver1.py ← 1地点の風向・風速の時間変化のグラフ
- 変更するところ
- 入力ファイル名の日時 syrmo = '20231013_14'
- graph では、地点名 stnamesel に代入する文字列 (ローマ字、頭が大文字)
- 都道府県名。
- map では regionlist の内容を、関東 7都県 にする (今回用意したデータはこの範囲なので)。
- graph では region に、選択した地点がふくまれる都道府県名 (こちらは全部小文字)。
- map では地図の緯度経度範囲 latmin, latmax, lonmin, lonmax
- map wind では 凡例の矢印をかく位置 (地図上の緯度経度) lata[0], lona[0]
- (今回の対象ではないが map contour では格子点の緯度経度範囲 latming, latmaxg, lonming, lonmaxg)
- graph では、横軸の目盛りの日時の書式 formatter の "%d %H:%M" を "%H" に
- graph では、横軸のラベル "June 2023" を "13-14 Oct. 2023" に
- 結果を画面に表示する plt.show() と、pngファイルに保存する plt.savefig をふくむ5行との選択
- 気象庁サイトからデータをダウンロードする手順
- 11月10日 アメダスデータの作図 (気温など)
- 11月17日 アメダスデータと 自分たちの観測データを同じ図に表示してみる (時系列グラフ)
- hamana で
cp -r /home/masudako/ex/20231117 .
cd 20231117
をして作業してほしい。
【[2023-12-01] このもとのディレクトリは /home/masudako/adv2023/20231117 に移動する。】
- 前回 Excel から書きだした 4つの csv ファイルを、ひとつにまとめた。
今後、ruobs.csv をつかう。
- 変換につかったプログラムは、edit_obs.py である。
これは同じ処理を 6回くりかえし書いており、プログラムの手本としてはうまくない。
【11月17日夜、Python の「関数」の機能をつかって同じ処理を1回書く形に整理した。
edit_obs_v2.py として /home/masudako/ex/20231117 に追加した。
関数のなかで書きかえられて関数の外でつかわれる変数 dfwrite を
global 宣言しておくことが重要である。】
- シフトJIS文字コードのエラーが生じたので、codecs をつかって回避した。
([トラブルシューティング事例集]のページ参照。)
- location という列に、観測場所 (霞が関のばあいは場所と観測班の略号) を入れた。
【loc と略そうと思いかけたのだが、loc は Pandas のデータフレームでは
きまった意味をもつので、列名につかってはまずいのだった。】
- 風速が 0 のところは、風向はデータ欠損 (CSV では空白もない空文字列) にそろえた。
- 気象庁の観測点「東京」と「熊谷」の 10月13-14日の時別値をダウンロードしなおした。
気温、降水量、風向・風速だけでなく、相対湿度、露点温度、海面気圧、現地気圧、日照時間、日射量などをふくめた。
(熊谷では日射量を観測していないので、その列は全部欠損値となる。)
そのデータは hamana では /home/masuda/advanced2023/ に置いた。
(WinSCP でパソコンに転送して Excel でみるとよい。)
- graph_amedas_hourly_wind_tokyo.py, graph_amedas_hourly_t_pr_tokyo.py は
前回つくったのと同じ 2日間の時系列グラフを、今回ダウンロードしたファイルでつくるものである。
地点表との照合の必要がないぶんだけ単純なプログラムになった。
- graph1.py も 2日間の時系列グラフで、
気温と同じわくに露点温度、その下に、相対湿度、日照時間、日射量を表示する。
日射量はCSVファイルでは 1時間の累積値 (単位 MJ/m2)ではいっているが、
単位時間あたりの値 (単位 W/m2) に換算する。
1か所「Tokyo」を「Kumagaya] に変えれば熊谷の観測値を表示するようにつくってある。
画像をファイルに保存するには、plt.show を plt.savefig をふくむ5行に変えればよい。
- graph2.py も 同じ量のグラフだが、日時の範囲を 14日の昼間にかぎったものである。
- [注意] 日照時間、日射量、(この図にはないが) 降水量の、
たとえば、13時のところに表示されている値は、12時から13時までの累積量である。
12時から13時まで横棒をひくように変えたかったのだが、11月17日の授業にはまにあわなかった。
【11月17日の夜に、全期間まとめてではなく1時間ごとに plot をよびだす方法で実現できた。
そのプログラム例を graph2a.py として /home/masudako/ex/20231117 に追加した。
graph3.py も同様に修正可能だが、まだつくっていない。】
- graph3.py では、graph2.py に、
自分たちの観測の値 (ruobs.csv から読む) をかさねて表示してみる。
dfruobs.location[krecru] の値 (文字列) によって分岐して、
marker, markersize, color (の いずれか) を変えて観測場所を区別できるようにしよう。
if dfruobs.location[krecru] == 'shinjuku':
....
elif dfruobs.location[krecru] == 'kanda':
....
else
....
- 横軸 (日時) の範囲 xlim0, xlim1 を、データを見やすいように調節しよう。
- 相対湿度 (dfruobs の 列名 rh) も表示してみよう。
- (つぎの部分は17日の授業では省略した。)
- 露点温度 (dfruobs の 列名 tdew) も表示してみよう。
- 必要に応じて、縦軸 (気温など) の範囲 ylim0, ylim1 を調節しよう。
(set_yticks で指定している目盛りの数値も調節する。)
- 図示する気象要素を変更して、気圧をふくめてみよう。
(気象庁観測値CSVで現地気圧の列番号は 20、dfruobs で気圧の列名は pres)
- 熊谷についてもやってみよう。(プログラムのコピー graph3kumagaya.py をつくってから編集するとよい。)
- 11月24日 アメダスデータと 自分たちの観測データを同じ図に表示してみる (時系列グラフのつづき、地図)
- hamana で
cp -r /home/masudako/ex/20231124 .
cd 20231124
をして作業してほしい。
【[2023-12-01] このもとのディレクトリは /home/masudako/adv2023/20231124 に移動する。】
前回の「20231117」のディレクトリを参照するので、それと横ならびの位置につくってほしい。
- graphpres.py ... 時系列グラフ。
前回の graph3.py の日射量のかわりに気圧をふくめた。
気象庁の観測については、東京の現地気圧を緑、海面気圧を青で表示した。
立正大学の観測については、熊谷以外の全地点を区別なく赤で表示した。
(結果をみると、残念ながら、センサー較正がうまくいっていなかったようだ。)
前回に気温についてやったのと同様に地点を区別してみてほしい。
また、熊谷についても作図してみてほしい。
- graphwind.py ... 時系列グラフ。
前回の graph_amedas_hourly_wind_tokyo.py の時間範囲を graph2.py と同様にしたうえで、
立正大学の観測をくわえた (熊谷以外の全地点を区別なく赤で表示した)。
ただし、地上での観測は、気象庁の(原則) 地上 10 m での観測よりも風速が弱いので、
factor (暫定的に 2.0 としてある) をかけて表示することにした。
前回に気温についてやったのと同様に地点を区別してみてほしい。
また、熊谷についても作図してみてほしい。
- ruloc.csv ... 立正大学の観測の地点の緯度・経度の表
(緯度・経度は「地理院地図」からおよその位置の値を読んだ。)
- mapwind_amedas.py ... 地図。
前々回の map_amedas_wind_hourly.py と同様に気象庁アメダスの風を矢印で表示する。
ただし、緯度経度範囲と期間をかぎった。
- mapwind.py ... 地図。
mapwind_amedas.py に 立正大学の観測をくわえた。
ただし、観測時刻は 1時間ごとにまとめた。
地上での観測には factor (暫定的に 5.0 としてある) をかけて表示することにした。
- maptmp_amedas.py ... 地図。
前々回の map_contour_amedas_hourly.py と同様に気象庁アメダスの気温を数字で表示する。
ただし、緯度経度範囲と期間をかぎった。
- maptmp.py ... 地図。
maptmp_amedas.py に 立正大学の観測をくわえた。
ただし、観測時刻は 1時間ごとにまとめた。
気象庁の観測値を緑、立正大学の観測値を赤で表示するようにしてみた。
ただし、このままでは、数字どうしがかさなって、ほとんど読めない。
どうしたら数値がわかるようにできるだろうか。
- 12月 1日 気象衛星ひまわりの赤外の輝度温度と可視の反射の分布を作図してみる
- 気象衛星ひまわりのデータ
- ある波長帯の放射輝度 (単位面積・単位時間あたりの電磁波のエネルギーの流れ)
- 赤外は、「輝度温度」または「相当黒体温度」 (Tbb) ... 同じ輝度値をもたらす黒体放射の温度。
(大気成分による吸収・射出のすくない 10～12 μm 波長帯では、雲頂温度または地表面温度に近い。)
- 可視は、仮称「反射の強さ」。太陽放射が真上からあたって完全に反射するばあいを 1.0 (100 %) とした相対値。
(衛星観測技術の専門家は「アルベド」というが、気象学でいうアルベドと同じではない。)
- 実際には整数でつめこまれている。整数値から物理量 (輝度温度あるいは反射の強さ) への換算表が用意されている。
- 千葉大学 環境リモートセンシング研究センターで、
データを受信し、緯度経度格子への位置合わせをしている。
ひまわり 9号のばあい、緯度経度範囲は、南緯60度～北緯60度、東経85度～205度 (西経155度)。
格子間隔は、可視は緯度経度 0.005 度 (約 500 m)、赤外は0.02 度 (約 2 km)。
- 増田が千葉大から ftp でデータを hamana に持ってきて、
整数値から物理量 (浮動小数点型) への変換と、領域の切り出しをした。
切り出された領域は、可視は 北緯30～50度、東経125～145度 (符丁 sjapan)。
赤外は 北緯20～50度、東経120～150度 (符丁 gjapan)。
- 作図例
北緯 34.5～38.5度、東経 136.5～141.5度 (関東・中部地方、符丁 gkanto) のひろがりで作図した。
階級わけと色は、暫定的に、大きい値を赤、小さい値を青、中間をうすい黄色としてみた。
(気象衛星の習慣は、可視は大きい値ほど白く、赤外は小さい値ほど白くする。)
- 実習作業
hamana で
cp -r /home/masudako/adv2023/himawari .
cd himawari
をして作業してほしい。
- map_tbb.py ... 1時間間隔の赤外画像(輝度温度)を作図。結果は image_tbb の中
- map_refl.py ... 1時間間隔の可視画像(反射の強さ)を作図。結果は image_refl の中
- 【[注] このプログラムをそのまま実行すると、
image_tbb, image_refl の下に「年・月」のディレクトリ「202310」をつくり、
その下に画像ファイルをつくります。
「# 図を画像ファイル (png) に保存」につづく imagefilepathx をつかう 3つの行をコメントにし、
imagefilepathname をくみたてるところの imagefilepathx を x のない imagefilepath に変えれば、
「年・月」のディレクトリをつくらなくなります。】
- 日時の範囲を変えるには daystart, hourstart, dayend, hourend にあたえる値を変える。
UTC (世界時) で指定している。
UTC の 0時は 日本標準時 (JST) の 9時。
- 緯度経度範囲を変えるには lonmin, lonmax, latmin, latmax にあたえる値を変える。
緯線・経線の間隔 dlat, dlon も変えたほうがよいかもしれない。
- 階級わけは levels[ ]、それに対応する色は colorf[ ] で指定している。
色の指定のしかたは、文字列として、'#' につづいて、R (赤), G (緑), B (青) 値をそれぞれ
2けたの16進数で指定している。黒が '#000000', 白が '#FFFFFF'。
グレースケールならば、十進小数を引用符でくくった文字列
(たとえば '0.00', '0.05', ... '0.95', '1.00' など) をつかうこともできる。
黒が 0、白が 1。
- この日に見られた 波状? ロール状? の雲をおいかけることは、
1時間間隔ではむずかしいが、10分間隔ではできるようだ。
- map_tbb_10min.py ... 10分間隔の赤外画像(輝度温度)を作図。結果は image_tbb_10min の中
- map_refl_10min.py ... 10分間隔の可視画像(反射の強さ)を作図。結果は image_refl_10min の中
- 【[注] 12月1日授業開始時のバージョンは、画像ファイルを保存するところでエラーになりました。
「# 図を画像ファイル (png) に保存」につづく imagefilepathx をつかう 3つの行をコメントにし、
imagefilepathname をくみたてるところの imagefilepathx を x のない imagefilepath に変えれば動きます。
/home/masudako/adv2023/himawari にあるものは授業終了後に修正しました。】
- 12月 8日 気象衛星ひまわりの赤外の輝度温度と可視の反射の分布を作図してみる (つづき):
- 10分間の変化 と もとの量 (輝度温度、反射の強さ) とを同じ大きさの地図でならべてみよう。
雲が移動したのか、あらたに出現したのか、判断できそうだ。
- hamana で
- cd himawari
- cp /home/masudako/adv2023/himawari/map_delta_refl_10min.py .
- cp /home/masudako/adv2023/himawari/map_delta_tbb_10min.py .
- mkdir image_delta_refl_10min
- mkdir image_delta_tbb_10min
- map_delta_refl_10min.py, map_delta_tbb_10min.py をそれぞれ編集して実行
- 変えるところ
- 緯度経度範囲 lonmin, lonmax, latmin, latmax
- 時刻の範囲 hourstart, hourend
- 12月 15日 気象衛星ひまわりの赤外の輝度温度と可視の反射の分布を作図してみる (3):
- 前回・前々回のプログラムの訂正
- 「xticks = np.arange(lonmin, lonmax+lolaeps, dlat)」
の dlat とあるところは、
経線の間隔なので dlon が正しい。
(dlon と dlat の値をちがえたばあいに修正が必要。同じならばさしつかえない。)
- 局所的極大・極小を検出・表示する
- ここで局所的極大とは、緯線番号 j、経線番号 i の格子点の値が、
緯線番号 j ± jdy, 経線番号 i ± jdx の長方形にふくまれる格子点のうちで最大であるばあいとする。
たとえば jdx = 5, jdy = 5 ならば、11 × 11 = 121 個の格子点のうちの最大である。
- プログラム例は /home/masudako/adv2023/himawari/ の下に用意した。
- 地図表示とテキストファイル書きだしを同じプログラムでやったほうが総作業時間は短くなるのだが、
プログラムがすこしわかりにくくなるので、ここでは別々にした。
- 実行前に、プログラムの jdx, jdy と
緯度経度範囲 (lonmin, lonmax, latmin, latmax)、
時間の範囲 (hourstart, hourend ... UTC で指定) を確認し、変えたければ変えること。
また、出力データが置かれるディレクトリをつくっておくこと。
- map_refl_10min_localmax.py ... 可視の反射の強さの極大を地図上に示す。
画像は image_refl_10min のディレクトリ (前々回の10分値と同じ) につくられる。
- refl_10min_localmax.py ... 可視の反射の強さの極大の経度、緯度、反射の値をテキストファイルに書きだす。
出力データは txt_refl_10min のディレクトリにつくられる。
- map_tbb_10min_localmin.py ... 赤外の輝度温度の極小を地図上に示す。
画像は image_tbb_10min のディレクトリ (前々回の10分値と同じ) につくられる。
- tbb_10min_localmin.py ... 赤外の輝度温度の極小の経度、緯度、反射の値をテキストファイルに書きだす。
出力データは txt_tbb_10min のディレクトリにつくられる。
- 東西断面、南北断面の値を抽出する
- プログラム例は /home/masudako/adv2023/himawari/ の下に用意した。
- いずれもテキストファイルを書きだすプログラムであり、作図はふくんでいない。
[2023-12-18] Pythonで東西断面のグラフをつくるプログラム例を追加した。
下の [2023-12-18 補足] のところに説明を書いた。
- refl_10min_ew.py ... 可視の反射の強さの、latsel で指定した緯度にいちばん近い格子緯度の、
lonmin から lonmax までの経度での値を、経度・緯度とともに書きだす。
出力データは txt_refl_10min のディレクトリにつくられる。
- refl_10min_ns.py ... 可視の反射の強さの、lonsel で指定した経度にいちばん近い格子緯度の、
latmin から latmax までの緯度での値を、経度・緯度とともに書きだす。
出力データは txt_refl_10min のディレクトリにつくられる。
- tbb_10min_ew.py ... 赤外の輝度温度の、latsel で指定した緯度にいちばん近い格子緯度の、
lonmin から lonmax までの経度での値を、経度・緯度とともに書きだす。
出力データは txt_tbb_10min のディレクトリにつくられる。
- tbb_10min_ns.py ... 赤外の輝度温度の、lonsel で指定した経度にいちばん近い格子緯度の、
latmin から latmax までの緯度での値を、経度・緯度とともに書きだす。
出力データは txt_tbb_10min のディレクトリにつくられる。
- 出力データのファイル名には、ew は日付と緯度、ns は日付と経度がふくまれる。
ただし、緯度 (または経度) は latsel (またはlonsel) で指定した値そのままではなく、
実際に選択された格子線の緯度 (または経度) を、小数点下 4 けたでふくむようにした。
- 南北にのびた雲が東西にならんでいたので、
各時刻の東西断面をとり、横軸をそろえ、縦軸を時刻ごとにずらして
(Python の matplotlib ならば、横長の subplot を縦にならべて) グラフにすると、
雲が形をたもちながら移動しているかどうかがわかると思う。
- 雲があるばあい、輝度温度は雲頂温度に近いので、
温度の鉛直分布がわかっていれば、雲頂高度を推定できる。
- [参考] ラジオゾンデによる気象要素の鉛直分布
- [高層気象観測データの鉛直プロファイルのグラフ (Pythonによる試作)]
- ↑のページで説明したプログラムと、
館野 (Tateno、つくば市) での 2023年10月14日 日本時間 9時と21時
(UTC では 0時と12時) の「指定気圧面」の観測データを
hamana の /home/masudako/ex/vertical/ の下に置いた。
ディレクトリごとコピーして、plotz4_temp_tdew.py などをためしてみてほしい。
- [2023-12-18 補足]
- 東西断面のグラフを作図するプログラム例を /home/masudako/adv2023/himawari/ の下に用意した。
- graph_refl_10min_ew.py
- refl_10min_ew.py の結果を読み、
可視の反射の強さの東西断面をグラフにする。
- 経度範囲を lonmin, lonmax (浮動小数点数) で指定する
(refl_10min_ew.py で指定した範囲またはその部分)。
- 緯度を slat (文字列) で、txt_refl_10min のディレクトリに置かれた
refl*utc_lat*.txt の形のファイル名の lat につづく文字列のとおりに指定する。
(たとえば '36.1475')。
- 開始時刻を UTC の「時」単位で hourstart (整数) で指定する。
- 現状のプログラムでは、画面を6段にわけて、hourstart が 5 ならば
上から順に 05:00 UTC から 05:50 UTC までの 6つのグラフをかく。
- 結果の画像ファイルが image_graph_refl_10min のディレクトリにつくられる。
(ディレクトリがまだなければ、python3 実行のまえにつくっておくこと。)
- graph_tbb_10min_ew.py
- tbb_10min_ew.py の結果を読み、
赤外の輝度温度の東西断面をグラフにする。
- 経度範囲を lonmin, lonmax (浮動小数点数) で指定する
(tbb_10min_ew.py で指定した範囲またはその部分)。
- 緯度を slat (文字列) で、txt_tbb_10min のディレクトリに置かれた
tbb*utc_lat*.txt の形のファイル名の lat につづく文字列のとおりに指定する。
(たとえば '36.1300')。
- 開始時刻を UTC の「時」単位で hourstart (整数) で指定する。
- 現状のプログラムでは、画面を6段にわけて、hourstart が 5 ならば
上から順に 05:00 UTC から 05:50 UTC までの 6つのグラフをかく。
- 結果の画像ファイルが image_graph_tbb_10min のディレクトリにつくられる。
(ディレクトリがまだなければ、python3 実行のまえにつくっておくこと。)
- (南北断面のグラフを作成するプログラムはまだつくっていないが、
東西断面のグラフのプログラムを改造してつくることができる。やりたい人は質問してほしい。)
- 12月 22日
- HTMLの表[ひょう]に画像をならべてブラウザで見る。
- これまでやったことをふりかえる。
- 1月 5日 この授業はありません。
- 年末年始中、増田はメールでの質問に (即日とはかぎりませんが) 応じ、
zoom による会話で相談にのることができます。
Teams のアドバンストのチームはファイルを置くのにつかうことができますが、増田は常にはチェックしていません。
- 1月 12日 各自 (10分ぐらい) 発表、議論
- 10月14日の東京・熊谷付近の気象について、
あるいは、気象観測、気象データ処理技術について、
わかったこと、疑問がのこったことなど
10月13日 (金) 夕方、わくをつくりました。
15日 (日) に、天気図、気象衛星ひまわり9号の画像、地上気象 (アメダス) の観測値データ を入れました。
2023年度「セミナーの基礎」の観測に関連して収集した
[2023年度 観測実習期間とその前後の気象庁による観測値・天気図]
と、まったく同じではありませんが同様の内容です。
17日 (火)、14日午後のひまわりの10分間隔の画像を追加しました。
計算環境 Linux
増田研究室の Linux マシン "hamana" に ログインして実行する。
つかいかたはだいたい "luce" と同じにした。
Windows のばあい、xming を起動したうえで、TeraTerm の ssh でログインする。
ファイル転送は WinSCP で可能。
IPアドレスが luce と (にているが) 最後の部分がちがい、大学ネット内と外でもちがうことに注意。
emacs の操作について
- 実習では原則として 「emacs -nw」 (no window) をつかってください。
- キー操作の表記
- c-x は control (ctrl) キーをおしながら x をおす。
- m-x は (ふつうのキーボードには meta キー がないので) escape (esc) キーをおしたあとで x をおす。
- 保存 c-s
- 終了 c-c
- (変更して保存しないで終了しようとすると質問が出る。
保存しなくてよいとき「n」でこたえるべき問いと「yes」でこたえるべき問いがある。)
- 入力途中のコマンドをとりけす c-g
- 検索 c-s
- cut and paste, copy and paste
- 対象のはじめの位置で c-space または c-@
- 対象のおわりの位置で
- [cut ならば] c-w
- [copy ならば] m-w (ESC w)
- はりつけるところで c-y
- [外] 慶応大学 湘南藤沢キャンパス SFC-CNSガイド 2001年版より [Emacs]
- [外] unixuser.org Yusuke Shinyama さん (2001年) [Emacs クイックリファレンス]
- [外] Emacs JP [Emacsチュートリアル 日本語訳] (GNU emacs 本家のドキュメントの日本語訳)
- [外] Emacs JP conao3 さん [2020年代の Emacs 入門] (やや専門的か?)
- わたしのLinux上での実習では日本語文字 (漢字・かな) をつかわないので、他のサイトの教材の 日本語入力関係の話題は ひとまず無視してください。
大気汚染観測点の情報
気象の参考書
- 藤部 文昭, 2012: 都市の気候変動と異常気象。朝倉書店。
[読書メモ]
- 浅井 冨雄, 1996: ローカル気象学。東京大学出版会。
(更新) 2023-12-22 朝
増田 耕一 (MASUDA Kooiti)