2地点あれば True North がわかる

たりめーだ。

誰の役に立つんだぁ、てのはいつもの話。ワタシの興味だけでやっておる。

「磁偏差」が「ワタシにどう効いておるか」ってのは、本当のところは自分で測ってみないとよくわからんのだよねぇ、を突き詰めて、「自分で測れないかなぁ」と思ったのであった。

マグネットコンパスは手にしている

True North を知るには、「Magnet North と磁偏差があればわかる」。いや…その磁偏差を測りたいんですってば。

True North を知るには、アゲイン:

まっすぐな道を見つけて、両端の緯度経度を記録、その二地点の azimuth を測るってわけである。これで「ワタシは True North に対して azimuth 度に進んできた」ことがわかるってわけだ。

2地点の関係を知るためのライブラリとして、GeographicLib は公式に pure Python 版をリリースしている。site-packages ライブラリを制御出来ない kivyLauncher で使うためには是が非でも pure Python でなければならないので、もっと使いやすい pyproj よりも今はありがたい。置き場所にはちょっと困るのだが、「オレオレ標準」として /sdcard/kivy-my-site-packages なんて場所を作ってここに放り込み、このように使うことに:

1 import sys
2 sys.path.append("/sdcard/kivy-my-site-packages")
3 from geographiclib.geodesic import Geodesic
4 WGS84 = Geodesic.WGS84

あとは「二地点」の管理が必要だ。地点Aと地点Bはアプリケーションをずっと起動しっぱなしに出来るとは限らない。というか 500m 歩きながら起動しっぱなしなんてやだしムリ。ので、sqlite に記録しよう:

loc_log.py
 1 # -*- coding: utf-8 -*-
 2 import os
 3 import sqlite3
 4 from math import degrees, radians, sin, cos
 5 import cmath
 6 
 7 _DBFILE = "./loc_log.dat"
 8 
 9 def mean_of_angle(a):
10     rho = sum(cmath.rect(1, radians(d)) for d in a)
11     if abs(rho) < 1.0e-5:
12         return float('nan')  # or None
13     rho /= len(a)
14     th = cmath.phase(rho)
15     return degrees(th)
16 
17 class LocLog(object):
18     def __init__(self):
19         self._conn = self._open_connection()
20 
21     def _open_connection(self):
22         conn = sqlite3.connect(_DBFILE)
23         cur = conn.cursor()
24         try:
25             cur.executescript('''\
26 CREATE TABLE t_loc_log (
27     mng_name TEXT,
28     lon FLOAT,
29     lat FLOAT,
30     compass_azimuth FLOAT);
31 ''')
32             conn.commit()
33         except:
34             pass
35         cur.close()
36 
37         return conn
38 
39     def insert(self, mng_name, lon, lat, compass_azimuth):
40         cur = self._conn.cursor()
41         cur.execute(
42         """INSERT INTO t_loc_log(mng_name, lon, lat, compass_azimuth)
43            VALUES (?, ?, ?, ?)""", (
44                 mng_name, lon, lat, compass_azimuth))
45         self._conn.commit()
46         cur.close()
47 
48     def select(self, mng_name):
49         cur = self._conn.cursor()
50         cur.execute("""SELECT lon, lat, compass_azimuth FROM t_loc_log 
51                        WHERE mng_name = ?""", (mng_name,))
52         lons = []
53         lats = []
54         azis = []
55         for row in cur:
56             lons.append(row[0])
57             lats.append(row[1])
58             azis.append(row[2])
59         cur.close()
60         if lons:
61             return (
62                 mean_of_angle(lons),
63                 mean_of_angle(lats),
64                 mean_of_angle(azis))
65 
66 if __name__ == '__main__':
67     import sys
68     log = LocLog()
69     print(log.select(sys.argv[1]))

管理は「地点Aと地点B」というセットをひとまとまりにしたいので、mng_name というカラムを持っているのだが、「A」と「B」を区別する必要もあって、これは mng_name の頭に Origin の「O」と Destination の「D」を付与するルールにした。

主キーなしで延々追記していくようになっていて、これは一定時間の平均を取りたかったため。

あとは UI は、

  1. mng_name を入力するテキストボックス
  2. 地点 A なのか地点 B なのかを特定するトグルボタン
  3. 記録開始終了ボタン
  4. 描画ボタン

として、地点A で記録、地点Bで記録、地点Bで描画。

main.py
  1 from math import radians, degrees, sin, cos
  2 from kivy.lang import Builder
  3 from plyer import gps
  4 from kivy.app import App
  5 from kivy.properties import StringProperty
  6 from kivy.clock import Clock, mainthread
  7 from kivy.graphics import Color, Line
  8 from kivy.logger import Logger
  9 from my_compass import AndroidCompass
 10 compass = AndroidCompass()
 11 from loc_log import LocLog
 12 import sys
 13 sys.path.append("/sdcard/kivy-my-site-packages")
 14 from geographiclib.geodesic import Geodesic
 15 WGS84 = Geodesic.WGS84
 16 
 17 
 18 kv = '''
 19 BoxLayout:
 20     orientation: 'vertical'
 21 
 22     face: face
 23     details: details
 24     bottom_pane: bottom_pane
 25 
 26     Label:
 27         id: details
 28         text: app.disp
 29 
 30     Widget:
 31         id: face
 32         r: min(root.size) * 0.9 / 2
 33 
 34     BoxLayout:
 35         size_hint_y: None
 36         height: '48dp'
 37         padding: '4dp'
 38         id: bottom_pane
 39 
 40         TextInput:
 41             on_text:
 42                 app.mng_name = args[1].strip()
 43 
 44         ToggleButton:
 45             text: 'Origin' if self.state == 'normal' else 'Destination'
 46             on_state:
 47                 app.mode = 'D' if self.state == 'down' else 'O'
 48 
 49         ToggleButton:
 50             text: 'Log Start' if self.state == 'normal' else 'Log Stop'
 51             on_state:
 52                 app.start() if self.state == 'down' else app.stop()
 53 
 54         ToggleButton:
 55             text: 'Draw Start' if self.state == 'normal' else 'Draw Stop'
 56             on_state:
 57                 app.draw_start() if self.state == 'down' else app.draw_stop()
 58 '''
 59 
 60 class GpsCompassLogging(App):
 61 
 62     disp = StringProperty()
 63     mng_name = StringProperty()
 64     mode = StringProperty('O')
 65 
 66     def build(self):
 67         self._locLog = LocLog()
 68 
 69         self.gps = gps
 70         self.gps.configure(on_location=self.on_location,
 71                            on_status=self.on_status)
 72 
 73         compass.enable()
 74 
 75         return Builder.load_string(kv)
 76 
 77     def start(self):
 78         self.lon = None
 79         self.lat = None
 80 
 81         self.gps.start()
 82         self._clock_event = Clock.schedule_interval(self._do_record, 1.)
 83 
 84     def stop(self):
 85         self._clock_event.cancel()
 86         self.gps.stop()
 87         self.lon = None
 88         self.lat = None
 89 
 90     def _do_record(self, *args, **kwarg):
 91 
 92         lon = self.lon
 93         lat = self.lat
 94         azimuth = compass.orientation[0]
 95         self.disp = "({}, {}, {})".format(lon, lat, azimuth)
 96         if lon is None or lat is None or azimuth is None:
 97             return
 98         self._locLog.insert(self.mode + self.mng_name, lon, lat, azimuth)
 99 
100     def draw_start(self):
101         self._clock_event2 = Clock.schedule_interval(self._do_draw, 1.)
102         ores = self._locLog.select("O" + self.mng_name)
103         dres = self._locLog.select("D" + self.mng_name)
104         self.invres = None
105         if not ores or not dres:
106             return
107 
108         self.invres = WGS84.Inverse(ores[1], ores[0], dres[1], dres[0])
109 
110     def draw_stop(self):
111         self._clock_event2.cancel()
112 
113     def _do_draw(self, *args, **kwarg):
114         if not self.invres:
115             return
116 
117         azimuth = compass.orientation[0]
118         #
119         # ------------------------------------------------
120         self.root.face.canvas.clear()
121         truenorth = radians(-self.invres['azi2'])
122         magnorth = radians(360 - azimuth)
123         tn_mn_diff = degrees(magnorth - truenorth) % 360
124         if tn_mn_diff > 180:
125             # 350 -> -10
126             # 190 -> -170
127             tn_mn_diff -= 360.
128         self.disp = "sl2 = {}[m]\nTN = {}\nMN = {}\nMN - TN= {}".format(
129             self.invres['s12'],
130             degrees(truenorth),
131             degrees(magnorth) % 360,
132             tn_mn_diff)
133             
134         with self.root.face.canvas:
135             Color(0.2, 0.2, 0.7)
136             center_x, center_y = self.root.face.center_x, self.root.face.center_y
137             Line(points=[
138                     center_x, center_y,
139                     center_x + 0.8 * self.root.face.r * sin(truenorth),
140                     center_y + 0.8 * self.root.face.r * cos(truenorth)
141                     ], width=5, cap="round")
142 
143             Color(0.7, 0.2, 0.2)
144             center_x, center_y = self.root.face.center_x, self.root.face.center_y
145             Line(points=[
146                     center_x, center_y,
147                     center_x + 0.8 * self.root.face.r * sin(magnorth),
148                     center_y + 0.8 * self.root.face.r * cos(magnorth)
149                     ], width=5, cap="round")
150         # ------------------------------------------------
151 
152     @mainthread
153     def on_location(self, **kwargs):
154         self.lon = float(kwargs['lon'])
155         self.lat = float(kwargs['lat'])
156         #self.alt = float(kwargs['altitude'])
157         # azimuth = compass.orientation[0]
158 
159     @mainthread
160     def on_status(self, stype, status):
161         pass
162 
163 if __name__ == '__main__':
164     GpsCompassLogging().run()

繰り返しになるけれど、my_compass はここ

android.txt はテキトー:

1 title=Simple GPS, Compass Logging
2 author=hhsprings
3 orientation=all

こんな:

青が True North。赤が Magnetic North。上のテキストの sl2 は地点 A から地点 B までの距離。てわけで、100m 歩けば問題なく計測出来る。

やってみてどうかって?

微妙だよね。そもそも実際にやってみればわかることだけれど、「デバイスをまっすぐに進行方向に向ける」ことだって結構むつかしいんであって。(これに神経質になる必要があるのは地点B。)

あと結局のところは、感覚的に「確かに True North と Magnetic North はズレてる」ということだけはわかるものの、「平均モデルの西偏 7 度はどの程度正しいのか?」を知れるほどのもんではない。ちなみに磁場は太陽風の影響とかも受けますよ。それが効いてるのかどうかだって、これでわかるようなもんではない。

ケド、「確かにやってみれば」てのは結構大事で、理想と現実がどれほど違うのかを感覚的に知っとくのって、オモシロイと思う。