GPVで沸点予報(誰得?)

ECMWF の GRIB-API がWindows 対応の実験版を提供し始めたので「喜んで」。

といっても、「Windows 対応の話」はおいおいで。やりました、できました、以上のことはここでは言わない。簡単といえば簡単だし、難しいといえば難しいのはbasemapのそれと同じく。それなりのスキルがあれば全く難しいところはないが、そうでないなら「破滅的に難し」く感じるかもしれない程度には難解、かも。GRIB-API そのものよりもね、pygribね、ハマるとすれば。

と、話をかなり唐突にはじめてしまいました、すいません。

  1. 「気象庁が提供してくれている、数値予報データ(GPV)」
  2. GPV は GRIB 形式
  3. ECMWF の GRIB-API は、GRIB を扱うインフラの中でも最も「清く正しい」もの
  4. つい最近までこれの Windows ビルドはほぼ不可能だった(Cygwin は可、MinGW では部分的に可)
  5. pygrib は GRIB-API の Python ラッパー
  6. 「気象庁が提供してくれている、数値予報データ(GPV)」を正式に入手するには、負担金が必要だがこれは本格的なサービス向け
  7. いくつかの研究機関が、「研究目的に」アーカイブを提供してくれている。その一つが京都大学

で、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、てわけです。