Python - 地球自転速度補正値 ΔT の取得(USNO から)!

Updated:


前回、 USNO(The United States Naval Observatory; アメリカ海軍天文台)から地球自転速度補正値 ΔT をする処理を Ruby で実装しました。

今回は、同じ処理を Python で実装してみました。

0. 前提条件

  • Python 3.6.5 での動作を想定。
  • ΔT 一覧を取得するページは「こちら」。

1. Python スクリプトの作成

  • Shebang ストリング(1行目)では、フルパスでコマンド指定している。(当方の慣習
  • User Agent UA の PG(Bot) 名やメールアドレスは自分のもので置き換えること。
    (User Agent の設定は Web スクレイプする者のマナーであると個人的に認識している)

File: get_delta_t_usno.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
#! /usr/local/bin/python3.6
"""
USNO(アメリカ海軍天文台) から delta T(ΔT) を取得
: 確定データは
  [Delta T determinations](http://maia.usno.navy.mil/ser7/deltat.data)
  から、推定データは
  [Delta T predictions](http://maia.usno.navy.mil/ser7/deltat.preds)
  から取得する。

  date          name            version
  2018.06.20    mk-mode.com     1.00 新規作成

Copyright(C) 2018 mk-mode.com All Rights Reserved.
---
引数 : UTC (YYYYMMDD 形式)
       * 無指定ならシステム日時を UTC とみなす
"""
from datetime import datetime
import re
import requests
import sys
import traceback


class GetDeltaTUsno:
    USAGE = "[USAGE] ./get_delta_t_usno.py [yyyymmdd]"
    URL_D = "http://maia.usno.navy.mil/ser7/deltat.data"
    URL_P = "http://maia.usno.navy.mil/ser7/deltat.preds"
    UA = (
        "mk-mode Bot (by Python/{}.{}.{}, "
        "Administrator: postmaster@mk-mode.com)"
    ).format(
        sys.version_info.major,
        sys.version_info.minor,
        sys.version_info.micro
    )
    PTN_0 = re.compile(r'[\r\n]+')
    PTN_1 = re.compile(r'^\s+')
    PTN_2 = re.compile(r'\s+')
    PTN_3 = re.compile(r'[\d\.]+')
    PTN_4 = re.compile(r'\.')

    def __init__(self):
        self.__get_arg()
        self.data = []    # delta t 一覧
        self.delta_t = 0  # delta t 初期値
        self.dp = ""      # 確定/推定値区分(推定値なら "*")

    def exec(self):
        try:
            self.__get_determination()  # 確定データ取得
            self.__get_prediction()     # 推定データ取得
            self.__get_delta_t()        # delta_t 取得
            self.__display()            # 結果出力
        except Exception as e:
            raise

    def __get_arg(self):
        """ argument getting """
        try:
            if len(sys.argv) > 1:
                ymd = sys.argv[1]
            else:
                ymd = datetime.now().strftime("%Y%m%d")
            if not(re.search(r"^\d{8}$", ymd)):
                print(self.USAGE)
                sys.exit(0)
            self.date = "{}-{}-{}".format(ymd[0:4], ymd[4:6], ymd[6:8])
        except Exception as e:
            raise

    def __get_determination(self):
        """ 確定データ取得 """
        try:
            headers = {'user-agent': self.UA}
            res = requests.get(self.URL_D, headers=headers)
            html = res.text.strip()
            for line in re.split(self.PTN_0, html):
                line = re.sub(self.PTN_1, "", line)
                items = re.split(self.PTN_2, re.sub(self.PTN_1, "", line))
                date = "{:04d}-{:02d}-{:02d}".format(
                    int(items[0]), int(items[1]), int(items[2])
                )
                self.date_max = date
                self.data.append([date, float(items[-1]), "D"])
        except Exception as e:
            raise

    def __get_prediction(self):
        """ 推定データ取得 """
        try:
            headers = {'user-agent': self.UA}
            res = requests.get(self.URL_P, headers=headers)
            html = res.text.strip()
            for line in re.split(self.PTN_0, html):
                line = re.sub(self.PTN_1, "", line)
                items = re.split(self.PTN_2, re.sub(self.PTN_1, "", line))
                if not(re.match(self.PTN_3, items[0])):
                    continue
                year, year_frac = re.split(self.PTN_4, items[0])
                year = int(year)
                if year_frac == "00":
                    month = 1
                elif year_frac == "25":
                    month = 4
                elif year_frac == "50":
                    month = 7
                elif year_frac == "75":
                    month = 10
                date = "{:04d}-{:02d}-01".format(year, month)
                if date <= self.date_max:
                    continue
                self.data.append([date, float(items[1]), "P"])
        except Exception as e:
            raise

    def __get_delta_t(self):
        """ delta t 取得 """
        try:
            for i, d in enumerate(reversed(self.data)):
                if d[0] <= self.date:
                    if i == 0:
                        break
                    self.delta_t = d[1]
                    if d[2] == "P":
                        self.dp = "*"
                    break
        except Exception as e:
            raise

    def __display(self):
        """ 結果出力 """
        try:
            print("[{}] delta_T = {} sec. {}".format(
                self.date, self.delta_t, self.dp)
            )
        except Exception as e:
            raise


if __name__ == '__main__':
    try:
        obj = GetDeltaTUsno()
        obj.exec()
    except Exception as e:
        traceback.print_exc()
        sys.exit(1)

2. Python スクリプトの実行

まず、実行権限を付与。

$ chmod +x get_delta_t_usno.py

そして、コマンドライン引数に UTC(協定世界時)を YYYYMMDD の書式で指定して実行する。(引数無指定なら、システム日付を UTC とみなす)

$ ./get_delta_t_usno.py 20180620
[2018-06-20] delta_T = 69.14 sec. *
  • 推定値の場合は行末に * が付与される。

人工衛星の正確な軌道計算等に利用できるでしょう。

ちなみに、今回の方法を応用して、 IERS(国際地球回転観測事業; International Earth Rotation and Reference systems Service)提供の Polar Motion(地球の極運動)のデータも取得できます。

以上。





 

Sponsored Link

 

Comments