BigQuery에서 설날 해돋이를 바라보다

이 글은 BigQuery Advent Calendar 2020 5일째 글(대타)이다.
https://qiita.com/advent-calendar/2020/bigquery
본고는 UDF가 위도 경도와 시간 스탬프를 이용하여 일출 시간을 얻기를 희망합니다!이것은 너에게 준 문장이다.(영구 UDF의 제작 방법과 일반 공개 데이터 세트의 사용 방법을 간단히 소개)

일출 시간을 가져오는 UDF 생성하기

Astronomy Answers의 계산 공식을 참고하여 UDF를 만들었다.
날짜와 시간에 대한 정보는 DATE형이 아니라 TIMESTAMP형이라는 점에 유의해야 한다.데이트형은 시간대에 대한 정보가 없기 때문에 저는 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))
);