続、てほどの目新しい話はないんだけれども。これの続きである。
前回のは「やってみたらどうなるんだろう」をやってみただけだったんだけれど、「ほんの少しの本気で結構まともな計測出来るんじゃないかしら」と思い始めてて。
おさらい。
- マグネットコンパス。
- 発想はこれだけ:
- GeographicLib の pure Python 版を利用して二点間経路の方位を求める。
ただ、前回のは「地点 B での計測」が「磁偏差計測の確定」にしてなかったので、使いにくかったのよね。あともっと刹那的でいい。「管理名」で二点のペア管理しようとしてたけど、「直前の地点 A 計測」だけが永続化されて、地点 B 計測、でいい。つまり前回あった「Draw」ボタンはいらなくて、計測終了とともにコンパス描画でいい。
てわけでイキナリ新版:
android.txt
1 title=Compass with Measuring Declination
2 author=hhsprings
3 orientation=all
measure_declination.py
1 # -*- coding: utf-8 -*-
2 import os
3 import pickle
4 from math import degrees, radians, sin, cos
5 import cmath
6 import sys
7 sys.path.append("/sdcard/kivy-my-site-packages")
8 from geographiclib.geodesic import Geodesic
9 WGS84 = Geodesic.WGS84
10
11 _DBFILE = "./_origin_loc.pickle"
12
13 def mean_of_angle(a):
14 rho = sum(cmath.rect(1, radians(d)) for d in a)
15 if abs(rho) < 1.0e-5:
16 return float('nan') # or None
17 rho /= len(a)
18 th = cmath.phase(rho)
19 return degrees(th)
20
21 def norm_deg(d):
22 d %= 360
23 if d > 180:
24 # 350 -> -10
25 # 190 -> -170
26 return d - 360.
27 return d
28
29 class MagneticDeclinationMeasure(object):
30 _record = []
31 mode = 'O'
32 s12 = 0.0
33 declination = 0.0
34
35 def start(self, mode):
36 self._record = []
37 self.mode = mode
38
39 def record(self, lon, lat, compass_azimuth):
40 self._record.append((lon, lat, compass_azimuth))
41
42 def finish(self):
43 if self.mode == 'O':
44 # dump means of lon and lat
45 result = (
46 mean_of_angle([_[0] for _ in self._record]), # lon
47 mean_of_angle([_[1] for _ in self._record]), # lat
48 )
49 pickle.dump(result, open(_DBFILE, "wb"))
50 else:
51 # means of lon, lat, compass_azimuth at current location
52 dlon = mean_of_angle([_[0] for _ in self._record])
53 dlat = mean_of_angle([_[1] for _ in self._record])
54 cazi = mean_of_angle([_[2] for _ in self._record])
55
56 # load means of lon and lat at origin
57 olon, olat = pickle.load(open(_DBFILE, "rb"))
58
59 # Inverse
60 iv = WGS84.Inverse(olat, olon, dlat, dlon)
61 azi2, self.s12 = iv["azi2"], iv["s12"]
62
63 # azi2: clockwise from TN [degrees]
64 # cazi: heading direction, clockwise from MN [degrees]
65 self.declination = norm_deg(azi2) - norm_deg(cazi)
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 measure_declination import MagneticDeclinationMeasure, norm_deg
12
13 kv = '''
14 BoxLayout:
15 orientation: 'vertical'
16
17 face: face
18 details: details
19 bottom_pane: bottom_pane
20
21 Label:
22 id: details
23 text: app.disp
24
25 Widget:
26 id: face
27 r: min(root.size) * 0.9 / 2
28
29 BoxLayout:
30 size_hint_y: None
31 height: '48dp'
32 padding: '4dp'
33 id: bottom_pane
34
35 ToggleButton:
36 text: 'Origin' if self.state == 'normal' else 'Destination'
37 on_state:
38 app.mode = 'D' if self.state == 'down' else 'O'
39
40 ToggleButton:
41 text: 'Measure Start' if self.state == 'normal' else 'Measure Stop'
42 on_state:
43 app.start() if self.state == 'down' else app.stop()
44 '''
45
46 class CompassWithMeasDecl(App):
47
48 disp = StringProperty()
49 mode = StringProperty('O')
50 _clock_event = None
51 _clock_event_2 = None
52
53 def build(self):
54 self._measure = MagneticDeclinationMeasure()
55
56 self.gps = gps
57 self.gps.configure(on_location=self.on_location,
58 on_status=self.on_status)
59
60 compass.enable()
61
62 return Builder.load_string(kv)
63
64 def start(self):
65 if self._clock_event_2:
66 self._clock_event_2.cancel()
67
68 self.lon = None
69 self.lat = None
70
71 self.gps.start()
72 self._measure.start(self.mode)
73 self._clock_event = Clock.schedule_interval(self._do_record, 1.)
74
75 def stop(self):
76 self._clock_event.cancel()
77 self.gps.stop()
78 self.lon = None
79 self.lat = None
80 self._measure.finish()
81 if self._measure.mode == 'D':
82 self._clock_event_2 = Clock.schedule_interval(self._do_draw, 1.)
83
84 def _do_record(self, *args, **kwarg):
85
86 lon = self.lon
87 lat = self.lat
88 azimuth = compass.orientation[0]
89 self.disp = "({}, {}, {})".format(lon, lat, azimuth)
90 if lon is None or lat is None or azimuth is None:
91 return
92 self._measure.record(lon, lat, azimuth)
93
94 def _do_draw(self, *args, **kwarg):
95 azimuth = compass.orientation[0]
96 #
97 # ------------------------------------------------
98 self.root.face.canvas.clear()
99 #
100 magnorth = radians(360 - azimuth)
101 truenorth = radians(360 - (azimuth + self._measure.declination))
102
103 self.disp = "TN = {}\nMN = {}\nDECL = {}".format(
104 norm_deg(degrees(truenorth)),
105 norm_deg(degrees(magnorth)),
106 self._measure.declination)
107
108 with self.root.face.canvas:
109 Color(0.2, 0.2, 0.7)
110 center_x, center_y = self.root.face.center_x, self.root.face.center_y
111 Line(points=[
112 center_x, center_y,
113 center_x + 0.8 * self.root.face.r * sin(truenorth),
114 center_y + 0.8 * self.root.face.r * cos(truenorth)
115 ], width=5, cap="round")
116
117 Color(0.7, 0.2, 0.2)
118 center_x, center_y = self.root.face.center_x, self.root.face.center_y
119 Line(points=[
120 center_x, center_y,
121 center_x + 0.8 * self.root.face.r * sin(magnorth),
122 center_y + 0.8 * self.root.face.r * cos(magnorth)
123 ], width=5, cap="round")
124 # ------------------------------------------------
125
126 @mainthread
127 def on_location(self, **kwargs):
128 self.lon = float(kwargs['lon'])
129 self.lat = float(kwargs['lat'])
130 #self.alt = float(kwargs['altitude'])
131 # azimuth = compass.orientation[0]
132
133 @mainthread
134 def on_status(self, stype, status):
135 pass
136
137 if __name__ == '__main__':
138 CompassWithMeasDecl().run()
地点B計測終了で、以下のように True North (青) 方位がくっついたコンパスとして動き始める:
こんなん誰が使いたいと思うんだ、とは思うけど一応操作説明:
- まっすぐな道を選び、静かに立って「Measure Start」、(None, None, …) でない表示は記録されないので、経度緯度がまともな値になってからしばらく記録。気が済んだら「Measure Stop」。これでこの地点の「経度緯度」がファイルに書き出される。
- そのまっすぐな道をある程度進む。GPS の誤差は確か10mくらいだったと思うので、これよりはずっと大きな距離を歩く。
- そのまっすぐな道を歩いた先で、デバイスの方位と道の方位を正確に合わせて、左側のトグル(「Origiin」表示)を押して「Destination」にする。
- そのまま静かに動かずに、最初の地点と同じ要領で「Measure Start」、同じく (None, None, …) 状態でなくなってしばらく記録したら「Measure Stop」。
- これで磁偏差計測が確定したので、あとは自由にぐるぐる回ってヨシ。
さっきやってみたら西偏6度くらいの値になった。うん、たぶんそんなもん。(キャプチャの絵は Origin のデータが残ってたので「道に合わせる」という正しい手順を取らずにいい加減にやったヤツ。)
8/4追記:
西偏6度くらいの値になった、よかったよかった、と、少し何回かやってみたら、「案の定値がぶれぶれ」であった。
多摩川大橋で試してみて西偏14度なんてのが出て、あぁうまくいかんなぁ、と思い、帰りに同じようにやってみても今度は東偏で出たりして。自宅前でも東偏3度とかそんな値にもなったりした。
True North そのものは(GPS 位置精度と GeographicLib 演算精度内において)確実なのだが、Magnetic North は、一つにはデバイスのセンサーの精度の問題があり(だから平均を採っているのだがそれでもなお)、もう一つは現実の磁場そのものが本当に刻々と変化しているということ。磁場は色んなものに左右されるが、むろん地上のちょっとした設備やらなんやらの影響も受ける。要するに本当にわからない。ので、こうやって測ってる値のもっともらしさを検証する術がないのね。
結局前回同様「やっぱ簡単にはいかんなぁ」てことが「わかった」に過ぎないのでした。おしまい。
