夏至の日なので（だったので）、pandasとかpvlibとかを使って夏至の日を調べてみた！
（座標系や時間系の知識が十分でないため、厳密では無いかもしれません。詳しい方に補足いただけると嬉しいです）
準備
環境は、Apple Silicon上のPython 3.9です。
$ uname -a
Darwin Iharas-Air 20.5.0 Darwin Kernel Version 20.5.0: Sat May 8 05:10:31 PDT 2021; root:xnu-7195.121.3~9/RELEASE_ARM64_T8101 arm64
$ python --version
Python 3.9.5+
pandasを入れます。
$ pip install pandas
matplotlibを入れます。
$ pip install matplotlib
pvlibを入れます。
$ pip install pvlib
importします。
>>> import pandas as pd
>>> import matplotlib.pyplot as plt
>>> import pvlib
Apple SiliconのMacだとpvlibのimportが(内部のscipy.optimizeのimportで)失敗するので、次の魔法のコマンドで入れてください。
$ OPENBLAS="/opt/homebrew/opt/openblas" CFLAGS="-falign-functions=8 ${CFLAGS}" pip install scikit-learn --no-binary :all:
時刻のインデックスを取得しよう！
タイムゾーン日本で2021年の間の10秒ごとの時刻を取得する。
>>> time = pd.date_range('2021-01-01 00:00:00', '2021-12-31 23:59:59', freq='10s', tz='Asia/Tokyo')
>>> time
DatetimeIndex(['2021-01-01 00:00:00+09:00', '2021-01-01 00:00:10+09:00',
'2021-01-01 00:00:20+09:00', '2021-01-01 00:00:30+09:00',
'2021-01-01 00:00:40+09:00', '2021-01-01 00:00:50+09:00',
'2021-01-01 00:01:00+09:00', '2021-01-01 00:01:10+09:00',
'2021-01-01 00:01:20+09:00', '2021-01-01 00:01:30+09:00',
...
'2021-12-31 23:58:20+09:00', '2021-12-31 23:58:30+09:00',
'2021-12-31 23:58:40+09:00', '2021-12-31 23:58:50+09:00',
'2021-12-31 23:59:00+09:00', '2021-12-31 23:59:10+09:00',
'2021-12-31 23:59:20+09:00', '2021-12-31 23:59:30+09:00',
'2021-12-31 23:59:40+09:00', '2021-12-31 23:59:50+09:00'],
dtype='datetime64[ns, Asia/Tokyo]', length=3153600, freq='10S')
太陽位置を取得しよう！
日本標準時、兵庫県明石市の場所を指定して、時刻に対応する太陽位置を取得します。
zenithが天頂角、azimuthが方位角です。
>>> solarposition = pvlib.solarposition.get_solarposition(time, 35, 135)
solarposition
apparent_zenith zenith apparent_elevation elevation azimuth equation_of_time
2021-01-01 00:00:00+09:00 168.009062 168.009062 -78.009062 -78.009062 356.390121 -3.254458
2021-01-01 00:00:10+09:00 168.011146 168.011146 -78.011146 -78.011146 356.574366 -3.254512
2021-01-01 00:00:20+09:00 168.013121 168.013121 -78.013121 -78.013121 356.758671 -3.254567
2021-01-01 00:00:30+09:00 168.014986 168.014986 -78.014986 -78.014986 356.943031 -3.254622
2021-01-01 00:00:40+09:00 168.016741 168.016741 -78.016741 -78.016741 357.127443 -3.254676
... ... ... ... ... ... ...
2021-12-31 23:59:10+09:00 168.019317 168.019317 -78.019317 -78.019317 355.612222 -3.119232
2021-12-31 23:59:20+09:00 168.021864 168.021864 -78.021864 -78.021864 355.796464 -3.119288
2021-12-31 23:59:30+09:00 168.024301 168.024301 -78.024301 -78.024301 355.980777 -3.119343
2021-12-31 23:59:40+09:00 168.026629 168.026629 -78.026629 -78.026629 356.165159 -3.119398
2021-12-31 23:59:50+09:00 168.028847 168.028847 -78.028847 -78.028847 356.349608 -3.119453
[3153600 rows x 6 columns]
試しに、お正月の太陽位置を抜き出してみます。indexでお正月の日付を指定するだけです。
>>> syogatsu = solarposition['2021-01-01']
>>> syogatsu
apparent_zenith zenith apparent_elevation elevation azimuth equation_of_time
2021-01-01 00:00:00+09:00 168.009062 168.009062 -78.009062 -78.009062 356.390121 -3.254458
2021-01-01 00:00:10+09:00 168.011146 168.011146 -78.011146 -78.011146 356.574366 -3.254512
2021-01-01 00:00:20+09:00 168.013121 168.013121 -78.013121 -78.013121 356.758671 -3.254567
2021-01-01 00:00:30+09:00 168.014986 168.014986 -78.014986 -78.014986 356.943031 -3.254622
2021-01-01 00:00:40+09:00 168.016741 168.016741 -78.016741 -78.016741 357.127443 -3.254676
... ... ... ... ... ... ...
2021-01-01 23:59:10+09:00 167.905916 167.905916 -77.905916 -77.905916 354.982736 -3.724131
2021-01-01 23:59:20+09:00 167.908836 167.908836 -77.908836 -77.908836 355.165293 -3.724185
2021-01-01 23:59:30+09:00 167.911647 167.911647 -77.911647 -77.911647 355.347931 -3.724239
2021-01-01 23:59:40+09:00 167.914350 167.914350 -77.914350 -77.914350 355.530647 -3.724293
2021-01-01 23:59:50+09:00 167.916945 167.916945 -77.916945 -77.916945 355.713439 -3.724347
[8640 rows x 6 columns]
最初が01-01、最後が12-31で終わっているようなので、抜き出せているようです。
せっかくなので天頂角を表示してみます。
天頂角は天頂からの太陽の角度です。つまり天頂角度0は真上で、天頂角度90度は真横になります。天頂角90度になるタイミングは二度あり、これは日の入りと日の出なわけです。
このグラフでは90度未満のときが日の出ている時間、つまり日中となります。
>>> syogatsu.zenith.plot()
>>> plt.show()
お昼の時間帯に下がっていって(=>太陽が上っていっている)、ちゃんと表示できているようです。
先程話したとおり、zenithが天頂角で、これが90度未満のときが日中です。
この条件のままそのまま抜き出します。
>>> daytime = solarposition[solarposition.zenith < 90]
>>> daytime
apparent_zenith zenith apparent_elevation elevation azimuth equation_of_time
2021-01-01 07:12:40+09:00 89.507743 89.987061 0.492257 0.012939 118.508416 -3.396202
2021-01-01 07:12:50+09:00 89.482053 89.957075 0.517947 0.042925 118.532309 -3.396257
2021-01-01 07:13:00+09:00 89.456325 89.927096 0.543675 0.072904 118.556210 -3.396311
2021-01-01 07:13:10+09:00 89.430556 89.897124 0.569444 0.102876 118.580119 -3.396366
2021-01-01 07:13:20+09:00 89.404749 89.867159 0.595251 0.132841 118.604037 -3.396420
... ... ... ... ... ... ...
2021-12-31 16:52:50+09:00 89.407694 89.870576 0.592306 0.129424 241.312706 -2.977910
2021-12-31 16:53:00+09:00 89.433466 89.900506 0.566534 0.099494 241.336637 -2.977965
2021-12-31 16:53:10+09:00 89.459199 89.930443 0.540801 0.069557 241.360560 -2.978021
2021-12-31 16:53:20+09:00 89.484893 89.960386 0.515107 0.039614 241.384474 -2.978076
2021-12-31 16:53:30+09:00 89.510547 89.990337 0.489453 0.009663 241.408380 -2.978131
[1581850 rows x 6 columns]
最初がお正月の日の出7:12, 最後が大晦日の日の入り16:53なので正しそうです。ちゃんと抜き出せているようです。
日の出、日の入り、日中時間を調べよう！
やり方は単純です。先程抜き出した日中データから、日毎に最小時刻、最大時刻、行数を計算するだけです。
日中は日の出ている時間ですから、最小時刻は日の出、最大時刻は日の入り、その間の時刻の数を数えると日中時間になります。
pandasでもこの条件に従い、かんたんに書き下すだけです。
>>> daily_summary = daytime.resample('D').apply(lambda df: pd.Series({'sunrise': df.index[0], 'sunset': df.index[-1], 'daytime': df['zenith'].count()/(60*6)}))
>>> daily_summary
sunrise sunset daytime
2021-01-01 00:00:00+09:00 2021-01-01 07:12:40+09:00 2021-01-01 16:54:30+09:00 9.700000
2021-01-02 00:00:00+09:00 2021-01-02 07:12:50+09:00 2021-01-02 16:55:10+09:00 9.708333
2021-01-03 00:00:00+09:00 2021-01-03 07:13:00+09:00 2021-01-03 16:56:00+09:00 9.719444
2021-01-04 00:00:00+09:00 2021-01-04 07:13:10+09:00 2021-01-04 16:56:50+09:00 9.730556
2021-01-05 00:00:00+09:00 2021-01-05 07:13:10+09:00 2021-01-05 16:57:40+09:00 9.744444
... ... ... ...
2021-12-27 00:00:00+09:00 2021-12-27 07:11:10+09:00 2021-12-27 16:50:40+09:00 9.661111
2021-12-28 00:00:00+09:00 2021-12-28 07:11:30+09:00 2021-12-28 16:51:20+09:00 9.666667
2021-12-29 00:00:00+09:00 2021-12-29 07:11:50+09:00 2021-12-29 16:52:00+09:00 9.672222
2021-12-30 00:00:00+09:00 2021-12-30 07:12:10+09:00 2021-12-30 16:52:40+09:00 9.677778
2021-12-31 00:00:00+09:00 2021-12-31 07:12:20+09:00 2021-12-31 16:53:30+09:00 9.688889
[365 rows x 3 columns]
いきなり夏至の日を計算するのもいいですが、せっかくなので1年間の傾向を見ておきましょう。
とりあえず何も考えず、日中時間を表示してみます。
>>> daily_summary.daytime.plot.bar()
>>> plt.show()
ちょっと目が細かすぎて見ずらいので、月の平均にして表示してみます。
>>> daily_summary.resample('M').mean().daytime.plot.bar()
>>> plt.show()
x軸のメモリの値が途切れてしまっていますが、どうやら、6月が一番日中時間が長くて、逆に12月が一番日中時間が短いようです。
これで傾向がわかったので、では夏至を計算してみましょう。
>>> geshi_day = daily_summary.daytime.idxmax()
>>> geshi_day
Timestamp('2021-06-21 00:00:00+0900', tz='Asia/Tokyo', freq='D')
どうやら夏至の日は6/21のようです。
Google先生もそうだと言っているので、正解のようです。
夏至の日がわかったので、一応日の出時刻、日の入り時刻、日中時間も調べておきましょう。
>>> daily_summary.loc[geshi_day]
sunrise 2021-06-21 04:51:10+09:00
sunset 2021-06-21 19:12:30+09:00
daytime 8.615
Name: 2021-06-21 00:00:00+09:00, dtype: object
どうやら午前4時51分日の出、午後7時12分日の入りで、14.35時間が日中のようです。
やったー計算できたー
イェーイ
おわり