夏至の日なので（だったので）、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時間が日中のようです。

やったー計算できたー
イェーイ

おわり