せっかくTYPE_PRESSURE 取れてるので少しは遊んでみる。
気圧が計れていれば Pressure Altitude を計算出来る。ので「色んな高度」を出して遊んでみる。
国際標準大気モデル(ISA)に基いた Pressure Altitude の計算式を NOAA が公開してるんですと:
\displaystyle (1-(\frac{P_a}{1013.25})^{.190284})*145366.45
\)
ISA の以下を逆演算するつもりでいた:
\displaystyle P_a = P_0 \cdot \left(1 – \frac{L \cdot h}{T_0} \right)^\frac{g \cdot M}{R \cdot L}
\)
逆演算になってんのかな、どうなのかな。確かめてない。(P_0 は言うまでもなく 101325 [Pa]、T_0 は 15 + 273.15 [K] である。L とか R とかわからん人にはわからんと思うが、L が temperature lapse rate、R が universal gas constant である…、ってますますわからんか。)
さて。これとは別に、当然スマホには GPS センサーがついており、この子も高度を検出してくれる。これは幾何高度。もうひとつ、緯度経度がわかれば、ジオイド高がわかる。
ジオイド高のデータは世界標準のグリッドモデルがある。GeographicLib はこのグリッドデータへのアクセスが出来る API があるが、シンプルな WEB API がある。プログラムからはとても使いやすいというシロモノではなくて、こんな結果で返って来る:
1 <html>
2 <head>
3 <title>
4 Online geoid calculator
5 </title>
6 <meta name="description" content="Online geoid calculator" />
7 <!-- ... snip -->
8 <p>
9 Geoid height:
10 <font size="4"><pre>
11 lat lon = <a href="http://tools.wmflabs.org/geohack/geohack.php?params=35.55819;139.70072" style="color:Black">35.55819 139.70072</a> (35°33'29"N 139°42'03"E)
12 geoid heights (m)
13 <a href="http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008">EGM2008</a> = <font color="blue">35.9787</font>
14 <a href="http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html">EGM96</a> = <font color="blue">36.0662</font>
15 <a href="http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html">EGM84</a> = <font color="blue">37.1507</font></pre></font>
16 </p>
17 </form>
18 <!-- ... snip -->
インチキ EAN lookupと同じノリで HTMLParser で書くのがさらにダルいが単純な構造ではあるので、面倒なのでもう正規表現で:
1 # -*- coding: utf-8 -*-
2 import re
3 import urllib
4 import urllib2
5
6
7 _RGX = re.compile(
8 r'<a href="http://earth\-info[^"]+">(EGM[0-9]+)</a>\s*=\s*<[^<>]+>([0-9.]+)')
9 _LOC_RESOL_SECS = 1 # seconds
10 _CACHE = {}
11
12 def query_geoid_eval(lon, lat):
13 resol = 3600. / (_LOC_RESOL_SECS * 60 * 60.)
14
15 result = {}
16 url = "http://geographiclib.sourceforge.net/cgi-bin/GeoidEval"
17 data = urllib.urlencode({
18 "input": "%s %s" % (
19 "%.5f" % (int(lat * resol) / resol),
20 "%.5f" % (int(lon * resol) / resol)),
21 "option": "Submit",
22 })
23 if data in _CACHE:
24 return _CACHE[data]
25 opener = urllib2.build_opener()
26 #opener.addheaders = [('User-Agent', user_agent)]
27 res = opener.open(url + "?" + data)
28 lines = res.read().split("\n")
29 for line in lines:
30 m = _RGX.search(line)
31 if m:
32 model, elev = m.group(1, 2)
33 result[model] = elev
34 _CACHE[data] = result
35 return result
environmental_sensors.py は一つ前の投稿のものと同じ。
あとは:
1 title=Various Altitude Demo
2 author=hhsprings
3 orientation=all
1 # -*- coding: utf-8 -*-
2 import math
3 from kivy.lang import Builder
4 from kivy.app import App
5 from kivy.clock import Clock, mainthread
6 from kivy.properties import StringProperty
7 from plyer import gps
8 from environmental_sensors import pressureSensor
9 from geographiclib_web import query_geoid_eval
10
11 kv = '''
12 BoxLayout:
13 RstDocument:
14 id: result
15 text: app.result
16 '''
17 class VariousAltitudeDemo(App):
18 result = StringProperty()
19 lon = None
20 lat = None
21 alt = None
22
23 def build(self):
24 pressureSensor.enable()
25 gps.configure(on_location=self.on_location,
26 on_status=self.on_status)
27 gps.start()
28 Clock.schedule_interval(self.do_interval, 1.)
29 return Builder.load_string(kv)
30
31 @mainthread
32 def on_location(self, **kwargs):
33 self.lon = float(kwargs['lon'])
34 self.lat = float(kwargs['lat'])
35 self.alt = float(kwargs['altitude'])
36
37 @mainthread
38 def on_status(self, stype, status):
39 pass
40
41 def do_interval(self, *args, **kwargs):
42 geoid = {}
43 if self.lon:
44 geoid = query_geoid_eval(self.lon, self.lat)
45
46 P = pressureSensor.value
47
48 # NOAA formula (in feet)
49 # (see https://en.wikipedia.org/wiki/Pressure_altitude)
50 pres_alt = (1.0 - math.pow((P / 1013.25), 0.190284)) * 145366.45
51 pres_alt *= 0.3048 # in metre
52
53 self.result = """
54 GPS
55 ===
56 * lon: {} [deg]
57 * lat: {} [deg]
58 * alt: {} [m]
59
60 World Geoid Model
61 =================
62 * EGM2008: {} [m]
63 * EGM96: {} [m]
64 * EGM84: {} [m]
65
66 Pressure Altitude
67 =================
68 * Pressure: {} [hPa]
69 * Pressure altitude: {} [m]
70
71 """.format(self.lon, self.lat, self.alt,
72 geoid.get('EGM2008'), geoid.get('EGM96'), geoid.get('EGM84'),
73 P, pres_alt)
74
75
76 if __name__ == '__main__':
77 VariousAltitudeDemo().run()
こんな:
幾何高度と地理高度の説明はいいかな? GPS がはじき出す幾何高度は「地球楕円体だけに基いた幾何学的高度」、地理高度は「およそ地球上の生きとし生けるものは重力には逆らえないのである」ことに従って、重力に逆らわない計測をするもの、です、簡単に言えば。「ワタシの真上」(鉛直)は、地球上にいる限りは重力の向き、よね? これが幾何学的な地球中心を向いてない。この重力場のモデルが EGM。「平均海水面(MSL)」は重力モデルに従ったもの。そして「標高」とはこの MSL からの高さね。
ここで計算した Pressure Altitude はいわゆる補正なしの高度で、「気圧高度計」として使うには不十分。登山で気圧高度計を使う場合は、高度が地図からわかっている地点で高度を補正する。つまり「海水面気圧が1気圧、海水面気温が15℃」に従ったのがこの計算なんだけれど、当たり前だけどいつでもこの状態じゃぁない。けど GPS 高度もジオイド高もセットでこうやって一緒にみるといいよね。たとえば「今ほんとは100mにいるのに標高500m相当の気圧だ!」と面白がることが出来る。
このアプリをさぁ、飛行機の中で使えたら面白かったんだろうけどねぇ…。民間航空機の中って、思いっきり気圧調整してあるからさ、外気圧はわからんのよね。面白いんだけどなぁ、「民間航空機は(巡航時)気圧面に沿って飛んでいる」ことなんかが測れて面白かったんじゃないかと思うんだけどね、これは出来ない。
