「大気上端」に入射する太陽放射

主題: 地球の公転・自転にともなって、 「大気上端」に入射する太陽放射 (単位時間・単位面積あたりのエネルギーの流れ)は どのように変化するか。

手段: Python (パイソン) をプログラム言語処理系として使い、数値計算と作図のプログラムを動かす。

(1) 日変化

緯度・経度をあたえられた地点での、 1日のうちの各瞬間の太陽の方位・高度角を計算し、 「太陽定数」を定数として、「大気上端」に入射する日射量を計算する。 現在の地球の軌道要素にもとづく「均時差」の計算がはいっている。

• [注] 均時差は、浜島書店 (2013) 『地学図表』では、49ページ [4 A] に出てくる。 ただし単位が、『図表』では時間の「分」、ここでしめす計算では角度のラジアンである。
• 第一学習社 『地学図説』 では、宇宙の章の17節 「天球の座標」の [2] 「均時差」 として出てくる。

計算式は、立正大学 地球環境科学部 (2021年度からは名誉教授) の 中川 清隆 先生 [中川 清隆 研究室]の [太陽方位、高度、大気外日射量の計算]の ページの記述にもとづいた。

増田が、プログラムを、2019年度に R 言語で書きかえ、グラフの作図を追加した。 さらに 2021年度に Python 言語で書きかえた。

プログラム例

• 「均時差」を計算してファイル "eqtime.csv"に書き出す。例文作成作業の途中の動作確認用。
  • test_kinzisa.py (文字コード UTF-8)
  • 出力CSV例 eqtime.csv
• 1つの時刻の数値を計算する。
  • calc1.py (文字コード UTF-8)
• 1時間きざみで24回の数値を計算し、CSVで書きだす。
  • calc_hourly.py (文字コード UTF-8)
  • 出力CSV例 hour_insol_day140.csv
• 1時間きざみのCSVファイルを読み、グラフにする。「#画面表示」と「#画像ファイル」が選択可能。
  • plothourly.py (文字コード UTF-8)
  • 出力画像例 hour_insol_day140.png
• 1時間きざみのCSVファイルを 2つ読み、同じ画面上にちがう記号でグラフにする。
  • plothourly2.py (文字コード UTF-8)
  • 出力画像例 hour_insol_day140_173.png
• 1時間きざみのCSVファイルを 2つ読み、両者の差の値をグラフにする。
  • plothourlydiff.py (文字コード UTF-8)
  • 出力画像例 hour_insol_day173_minus_day140.png

注

• 日付は1月1日を1とした年ごとのとおしの日付を使っている。 「月・日」とのあいだの相互変換は、Python の datetime パッケージを使う方法も、 自まえでプログラムする方法もあるが、例文にはまだくみこんでいない。
• 日本語文字をふくむ Python プログラムの文字コードは、 「UTF-8」が標準になっている。 Windows上のほかのソフトウェア (たとえば Excel) では「Shift_JIS」が標準になっているので、 文字化けのおそれがある。 文字化けしたら文字コード変換をこころみてほしい。
• 図に日本語文字列をふくめるために、plt.text や plt.title などのよびだしで、 「fontname='MS Gothic'」という形で MS Windows にあるフォントを指定している。 他の OS で実行するときは、その OS の日本語フォント名を指定すればよいかもしれないが、 それでうまくいかなければ、 図にふくめる文字列をアルファベットからくみたてられたものに変えてほしい。
• 例文のプログラムを画面に表示させて保存すると、 ブラウザによっては、ファイル名の「拡張子」部分が「.txt」になったり、 文字コードが変わったりする。 プログラムを表示させず、ファイルへのリンクが書かれているページで、 リンクを(たいていのブラウザでは)右クリックして、 (Firefoxのばあいでいえば)「名まえを付けてリンク先を保存」の機能で保存するのがよい。

演習課題

• calc_hourly.py をもとに、 次の変数に代入される数値をひとつずつ、(プラスマイナス10日以上・10度以上の差で)かえて、結果を比較し、結果の特徴を簡単にのべてください。
  • dayofyear 元旦からの通し日数 [日] (例文では5月20日にあたる数値を 31+28+31+30+20 としてあたえている)
  • geolatdeg 緯度 [度、北緯が正] (例文の値は熊谷地方気象台の緯度)
  • geolondeg 経度 [度、東経が正] (例文の値は熊谷地方気象台の経度)
  • 結果をファイルに書き出すばあいは、値を変えたものを区別できるようにファイル名をつけること。 プログラム calc_hourly.py では、 テキストファイルに書き出すところで、dayofyearの値がファイル名にくみこまれるようにしたが、 緯度・経度はくみこんでいない (自動化してうまくいくか不安なので)。 必要に応じて、calc_hourly.py のおわり付近の「fname = ...」の文の 右辺の文字列をエディタ作業で変更してほしい。
  • solconst いわゆる太陽定数 [W/m²] をかえてみてもよいが、定数倍されるだけで、あまりおもしろくない。

結果の例

緯度・経度・日付を指定し、1日ぶんの1時間きざみの値を、CSVファイルに書きだす。

• 緯度・経度は熊谷のもの、日付は5月20日。
  • calc_hourly.py
  • 出力CSV hour_insol_day140.csv
• 日付を変える
  • 熊谷の緯度経度、6月22日 (ほぼ夏至)
    • calc_hourly_day173.py
    • 出力CSV hour_insol_day173.csv
  • 熊谷の緯度経度、12月22日 (ほぼ冬至)
    • calc_hourly_day356.py
    • 出力CSV hour_insol_day356.csv
• 経度を変える
  • 緯度は熊谷と同じ、東経120度、5月20日
    • calc_hourly_lon120.py
    • 出力CSV hour_insol_lon120_day140.csv
  • 緯度は熊谷と同じ、東経150度、5月20日
    • calc_hourly_lon150.py
    • 出力CSV hour_insol_lon150_day140.csv
• 緯度を変える
  • 北緯 23.4度(北回帰線)、東経 135度、6月22日
    • calc_hourly_lat234_day173.py
    • 出力CSV hour_insol_lat234_day173.csv
  • 北緯 66.6度(北極圏の緯線)、東経 135度、6月22日
    • calc_hourly_lat666_day173.py
    • 出力CSV hour_insol_lat666_day173.csv

1時間きざみのCSVファイルを2つ読み、グラフにする。PNG画像ファイルに書きだす。

• 緯度・経度は熊谷のもの、5月20日(140)と6月22日(173)。
  • plothourly2.py
  • 出力画像 hour_insol_day140_173.png
• 日付を(大きく)変える。6月22日(173, ほぼ夏至)と12月22日(356, ほぼ冬至)。緯度・経度は熊谷のもの。
  • plothourly2_day173_356.py
  • 出力画像 hour_insol_day173_356.png
  • 結果の特徴: 振幅が夏至に大きく冬至に小さい。
• 経度を変える。東経120度と東経150度。緯度は熊谷のもの、日付は5月20日(140)。
  • plothourly2_lon120_150.py
  • 出力画像 hour_insol_day120_150.png
  • 結果の特徴: 位相がちがう。東経120度では日本標準時(東経135度の時刻)で13時、東経150度では11時に極大となる。
• 緯度を変える。北緯23.4度と北緯66.6度。経度は東経135度、日付は6月22日 (173, ほぼ夏至)。
  • plothourly2_lat234_666.py
  • 出力画像 hour_insol_lat234_666.png
  • 結果の特徴: 極大値は低緯度のほうが大きいが、日射のある時間は高緯度のほうが長い。