Zenn
🌅

BigQueryで初日の出を拝む

2020/12/23に公開
GCP
BigQuery
#
アドベントカレンダー
#
udf
#
astronomy
tech

https://qiita.com/advent-calendar/2020/bigquery

本稿は、緯度経度とタイムスタンプを使って日の出時刻を取得する UDF がほしい！というあなたに送る記事です。（永続 UDF の作り方と一般公開データセットの使い方を簡単に紹介します）

日の出時刻を取得する UDF を作成する

Astronomy Answersの計算式を参考に UDF を作成しました。
注意点としては、日時に関する情報は DATE 型ではなく TIMESTAMP 型で扱うことです。DATE 型はタイムゾーンの情報を持たないので TIMESTAMP で渡して UTC の TIMESTAMP で返すのが潔いかと思いました(よしなに変換しましょう)。

日の出時刻

CREATE OR REPLACE FUNCTION udf.SUNRISE_TIMESTAMP(ts TIMESTAMP, latitude FLOAT64, longitude FLOAT64) AS (
  TIMESTAMP_MILLIS(CAST(ROUND((((2451545 + (0.0009 + (0 + (ACOS(-1) / 180 * - longitude)) / (2 * (ACOS(-1))) + (ROUND((UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545) - 0.0009 - (ACOS(-1) / 180 * - longitude) / (2 * (ACOS(-1)))))) + 0.0053 * sin(((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) - 0.0069 * sin(2 * (((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545))) + ((ACOS(-1) / 180) * (1.9148 * sin(((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) + 0.02 * sin(2 * ((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) + 0.0003 * sin(3 * ((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))))) + ((ACOS(-1) / 180) * 102.9372) + (ACOS(-1))))) - ((2451545 + (0.0009 + ((acos((sin((-0.833 + 0) * (ACOS(-1) / 180)) - sin((ACOS(-1) / 180 * latitude)) * sin((asin(sin(0) * cos((ACOS(-1) / 180 * 23.4397)) + cos(0) * sin((ACOS(-1) / 180 * 23.4397)) * sin((((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545))) + ((ACOS(-1) / 180) * (1.9148 * sin(((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) + 0.02 * sin(2 * ((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) + 0.0003 * sin(3 * ((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))))) + ((ACOS(-1) / 180) * 102.9372) + (ACOS(-1)))))))) / (cos((ACOS(-1) / 180 * latitude)) * cos((asin(sin(0) * cos((ACOS(-1) / 180 * 23.4397)) + cos(0) * sin((ACOS(-1) / 180 * 23.4397)) * sin((((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545))) + ((ACOS(-1) / 180) * (1.9148 * sin(((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) + 0.02 * sin(2 * ((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) + 0.0003 * sin(3 * ((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))))) + ((ACOS(-1) / 180) * 102.9372) + (ACOS(-1)))))))))) + (ACOS(-1) / 180 * - longitude)) / (2 * (ACOS(-1))) + (ROUND((UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545) - 0.0009 - (ACOS(-1) / 180 * - longitude) / (2 * (ACOS(-1)))))) + 0.0053 * sin(((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) - 0.0069 * sin(2 * (((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545))) + ((ACOS(-1) / 180) * (1.9148 * sin(((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) + 0.02 * sin(2 * ((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) + 0.0003 * sin(3 * ((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))))) + ((ACOS(-1) / 180) * 102.9372) + (ACOS(-1))))) - (2451545 + (0.0009 + (0 + (ACOS(-1) / 180 * - longitude)) / (2 * (ACOS(-1))) + (ROUND((UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545) - 0.0009 - (ACOS(-1) / 180 * - longitude) / (2 * (ACOS(-1)))))) + 0.0053 * sin(((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) - 0.0069 * sin(2 * (((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545))) + ((ACOS(-1) / 180) * (1.9148 * sin(((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) + 0.02 * sin(2 * ((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) + 0.0003 * sin(3 * ((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))))) + ((ACOS(-1) / 180) * 102.9372) + (ACOS(-1))))))) + 0.5 - 2440588) * 1000 * 60 * 60 * 24) AS INT64))
);

日本各地の 2021 年の初日の出時刻を求める

BigQuery の一般公開データセットには気象データがあります。なかでもghcnd_stationsというテーブルには世界各地の気象観測ポイントの緯度経度などの情報が含まれています。これを利用して、全国各地の 2021 年の初日の出時刻を求めてみたいと思います。
国際地点番号表を見ると日本の wmoid は 47401〜47991 に割り当てられていることが分かります。

SELECT
  name,
  TIME(udf.SUNRISE_TIMESTAMP(TIMESTAMP("2021-01-01 15:30:00+09"), latitude, longitude), "Asia/Tokyo") AS sunrise
FROM
  `bigquery-public-data.ghcn_d.ghcnd_stations`
WHERE
  wmoid BETWEEN 47401 AND 47991
ORDER BY
  sunrise

このクエリを実行すると、南鳥島が日本国内最速で初日の出が見られることが分かります(本州では銚子となっていますが、高度などを考慮していないのでおそらく富士山山頂が一番早いでしょう)。
クエリ結果

ちなみに最も遅いのは与那国島で7:33でした。寝坊しても急いで与那国島に行けば初日の出に間に合うかもしれません。

2021 年は家でおとなしく正月を迎えますが、2022 年は南鳥島で初日の出を見たいですね。

おまけ

日の入り時刻

CREATE OR REPLACE FUNCTION udf.SUNSET_TIMESTAMP(ts TIMESTAMP, latitude FLOAT64, longitude FLOAT64) AS (
  TIMESTAMP_MILLIS(CAST(ROUND(((2451545 + (0.0009 + ((acos((sin((-0.833 + 0) * (ACOS(-1) / 180)) - sin((ACOS(-1) / 180 * latitude)) * sin((asin(sin(0) * cos((ACOS(-1) / 180 * 23.4397)) + cos(0) * sin((ACOS(-1) / 180 * 23.4397)) * sin((((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545))) + ((ACOS(-1) / 180) * (1.9148 * sin(((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) + 0.02 * sin(2 * ((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) + 0.0003 * sin(3 * ((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))))) + ((ACOS(-1) / 180) * 102.9372) + (ACOS(-1)))))))) / (cos((ACOS(-1) / 180 * latitude)) * cos((asin(sin(0) * cos((ACOS(-1) / 180 * 23.4397)) + cos(0) * sin((ACOS(-1) / 180 * 23.4397)) * sin((((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545))) + ((ACOS(-1) / 180) * (1.9148 * sin(((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) + 0.02 * sin(2 * ((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) + 0.0003 * sin(3 * ((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))))) + ((ACOS(-1) / 180) * 102.9372) + (ACOS(-1)))))))))) + (ACOS(-1) / 180 * - longitude)) / (2 * (ACOS(-1))) + (ROUND((UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545) - 0.0009 - (ACOS(-1) / 180 * - longitude) / (2 * (ACOS(-1)))))) + 0.0053 * sin(((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) - 0.0069 * sin(2 * (((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545))) + ((ACOS(-1) / 180) * (1.9148 * sin(((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) + 0.02 * sin(2 * ((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))) + 0.0003 * sin(3 * ((ACOS(-1) / 180) * (357.5291 + 0.98560028 * (UNIX_MILLIS(ts) / (1000 * 60 * 60 * 24) - 0.5 + 2440588 - 2451545)))))) + ((ACOS(-1) / 180) * 102.9372) + (ACOS(-1))))) + 0.5 - 2440588) * 1000 * 60 * 60 * 24) AS INT64))
);

Discussion