ECMWF の GRIB-API がWindows 対応の実験版を提供し始めたので「喜んで」。
といっても、「Windows 対応の話」はおいおいで。やりました、できました、以上のことはここでは言わない。簡単といえば簡単だし、難しいといえば難しいのはbasemapのそれと同じく。それなりのスキルがあれば全く難しいところはないが、そうでないなら「破滅的に難し」く感じるかもしれない程度には難解、かも。GRIB-API そのものよりもね、pygribね、ハマるとすれば。
と、話をかなり唐突にはじめてしまいました、すいません。
- 「気象庁が提供してくれている、数値予報データ(GPV)」
- GPV は GRIB 形式
- ECMWF の GRIB-API は、GRIB を扱うインフラの中でも最も「清く正しい」もの
- つい最近までこれの Windows ビルドはほぼ不可能だった(Cygwin は可、MinGW では部分的に可)
- pygrib は GRIB-API の Python ラッパー
- 「気象庁が提供してくれている、数値予報データ(GPV)」を正式に入手するには、負担金が必要だがこれは本格的なサービス向け
- いくつかの研究機関が、「研究目的に」アーカイブを提供してくれている。その一つが京都大学
で、GPV の面白い遊び方なんかいっくらでもあるんだけれども、数ある遊びの中でも最も価値がなさそうな、「沸点予報」してみようと思った。なんで、って聞かないで。やりたくなっただけ。
「飽和蒸気圧が大気圧と等しくなる温度で液体は沸騰する」ということは…、GPV から得られる「Surface pressure」を飽和水蒸気圧として、その飽和水蒸気圧となる温度を求めれば良い、はず。Murray の逆関数を使えばよいか。
ひとまずちょっと基礎を知らない人のために、「Surface pressure」と「Pressure reduced to MSL」を可視化したものを先に見せておく：
1 # -*- coding: utf-8 -*-
2 import sys
3 import pygrib
4 import numpy as np
5
6
7 def _read(srcfn):
8 # keys in Surfece:
9 # Pressure reduced to MSL
10 # Surface pressure
11 # 10 metre U wind component
12 # 10 metre V wind component
13 # Temperature
14 # Relative humidity
15 # Low cloud cover
16 # Medium cloud cover
17 # High cloud cover
18 # Total cloud cover
19 # Total precipitation
20 #
21 src = pygrib.open(srcfn)
22 for s in src:
23 yield (s.name, s.parameterName), s.validDate, s.latlons(), s.values
24 src.close()
25
26
27 data = {}
28 for yk, vd, latlons, values in _read(sys.argv[1]):
29 if "Surface pressure" in yk and "Surface pressure" not in data:
30 data["Surface pressure"] = latlons, values
31 elif "Pressure reduced to MSL" in yk and "Pressure reduced to MSL" not in data:
32 data["Pressure reduced to MSL"] = latlons, values
33 if len(data) == 2:
34 break
35
36
37 # -----------------------------------------------------
38 #
39 from mpl_toolkits.basemap import Basemap
40 import matplotlib.pyplot as plt
41 import matplotlib.cm as cm
42 import matplotlib.ticker as ticker
43
44
45 bmparams = dict(
46 llcrnrlat=34.25,
47 urcrnrlat=37.5,
48 llcrnrlon=137.,
49 urcrnrlon=140.825,
50 resolution='i',
51 projection='merc')
52
53 fig, axes = plt.subplots(1, 2)
54 fig.set_size_inches(16.53 * 2, 11.69 * 2)
55
56 circles = np.arange(10, 90, 5)
57 meridians = np.arange(-180, 180, 5)
58 # -----------------------------------------------------
59 latlons, values = data["Surface pressure"]
60 lons = latlons[1]
61 lats = latlons[0]
62 #
63 X = latlons[1][0,:]
64 Y = latlons[0][:,0]
65 # -----------------------------------------------------
66 bmparams.update(dict(
67 llcrnrlat=Y[-1],
68 urcrnrlat=Y[0],
69 llcrnrlon=X[0],
70 urcrnrlon=X[-1]))
71
72 bmparams["ax"] = axes[0]
73 m = Basemap(**bmparams)
74 x, y = m(lons, lats)
75 m.drawcoastlines(color='white', linewidth=1)
76 m.drawparallels(circles, labels=[1, 0, 0, 0])
77 m.drawmeridians(meridians, labels=[0, 0, 0, 1])
78
79 #
80 cs = m.contourf(
81 x, y,
82 values / 100., # Pa to hPa
83 50,
84 cmap=cm.coolwarm_r)
85 m.colorbar(
86 cs,
87 format=ticker.FuncFormatter(lambda x, pos: "%.1f" % x))
88 axes[0].set_title("Surface pressure\n[" + str(vd) + "]")
89 # -----------------------------------------------------
90 latlons, values = data["Pressure reduced to MSL"]
91 #
92 bmparams["ax"] = axes[1]
93 m = Basemap(**bmparams)
94 x, y = m(lons, lats)
95 m.drawcoastlines(color='black', linewidth=2)
96 m.drawparallels(circles, labels=[1, 0, 0, 0])
97 m.drawmeridians(meridians, labels=[0, 0, 0, 1])
98
99 #
100 cs = m.contourf(
101 x, y,
102 values / 100., # Pa to hPa
103 20,
104 cmap=cm.coolwarm_r)
105 m.colorbar(
106 cs,
107 format=ticker.FuncFormatter(lambda x, pos: "%.1f" % x))
108 axes[1].set_title("Pressure reduced to MSL\n[" + str(vd) + "]")
109 # -----------------------------------------------------
110 #plt.show()
111 plt.savefig("pressure.png", bbox_inches="tight")
「Surface pressure」が文字通りその場所の気圧、なので、高い山ではいつでも「低気圧」になる。「Pressure reduced to MSL」がこれを海水面レベルに補正したもので、天気予報でみる気圧配置は常にこちらである。つまり前者は天気予報の役には全く立たない。けれども、現地で実測したとしたら、実測値と合うのは前者の方で、今回の目的に使えるのも前者。
さて。沸点温度は、一年通して毎日同じような絵しか描けないはずである。高山がいつでも沸点温度が低く、平地はいつでも高い。「Surface pressure」と傾向は同じはずで、似た絵になるはずだ。これでは面白くない。少なくとも「いつもと較べてどうか」を、カラースケールをみなくてもなんとなくわかるようにしたい。ので、「沸点温度予報」の絵と「Pressure reduced to MSL」の絵を並べるのがいい。こうしておけば「普段よりもお湯が沸くのが早いはず」なのかそうでないのかが、視覚的にわかりやすいであろう：
1 # -*- coding: utf-8 -*-
2 import sys
3 import pygrib
4 import numpy as np
5
6
7 def _read(srcfn):
8 # keys in Surfece:
9 # Pressure reduced to MSL
10 # Surface pressure
11 # 10 metre U wind component
12 # 10 metre V wind component
13 # Temperature
14 # Relative humidity
15 # Low cloud cover
16 # Medium cloud cover
17 # High cloud cover
18 # Total cloud cover
19 # Total precipitation
20 #
21 src = pygrib.open(srcfn)
22 for s in src:
23 yield (s.name, s.parameterName), s.validDate, s.latlons(), s.values
24 src.close()
25
26
27 data = {}
28 for yk, vd, latlons, values in _read(sys.argv[1]):
29 if "Pressure reduced to MSL" in yk and "Pressure reduced to MSL" not in data:
30 data["Pressure reduced to MSL"] = latlons, values
31 elif "Surface pressure" in yk and "Surface pressure" not in data:
32 data["Surface pressure"] = latlons, values
33 if len(data) == 2:
34 break
35
36
37 # -----------------------------------------------------
38 #
39 from mpl_toolkits.basemap import Basemap
40 import matplotlib.pyplot as plt
41 import matplotlib.cm as cm
42 import matplotlib.ticker as ticker
43
44
45 bmparams = dict(
46 llcrnrlat=34.25,
47 urcrnrlat=37.5,
48 llcrnrlon=137.,
49 urcrnrlon=140.825,
50 resolution='i',
51 projection='merc')
52
53 fig, axes = plt.subplots(1, 2)
54 fig.set_size_inches(16.53 * 2, 11.69 * 2)
55
56 circles = np.arange(10, 90, 5)
57 meridians = np.arange(-180, 180, 5)
58 # -----------------------------------------------------
59 # surface pressure as saturated vapour pressure
60 latlons, e_s = data["Surface pressure"]
61
62 # inverse of Murray's equation
63 ln_es = np.log(e_s / 100. / 6.1078)
64 T_B = 273.15 + (-ln_es * 237.3) / (ln_es - 17.2693882)
65
66
67 lons = latlons[1]
68 lats = latlons[0]
69 #
70 X = latlons[1][0,:]
71 Y = latlons[0][:,0]
72 bmparams["ax"] = axes[0]
73 m = Basemap(**bmparams)
74 x, y = m(lons, lats)
75 m.drawcoastlines(color='white', linewidth=1)
76 m.drawparallels(circles, labels=[1, 0, 0, 0])
77 m.drawmeridians(meridians, labels=[0, 0, 0, 1])
78
79 cs = m.contourf(
80 x, y,
81 T_B - 273.15, # K to Celsius
82 20,
83 cmap=cm.coolwarm_r)
84 m.colorbar(cs)
85 axes[0].set_title("Boiling point of water\n[" + str(vd) + "]")
86 # -----------------------------------------------------
87 bmparams.update(dict(
88 llcrnrlat=Y[-1],
89 urcrnrlat=Y[0],
90 llcrnrlon=X[0],
91 urcrnrlon=X[-1]))
92 # -----------------------------------------------------
93 latlons, values = data["Pressure reduced to MSL"]
94 #
95 bmparams["ax"] = axes[1]
96 m = Basemap(**bmparams)
97 x, y = m(lons, lats)
98 m.drawcoastlines(color='black', linewidth=2)
99 m.drawparallels(circles, labels=[1, 0, 0, 0])
100 m.drawmeridians(meridians, labels=[0, 0, 0, 1])
101
102 #
103 cs = m.contourf(
104 x, y,
105 values / 100., # Pa to hPa
106 20,
107 cmap=cm.coolwarm_r)
108 m.colorbar(
109 cs,
110 format=ticker.FuncFormatter(lambda x, pos: "%.1f" % x))
111 axes[1].set_title("Pressure reduced to MSL\n[" + str(vd) + "]")
112 # -----------------------------------------------------
113 #plt.show()
114 plt.savefig("boiling_point_of_water.png", bbox_inches="tight")
今日は高気圧に覆われてるので富士山でも比較的沸点が高いのですな。
正確性は求めないでくださいな。ワタシは専門家ちゃうので、精度についての詳しいことは言えないけれど、最低でも「1気圧」時の計算は、小数点以下は合ってない。
なんのこっちゃ? はいそうですね、すいません。
いや、立体地形図遊びで basemap 紹介出来たのはいいんだけど、ちょっと「Basemapらしい」ことみせてあげたいなと思ってさ。標高データなんかは「地図に地図を重ねる」ようなもんであって、Basemap の価値がわかりにくいんだけれど、こういう「観測データ」とか「予測データ」とともに使ってこその basemap、てわけです。
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 はワタシ自身にとって完全に未知なのでまだ何も言わないけれど、余力があったらそのうち何か書くかも。
