思いついたのですぐに試してみたくなり。
国土地理院の基盤地図情報を持ち歩けないか、それで気圧補正すれば気圧高度計が結構精度よくなるんじゃないのか、のお試しをどうしてもすぐにやりたくなる。
発想のおさらい:
- スマホの気圧センサーは結構安定して値が取れる(大きく暴れることがない)
- やってみたらかなり細かい地形の起伏に繊細に反応してることがわかった
- 標準大気モデルに従った気圧高度を計算してみて、国土地理院の基盤地図情報標高データ(DEM5は5mメッシュ)と合わせてみると、かなり精密に起伏には反応していた
- 市販の気圧高度計のように「高度がわかっている地点で気圧補正」と気温補正をすると、結構一致した。
- 基盤地図情報標高データは「バカでかい」のでスマホに全部入れて持ってくなんて論外だが、「高度がわかっている地点で気圧補正」のためにピンポイントで抽出してスマホに突っ込んどくのはアリ。
- GPS センサでとらえた現在位置が抽出して持っている標高データの範囲にヒットしたら、その正確な標高とその時点での気圧を記憶し、以後の気圧高度計算の補正値に使う。
DEM5 の抽出、だけれども、「国土地理院からダウンロードしてきた XML を…」からの説明はかなり煩わしいので省略。とにかく「オレオレ的に加工済み」のものからさらに numpy の ndarray のまんま npz に書き出しておき、このようにロード出来る形にした:
1 # -*- coding: utf-8 -*-
2 import os
3 from glob import glob
4 import numpy as np
5
6 # ... (省略) ...
7 for fn in glob("dems/*.npz"):
8 bn = os.path.splitext(os.path.basename(fn))[0]
9 Logger.info("processing npz: {}".format(bn))
10 _ = np.load(fn)
11 X, Y, elevs = _["X"], _["Y"], _["elevs"]
12 self._known_elevs[bn] = (X, Y, elevs)
この X, Y, elevs はオリジナルの基盤地図情報とは別に関係はなくてあくまでもオレ流で、たとえば basemap で:
1 from mpl_toolkits.basemap import Basemap
2 import matplotlib.pyplot as plt
3 import matplotlib.cm as cm
4
5 fig, ax = plt.subplots()
6 m = Basemap(
7 llcrnrlon=X[0] - (X[10] - X[0]) / 2.,
8 llcrnrlat=Y[0] - (Y[10] - Y[0]) / 2.,
9 urcrnrlon=X[-1] + (X[2 * 10] - X[0]),
10 urcrnrlat=Y[-1],
11 resolution='i',
12 lat_0=(Y[0] + Y[-1]) / 2.,
13 lon_0=(X[0] + X[-1]) / 2.,
14 projection='cass')
15 lons, lats = np.meshgrid(X, Y)
16 x, y = m(lons, lats)
17 m.drawcoastlines(color='0.9', linewidth=1)
18 circles = np.arange(34, 36, 1./3600*10)
19 m.drawparallels(circles, labels=[1, 0, 0, 0])
20 meridians = np.arange(139, 141, 1./3600*10)
21 m.drawmeridians(meridians, labels=[0, 0, 0, 1])
22 m.pcolor(x, y, elevs) #, cmap=cm.summer)
23 plt.savefig("ofuna_st.png", bbox_inches="tight")
なんてすれば:
となる、てこと。print(X), print(Y), print(elevs) するとつまりこんなデータ:
1 [ 139.52706473 139.52712054 139.52717634 139.52723214 139.52728795
2 139.52734375 139.52739955 139.52745536 139.52751116 139.52756696
3 139.52762277 139.52767857 139.52773437 139.52779018 139.52784598
4 139.52790179 139.52795759 139.52801339 139.5280692 139.528125
5 ...省略...
6 139.53487723 139.53493304 139.53498884 139.53504464 139.53510045
7 139.53515625 139.53521205 139.53526786 139.53532366 139.53537946
8 139.53543527]
9 [ 35.34963768 35.34969342 35.34974916 35.34980491 35.34986065
10 35.34991639 35.34997213 35.35002787 35.35008361 35.35013935
11 ...省略...
12 35.35744147 35.35749721 35.35755295 35.3576087 35.35766444
13 35.35772018 35.35777592 35.35783166 35.3578874 35.35794314
14 35.35799888]
15 [[ 8.31 8.26 8.41 ..., 12.81 12.82 12.92]
16 [ 8.24 8.15 8.24 ..., 12.94 12.94 13.08]
17 [ 8.16 8.25 8.4 ..., 13. 13.06 13.19]
18 ...,
19 [ 10.11 10.11 10.23 ..., 12.01 12.07 12.47]
20 [ 10.2 10.16 10.06 ..., 11.49 11.81 12.01]
21 [ 9.83 9.8 9.62 ..., 11.28 11.56 11.64]]
そもそもグリッドデータなので、この X, Y の外部化は合理的ではないしメモリに優しくもないんだけれど、今回はお試しということで、サボった。
さて。「既知の場所(DEM5を持ってる場所)にヒットしたら記憶」という都合、プログラムはサービスとして作らねばならぬ。そこらの詳細は今回は省略。核心部はこんななの:
observer.py
1 # -*- coding: utf-8 -*-
2 import math
3 import os
4 import sys
5 from datetime import datetime
6 from glob import glob
7 import numpy as np
8 from kivy.clock import Clock, mainthread
9 from kivy.logger import Logger
10 from plyer import gps
11 from gravity_sensor import GravitySensor
12 from environmental_sensors import pressureSensor
13
14
15 grav_sensor = GravitySensor()
16 _ERA = datetime(1970, 1, 1)
17
18
19 # values bellow are at base=0 (<11km)
20 _P0 = 1013.25 # static pressure (Pa) at MSL
21 _T0 = 273.15 + 15 # standard temperature (K) at MSL
22 _L0 = 6.49 / 1000. # standard temperature lapse rate (K/m) in ISA
23 _R = 8.31432 # universal gas constant in N·m /(mol·K)
24 _g0 = 9.80665 # gravitational acceleration in m/s**2
25 _M = 0.0289644 # molar mass of Earth's air in kg/mol
26
27 #
28 def a2p(h, delta_T=0):
29 t = _T0 + delta_T
30 return _P0 * math.pow((t / (t + _L0 * h)), _g0*_M/(_R*_L0))
31 #
32 def p2a(p, hA, pA, delta_T=0):
33 delta_p = pA - a2p(hA, delta_T)
34
35 t = _T0 + delta_T
36 return (t / _L0) * (math.pow((p - delta_p) / _P0, -_R*_L0/(_g0*_M)) - 1)
37 #
38
39 class Observer(object):
40 lon = None
41 lat = None
42 alt = None
43 delta_T = 0.0 # public property
44 _known_elevs = {}
45 _last_known = {
46 "time": datetime(1970, 1, 1),
47 "altitude": 0,
48 "pressure": 1013.25,
49 "lon": 0.0,
50 "lat": 0.0
51 }
52
53 def __init__(self):
54 for fn in glob("dems/*.npz"):
55 bn = os.path.splitext(os.path.basename(fn))[0]
56 Logger.info("processing npz: {}".format(bn))
57 _ = np.load(fn)
58 X, Y, elevs = _["X"], _["Y"], _["elevs"]
59 self._known_elevs[bn] = (X, Y, elevs)
60
61 pressureSensor.enable()
62 grav_sensor.enable()
63 gps.configure(on_location=self.on_location,
64 on_status=self.on_status)
65 gps.start()
66
67 def on_location(self, **kwargs):
68 lon = self.lon = float(kwargs['lon'])
69 lat = self.lat = float(kwargs['lat'])
70 self.alt = float(kwargs['altitude'])
71 now = datetime.now()
72 p = pressureSensor.value
73 if p and (now - self._last_known["time"]).total_seconds() > 10 * 60:
74
75 for bn in self._known_elevs:
76 X, Y, elevs = self._known_elevs[bn]
77 found_lat = np.where(Y <= lat)
78 found_lon = np.where(X <= lon)
79 if not len(found_lat[0]) or not len(found_lon[0]):
80 continue
81 i = found_lat[0][-1]
82 j = found_lon[0][-1]
83 elev = elevs[i, j]
84 self._last_known = {
85 "time": now,
86 "altitude": elev,
87 "pressure": p,
88 "lon": lon,
89 "lat": lat
90 }
91 #Logger.info(str(self._last_known))
92 break
93
94 @property
95 def current(self):
96 h_A = self._last_known["altitude"]
97 p_A = self._last_known["pressure"]
98 p = pressureSensor.value
99 return {
100 "time": (datetime.utcnow() - _ERA).total_seconds(),
101 "lon": self.lon,
102 "lat": self.lat,
103 "GPS altitude": self.alt,
104 "delta_T": self.delta_T,
105 "pressure": p,
106 "pressure altitude":
107 p2a(p, h_A, p_A, delta_T=self.delta_T),
108 "pressure altitude without correction":
109 p2a(p, 0, 1013.25, delta_T=0),
110 "last known altitude": h_A,
111 "last known pressure": p_A,
112 "last known lon": self._last_known["lon"],
113 "last known lat": self._last_known["lat"],
114 "gravity_X": grav_sensor.values[0],
115 "gravity_Y": grav_sensor.values[1],
116 "gravity_Z": grav_sensor.values[2],
117 "gravity": math.sqrt(
118 grav_sensor.values[0]**2 + grav_sensor.values[1]**2 + grav_sensor.values[2]**2),
119 }
120
121 def on_status(self, stype, status):
122 pass
本題とは無関係のデータも取ってるのは、ここらへんでやってるのと同じく重力にもやっぱし興味があって記録しておきたいけれども、2つサービス動かすのもやなので、てノリ。
よさげかな、と思い出かけた。ここ:
どこ? うん、ここ:
てわけで枡形山頂上は 84m、気圧高度の計算では 65.2m になっている。うーん、いいんだか悪いんだかわからんなぁ、と、「やや良いともやや悪いとも」の感触を持っていたのも束の間、ΔT (平均海水面気温15℃からのズレ)をちょっと変えてみるか、とした次の瞬間こんなことに:
うぞーん。ナニが起こったの??
数分は頭がハングアップしてたが、理解出来た。あぁ…しまった…、GPS 位置情報がとんでもない場所指すことへの配慮、やっぱなきゃダメだったか、と。
DEM データは結構狭い範囲でピンポイントでしか持ってなかったので、自分的には川崎駅での補正値が最後のもの、というつもりだったのね。それがまさかそれらにヒットしてしまうほどぶっ飛ぶとは。しかもそんな狭い範囲に器用にあたったもんだな…。帰ってからみておこう。
と思って調べてみたら実際は違った。あぅ…ロジック間違えてら…:
ダメじゃん…
1 X, Y, elevs = self._known_elevs[bn]
2 found_lat = np.where(Y <= lat)
3 found_lon = np.where(X <= lon)
4 if not len(found_lat[0]) or not len(found_lon[0]):
5 continue
6 i = found_lat[0][-1]
7 j = found_lon[0][-1]
8 elev = elevs[i, j]
無論 GPS 位置情報がとんでもない場所を指す問題は別途措置がいるけれども。うーん、もっと愚直に書くか…自分でわかんなくなってる:
スマートでなくても…
1 X, Y, elevs = self._known_elevs[bn]
2 if not (
3 X[0] <= lon and lon <= X[-1] and \
4 Y[0] <= lat and lat <= Y[-1]):
5 continue
6
7 i = np.where(Y <= lat)[0][-1]
8 j = np.where(X <= lon)[0][-1]
9 elev = elevs[i, j]
後日再チャレンジだな。つーかもっと遠出して高いとこに行きたいの。
2021-05-07追記:
紹介した basemap だが、作者である Jeff Whitaker が「basemap プロジェクト」としての保守を完全にやめ、pip でのインストールも不可能になった。自身による説明:
Deprecation Notice
Basemap is deprecated in favor of the Cartopy project. See notes in Cartopy, New Management, and EoL Announcement for more details.
リンク先をちょっと読むに、「より大きなプロジェクトに取り込まれた」ということに見える。cartopy はワタシ自身にとって完全に未知なのでまだ何も言わないけれど、余力があったらそのうち何か書くかも。
