曇天デバック 兼 DEM5 持ち歩きとそれで気圧補正したるぜ、のお試しが成功しなかったハナシ

思いついたのですぐに試してみたくなり。

国土地理院の基盤地図情報を持ち歩けないか、それで気圧補正すれば気圧高度計が結構精度よくなるんじゃないのか、のお試しをどうしてもすぐにやりたくなる。

発想のおさらい:

  1. スマホの気圧センサーは結構安定して値が取れる(大きく暴れることがない)
  2. やってみたらかなり細かい地形の起伏に繊細に反応してることがわかった
  3. 標準大気モデルに従った気圧高度を計算してみて、国土地理院の基盤地図情報標高データ(DEM5は5mメッシュ)と合わせてみると、かなり精密に起伏には反応していた
  4. 市販の気圧高度計のように「高度がわかっている地点で気圧補正」と気温補正をすると、結構一致した。
  5. 基盤地図情報標高データは「バカでかい」のでスマホに全部入れて持ってくなんて論外だが、「高度がわかっている地点で気圧補正」のためにピンポイントで抽出してスマホに突っ込んどくのはアリ。
  6. 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]

後日再チャレンジだな。つーかもっと遠出して高いとこに行きたいの。