本質的にはタイルそのものとは関係ないんだけれど、「典型的に」タイルを使う場合に悩みがち、なある事柄について。
(3)では、alpha を使って重ねる方法だけをやってみせたのだけれど、まぁ気になる人はすぐに気になると思うのだけど、「読みにくい」絵になりがちなわけね。
今回やってみせるこれは、本質的にはもう全然 cartopy と関係ないんだけれど、ただ、タイルを使わず海岸線だけ描画してくれるのを期待するだけの場合はこの悩みは発生しにくく、cartopy でタイルを使おうとして初めて悩む人もいるかもしらんので、一応紹介しとこうかと思った次第。
以下のように何も工夫しない K-Index ビジュアライズ:
1 # -*- coding: utf-8 -*-
2 import os
3 import sys
4 import datetime
5
6 import pygrib
7 import numpy as np
8 import matplotlib.pyplot as plt
9 import cartopy.crs as ccrs
10 from cartopy.io.img_tiles import GoogleTiles
11
12 # ...(省略)...
13
14
15 def main(gpv):
16 _alldata = {}
17 # ...(_alldata を gpv より。詳細は省略)...
18 #
19 # _alldata は、トップレベルのキーが MSM GPV データの日付時刻、
20 # 値がさらに辞書で、その辞書はキーが level (気圧面) で、値が
21 # 「GPVData」という自前 class。GPVData の中身は割愛。気象庁
22 # 提供のオリジナルデータとともに、それに基づくいくつかの計算結果
23 # を保持しているだけ。
24 X, Y = lons[0,:], lats[:,0]
25 tiler = GoogleTiles(style="satellite")
26 for vd in _alldata.keys():
27 fig = plt.figure()
28 fig.set_size_inches(16.53, 11.69)
29 ax1 = fig.add_subplot(projection=tiler.crs)
30 ext = [X.min(), X.max(), Y.min(), Y.max()]
31 ax1.add_image(tiler, 6)
32 ax1.set_extent(ext, ccrs.PlateCarree())
33
34 # _alldata[vd] を渡すと K-Index を計算する
35 # _K_index。これも省略。
36 george = _K_index(_alldata[vd])
37 CS = ax1.pcolormesh(
38 X, Y, george,
39 transform=ccrs.PlateCarree(),
40 shading="auto")
41 ax1.coastlines(resolution="10m")
42 ax1.gridlines()
43 plt.colorbar(CS)
44 ax1.set_title("K-index on {} JST".format(
45 str(vd + datetime.timedelta(hours=9))))
46 plt.savefig("{}_K-Index_{}_0.png".format(
47 os.path.basename(gpv),
48 str(vd).partition(":")[0]), bbox_inches="tight")
49 plt.close(fig)
50
51
52 if __name__ == '__main__':
53 # 入力は "Z__C_RJTD_20210514000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin" など。
54 main(sys.argv[1])
これは例えばこんな絵になる:
一つには「せっかくの GoogleTiles が台無しだよ参ったな」に目が行くと思うが、まぁ焦らない。「K-Index」のビジュアライズにおいて肝要なのは「値がどれだけなのか」という直値そのものではない。COVID-19 新規感染者数を、感染者数そのもので見るのではなく「ステージ」に分類するのに似た話。K-Index の場合は「20未満の値は「存在しないのと一緒」」「35超えるのがやっべーんだぜ」という色分け、これをして初めて K-Index のビジュアライズとして実用になる。ので、まずは例えばこうする:
1 # -*- coding: utf-8 -*-
2 import os
3 import sys
4 import datetime
5
6 import pygrib
7 import numpy as np
8 import matplotlib.pyplot as plt
9 from matplotlib.colors import Normalize
10 import cartopy.crs as ccrs
11 from cartopy.io.img_tiles import GoogleTiles
12
13 class MidpointNormalize(Normalize):
14 def __init__(self, vmin=None, vmax=None, vcenter=None, clip=False):
15 self.vcenter = vcenter
16 Normalize.__init__(self, vmin, vmax, clip)
17
18 def __call__(self, value, clip=None):
19 # I'm ignoring masked values and all kinds of edge cases to make a
20 # simple example...
21 x, y = [self.vmin, self.vcenter, self.vmax], [0, 0.5, 1]
22 return np.ma.masked_array(np.interp(value, x, y))
23
24 # ...(省略)...
25
26
27 def main(gpv):
28 _alldata = {}
29 # ...(_alldata を gpv より。詳細は省略)...
30 #
31 # _alldata は、トップレベルのキーが MSM GPV データの日付時刻、
32 # 値がさらに辞書で、その辞書はキーが level (気圧面) で、値が
33 # 「GPVData」という自前 class。GPVData の中身は割愛。気象庁
34 # 提供のオリジナルデータとともに、それに基づくいくつかの計算結果
35 # を保持しているだけ。
36 X, Y = lons[0,:], lats[:,0]
37 norm = MidpointNormalize(vcenter=30.0, vmin=20.0)
38 tiler = GoogleTiles(style="satellite")
39 for vd in _alldata.keys():
40 fig = plt.figure()
41 fig.set_size_inches(16.53, 11.69)
42 ax1 = fig.add_subplot(projection=tiler.crs)
43 ext = [X.min(), X.max(), Y.min(), Y.max()]
44 ax1.add_image(tiler, 6)
45 ax1.set_extent(ext, ccrs.PlateCarree())
46
47 # _alldata[vd] を渡すと K-Index を計算する
48 # _K_index。これも省略。
49 george = _K_index(_alldata[vd])
50 CS = ax1.pcolormesh(
51 X, Y, george,
52 norm=norm,
53 transform=ccrs.PlateCarree(),
54 shading="auto")
55 ax1.coastlines(resolution="10m")
56 ax1.gridlines()
57 plt.colorbar(CS)
58 ax1.set_title("K-index on {} JST".format(
59 str(vd + datetime.timedelta(hours=9))))
60 plt.savefig("{}_K-Index_{}_1.png".format(
61 os.path.basename(gpv),
62 str(vd).partition(":")[0]), bbox_inches="tight")
63 plt.close(fig)
64
65
66 if __name__ == '__main__':
67 # 入力は "Z__C_RJTD_20210514000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin" など。
68 main(sys.argv[1])
こうなる:
MidpointNormalize は古い matplotlib であっても使えるが、新しい matplotlib を前提にするなら、matplotlib.colors.TwoSlopeNorm や matplotlib.colors.CenteredNorm も使えるので、使えるならそちらを使えばいい。
GoogleTiles 台無し問題を度外視すれば、まぁこれでもいいことはいいのだが、「せっかくなので」とこうするのが…:
1 # -*- coding: utf-8 -*-
2 import os
3 import sys
4 import datetime
5
6 import pygrib
7 import numpy as np
8 import matplotlib.pyplot as plt
9 from matplotlib.colors import Normalize
10 import cartopy.crs as ccrs
11 from cartopy.io.img_tiles import GoogleTiles
12
13 class MidpointNormalize(Normalize):
14 def __init__(self, vmin=None, vmax=None, vcenter=None, clip=False):
15 self.vcenter = vcenter
16 Normalize.__init__(self, vmin, vmax, clip)
17
18 def __call__(self, value, clip=None):
19 # I'm ignoring masked values and all kinds of edge cases to make a
20 # simple example...
21 x, y = [self.vmin, self.vcenter, self.vmax], [0, 0.5, 1]
22 return np.ma.masked_array(np.interp(value, x, y))
23
24 # ...(省略)...
25
26
27 def main(gpv):
28 _alldata = {}
29 # ...(_alldata を gpv より。詳細は省略)...
30 #
31 # _alldata は、トップレベルのキーが MSM GPV データの日付時刻、
32 # 値がさらに辞書で、その辞書はキーが level (気圧面) で、値が
33 # 「GPVData」という自前 class。GPVData の中身は割愛。気象庁
34 # 提供のオリジナルデータとともに、それに基づくいくつかの計算結果
35 # を保持しているだけ。
36 X, Y = lons[0,:], lats[:,0]
37 norm = MidpointNormalize(vcenter=30.0, vmin=20.0)
38 tiler = GoogleTiles(style="satellite")
39 for vd in _alldata.keys():
40 fig = plt.figure()
41 fig.set_size_inches(16.53, 11.69)
42 ax1 = fig.add_subplot(projection=tiler.crs)
43 ext = [X.min(), X.max(), Y.min(), Y.max()]
44 ax1.add_image(tiler, 6)
45 ax1.set_extent(ext, ccrs.PlateCarree())
46
47 # _alldata[vd] を渡すと K-Index を計算する
48 # _K_index。これも省略。
49 george = _K_index(_alldata[vd])
50 CS = ax1.pcolormesh(
51 X, Y, george,
52 norm=norm,
53 transform=ccrs.PlateCarree(),
54 shading="auto",
55 alpha=0.5)
56 ax1.coastlines(resolution="10m")
57 ax1.gridlines()
58 plt.colorbar(CS)
59 ax1.set_title("K-index on {} JST".format(
60 str(vd + datetime.timedelta(hours=9))))
61 plt.savefig("{}_K-Index_{}_2.png".format(
62 os.path.basename(gpv),
63 str(vd).partition(":")[0]), bbox_inches="tight")
64 plt.close(fig)
65
66
67 if __name__ == '__main__':
68 # 入力は "Z__C_RJTD_20210514000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin" など。
69 main(sys.argv[1])
これはやだなぁてことよ:
かっこわるい、ということだけでなくて、そもそも非常に読み取りにくい。
こうした場合に何を考えればいいかと言えば、問題を「全ての値域に対して色を使おうとしてることが悪いのだ」とすり替えてしまうと良い。事実今の例の場合、K-Index が 20 未満のどぎつい紫色が美しい GoogleTiles の衛星画像を覆い隠してしまっている、ということなのだから。つまり:
1 # -*- coding: utf-8 -*-
2 import os
3 import sys
4 import datetime
5
6 import pygrib
7 import numpy as np
8 import matplotlib.pyplot as plt
9 from matplotlib.colors import Normalize
10 import cartopy.crs as ccrs
11 from cartopy.io.img_tiles import GoogleTiles
12
13 class MidpointNormalize(Normalize):
14 def __init__(self, vmin=None, vmax=None, vcenter=None, clip=False):
15 self.vcenter = vcenter
16 Normalize.__init__(self, vmin, vmax, clip)
17
18 def __call__(self, value, clip=None):
19 # I'm ignoring masked values and all kinds of edge cases to make a
20 # simple example...
21 x, y = [self.vmin, self.vcenter, self.vmax], [0, 0.5, 1]
22 return np.ma.masked_array(np.interp(value, x, y))
23
24 # ...(省略)...
25
26
27 def main(gpv):
28 _alldata = {}
29 # ...(_alldata を gpv より。詳細は省略)...
30 #
31 # _alldata は、トップレベルのキーが MSM GPV データの日付時刻、
32 # 値がさらに辞書で、その辞書はキーが level (気圧面) で、値が
33 # 「GPVData」という自前 class。GPVData の中身は割愛。気象庁
34 # 提供のオリジナルデータとともに、それに基づくいくつかの計算結果
35 # を保持しているだけ。
36 X, Y = lons[0,:], lats[:,0]
37 norm = MidpointNormalize(vcenter=30.0, vmin=20.0)
38 tiler = GoogleTiles(style="satellite")
39 for vd in _alldata.keys():
40 fig = plt.figure()
41 fig.set_size_inches(16.53, 11.69)
42 ax1 = fig.add_subplot(projection=tiler.crs)
43 ext = [X.min(), X.max(), Y.min(), Y.max()]
44 ax1.add_image(tiler, 6)
45 ax1.set_extent(ext, ccrs.PlateCarree())
46
47 # _alldata[vd] を渡すと K-Index を計算する
48 # _K_index。これも省略。
49 george = _K_index(_alldata[vd])
50 # ↓20以上にあらずんば人にあらず。
51 george[np.where(george < 20)] = np.nan
52
53 CS = ax1.pcolormesh(
54 X, Y, george,
55 norm=norm,
56 transform=ccrs.PlateCarree(),
57 shading="auto")
58 ax1.coastlines(resolution="10m")
59 ax1.gridlines()
60 plt.colorbar(CS)
61 ax1.set_title("K-index on {} JST".format(
62 str(vd + datetime.timedelta(hours=9))))
63 plt.savefig("{}_K-Index_{}_3.png".format(
64 os.path.basename(gpv),
65 str(vd).partition(":")[0]), bbox_inches="tight")
66 plt.close(fig)
67
68
69 if __name__ == '__main__':
70 # 入力は "Z__C_RJTD_20210514000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin" など。
71 main(sys.argv[1])
正直ちゃんと考えて設計された結果こうなってるのかについてはわからんのだけれど、結果論としては NaN を colormap は「値域外につき色なし」として扱った結果が、タイルとの重ね合わせ時には透ける:
これが唯一の解かどうかはわからないが、多くのケースでは、これだけの措置で十分と思ってる。(NaN になる領域の広さによっては結局 alpha も駆使することにはなると思うけれど、NaN を活用する方がより穏健な alpha でいい、となる、はず。)
この措置なしの場合 cmap をあれやこれや変えても視覚的効果としての改善にはまったく寄与しないが、この措置をした状態なら例えば:
1 # ...
2 CS = ax1.pcolormesh(
3 X, Y, george,
4 norm=norm,
5 transform=ccrs.PlateCarree(),
6 shading="auto",
7 cmap="gist_heat")
8 # ...
1 # ...
2 CS = ax1.pcolormesh(
3 X, Y, george,
4 norm=norm,
5 transform=ccrs.PlateCarree(),
6 shading="auto",
7 cmap="hsv_r")
8 # ...
1 # ...
2 CS = ax1.pcolormesh(
3 X, Y, george,
4 norm=norm,
5 transform=ccrs.PlateCarree(),
6 shading="auto",
7 cmap="jet")
8 # ...
みたいに望みのものになる、はず。
ところで、「GoogleTiles(style='terrain')」がそもそも実用にならんのだが…:
1 # -*- coding: utf-8 -*-
2 import os
3 import sys
4 import datetime
5
6 import pygrib
7 import numpy as np
8 import matplotlib.pyplot as plt
9 import cartopy.crs as ccrs
10 from cartopy.io.img_tiles import GoogleTiles
11
12 def _read(srcfn):
13 src = pygrib.open(srcfn)
14 for s in src:
15 yield (s.name, s.parameterName), s.validDate, s.latlons(), s.values
16 src.close()
17
18
19 def main(gpv):
20 key = "Total precipitation"
21 _alldata = {}
22 (lats, lons) = None, None
23 for yk, vd, latlons, values in _read(gpv):
24 (lats, lons) = latlons
25 if key in yk:
26 _alldata[vd] = values
27 X, Y = lons[0,:], lats[:,0]
28 tiler = GoogleTiles(style='terrain')
29 for vd in _alldata.keys():
30 fig = plt.figure()
31 fig.set_size_inches(16.53, 11.69)
32 ax1 = fig.add_subplot(projection=tiler.crs)
33 ext = [X.min(), X.max(), Y.min(), Y.max()]
34 ax1.add_image(tiler, 6)
35 ax1.set_extent(ext, ccrs.PlateCarree())
36 cmap = "turbo"
37 v = _alldata[vd]
38 v[np.where(v == 0)] = np.nan
39 # pcolormesh しないのも試してみるといい。そもそも
40 # 「terrain なタイルそのもの」からして機能してない。
41 CS = ax1.pcolormesh(
42 X, Y, v,
43 transform=ccrs.PlateCarree(),
44 shading="auto",
45 cmap=cmap)
46 ax1.coastlines(resolution="10m")
47 ax1.gridlines()
48 plt.colorbar(CS)
49 ax1.set_title("{} on {} JST".format(
50 key, str(vd + datetime.timedelta(hours=9))))
51 plt.savefig("{}_{}_0_{}.png".format(
52 os.path.basename(gpv), key,
53 str(vd).partition(":")[0], cmap), bbox_inches="tight")
54 plt.close(fig)
55
56
57 if __name__ == '__main__':
58 # 入力は "Z__C_RJTD_20210520030000_MSM_GPV_Rjp_Lsurf_FH34-39_grib2.bin" など。
59 main(sys.argv[1])
色々試行錯誤してみたが、全然ダメ。なんだろうなぁこれ。
ところで、「ついでなので」。一つ前の「ダメでした」で雨量を例にしたのは実は意味があって。FAQ に近いネタなんだけれど実はこの雨量の可視化、見かけよりむつかしい。
天気予報がいとも簡単にやってのけてるのでなかなか気付かないのだが、是非天気予報で採用しているカラースケールの目盛りをじっくり観察してみてほしい。どういう目盛りになってる? テレビの場合だと、各社少しずつ違ってはいるものの、等幅ではないことがわかる。
なんでそうしてるかは、多分人間側の感度にある程度対応付けたいんだと思う。「微量の降水」と「微々量の降水」をどの程度区別したいか、「強烈な降水」と「激烈な降水」をどの程度区別したいか、みたいなことね。どういうのが日常生活で便利か、という試行錯誤なんだろうと思う。
この件は実は COVID-19 関連可視化ネタで触れたことがあるし、もちろん上でやった「MidpointNormalize」もその一つである。matplotlib.colors.Normalize ファミリで実現することが出来る。
雨量で使いやすいのは matplotlib.colors.TwoSlopeNorm かもしれない:
1 # -*- coding: utf-8 -*-
2 import os
3 import sys
4 import datetime
5
6 import pygrib
7 import numpy as np
8 import matplotlib.pyplot as plt
9 from matplotlib.colors import TwoSlopeNorm
10 import cartopy.crs as ccrs
11 from cartopy.io.img_tiles import GoogleTiles
12
13
14 def _read(srcfn):
15 src = pygrib.open(srcfn)
16 for s in src:
17 yield (s.name, s.parameterName), s.validDate, s.latlons(), s.values
18 src.close()
19
20
21 def main(gpv):
22 key = "Total precipitation"
23 _alldata = {}
24 (lats, lons) = None, None
25 for yk, vd, latlons, values in _read(gpv):
26 (lats, lons) = latlons
27 if key in yk:
28 _alldata[vd] = values
29 X, Y = lons[0,:], lats[:,0]
30 norm = TwoSlopeNorm(vcenter=5)
31 tiler = GoogleTiles(style='satellite')
32 for vd in _alldata.keys():
33 fig = plt.figure()
34 fig.set_size_inches(16.53, 11.69)
35 ax1 = fig.add_subplot(projection=tiler.crs)
36 ext = [X.min() + 1, X.max() - 1, Y.min() + 3, Y.max() - 3]
37 ax1.add_image(tiler, 6)
38 ax1.set_extent(ext, ccrs.PlateCarree())
39 cmap = "turbo"
40 v = _alldata[vd]
41 v[np.where(np.isclose(v, 0))] = np.nan
42 CS = ax1.pcolormesh(
43 X, Y, v,
44 norm=norm,
45 transform=ccrs.PlateCarree(),
46 shading="auto",
47 cmap=cmap)
48 ax1.coastlines(resolution="10m")
49 ax1.gridlines()
50 plt.colorbar(CS)
51 ax1.set_title("{} on {} JST".format(
52 key, str(vd + datetime.timedelta(hours=9))))
53 plt.savefig("{}_{}_1_{}.png".format(
54 os.path.basename(gpv), key,
55 str(vd).partition(":")[0], cmap), bbox_inches="tight")
56 plt.close(fig)
57
58
59 if __name__ == '__main__':
60 # 入力は "Z__C_RJTD_20210520030000_MSM_GPV_Rjp_Lsurf_FH34-39_grib2.bin" など。
61 main(sys.argv[1])
「弱い雨の強弱にはあまり興味がない」とかあるいはその逆、みたいなことで vcenter を調整すればいい。たとえば同じ入力に対してTwoSlopeNorm(vcenter=1)とすればこうなる:
TwoSlopeNorm(vcenter=10)とすればこうなる:
完全にテレビの天気予報と同じにするのだと、BoundaryNorm と ListedColormap を組み合わせることで可能。NHK 方式のであればこう:
1 # -*- coding: utf-8 -*-
2 import os
3 import sys
4 import datetime
5
6 import pygrib
7 import numpy as np
8 import matplotlib.pyplot as plt
9 from matplotlib.colors import BoundaryNorm, ListedColormap
10 import cartopy.crs as ccrs
11 from cartopy.io.img_tiles import GoogleTiles
12
13
14 def _read(srcfn):
15 src = pygrib.open(srcfn)
16 for s in src:
17 yield (s.name, s.parameterName), s.validDate, s.latlons(), s.values
18 src.close()
19
20
21 def main(gpv):
22 key = "Total precipitation"
23 _alldata = {}
24 (lats, lons) = None, None
25 for yk, vd, latlons, values in _read(gpv):
26 (lats, lons) = latlons
27 if key in yk:
28 _alldata[vd] = values
29 X, Y = lons[0,:], lats[:,0]
30 #
31 # NHK方式の [0, 1, 5, 10, 20, 30, 50, 80]
32 cmap = ListedColormap(
33 ['white', 'cyan', 'skyblue', 'blue', 'yellow', 'orange', 'red', 'magenta'])
34 norm = BoundaryNorm(
35 [0, 1, 5, 10, 20, 30, 50, 80], cmap.N, extend="max")
36 #
37 tiler = GoogleTiles(style='satellite')
38 for vd in _alldata.keys():
39 fig = plt.figure()
40 fig.set_size_inches(16.53, 11.69)
41 ax1 = fig.add_subplot(projection=tiler.crs)
42 ext = [X.min(), X.max(), Y.min(), Y.max()]
43 ax1.add_image(tiler, 6)
44 ax1.set_extent(ext, ccrs.PlateCarree())
45 v = _alldata[vd]
46 v[np.where(np.isclose(v, 0))] = np.nan
47 CS = ax1.pcolormesh(
48 X, Y, v,
49 norm=norm,
50 transform=ccrs.PlateCarree(),
51 shading="auto",
52 cmap=cmap)
53 ax1.coastlines(resolution="10m")
54 ax1.gridlines()
55 plt.colorbar(CS)
56 ax1.set_title("{} on {} JST".format(
57 key, str(vd + datetime.timedelta(hours=9))))
58 plt.savefig("{}_{}_1_{}_bn.png".format(
59 os.path.basename(gpv), key,
60 str(vd).partition(":")[0], cmap), bbox_inches="tight")
61 plt.close(fig)
62
63
64 if __name__ == '__main__':
65 # 入力は "Z__C_RJTD_20210520030000_MSM_GPV_Rjp_Lsurf_FH34-39_grib2.bin" など。
66 main(sys.argv[1])
以下の絵は上のスクリプトそのものの結果ではなくて、九州付近をクリップしたもの:
2022-03-26 AM2時頃追記:
明日雨が降るということだが、その前に行動すれば降られないかしら、と地上波デジタルの予報を見ようとしたんだけれど、どの局も、三時間ごと予報は明日の夜までの予報なのに、詳細な雨量予測は6時間先までしか見せてくれないんだよね、明日の10時くらいから昼過ぎくらいまでの状況が知りたいのに。
というわけで、じゃぁ GPV のデータを可視化しちゃるかと:
上でやったスクリプトから変更してるのは「alpha=0.7」で可視化してること、だけ。そう、v[np.where(np.isclose(v, 0))] = np.nan と「穏健な alpha」ね、上でちょろっと言ったそれを実際にやった。cartopy ネタとしての追記の「テクニカルな部分の本題」はそれだけ。
というかね、「明日のワタシの行動」のために、実用のために、とこの絵を読むと悩ましくてな。まず、地上波デジタルの予報(おそらく気象庁やウェザーニュースのサイトを直接訪問しても同じハズだが)と、この GPV 可視化が微妙に違っていて。地上波デジタルの予報だと、大田区で明日降り始めるのは15時頃から、てことになってるわけね。けどこの絵、10時のものよ、降る予測になってる。うーむ。
そもそも「雨量≦1mm ってどんくらい?」てのを、ワタシはこれまで意識してきたことがなかったんだよな。当然そういうのを書いてくれる人は必ずいる。これとかこれとか。そうか、1mm以下って「傘がいるか悩むレベル」なのか。つまり、この「白」は、あまり行動に支障がない程度の降水。なるほど。でも15時でもこんなよ:
うーん、これは v[np.where(np.isclose(v, 0))] = np.nan ではなくて、もう少し大きな値で切り捨てたほうがいいかも。四捨五入のノリで v[np.where(v < 0.5)] = np.nan としたとするなら:
これだとやりすぎか。v[np.where(v < 0.1)] = np.nan ならどうだ?:
1mm以下の線状帯は10時のからは消えたが、11時～14時くらいには出たまま。そして15時にはないのよね…。うん、どうしたもんか。明日いつ行動すれば降られないで済むか、「天気予報」に基づけば15時より前だが、この GPV での予測の可視化を信じるなら「11～14時に行動するくらいなら15時のほうがマシ」てなことに。む。まぁ「10時くらいまでの用事を済ませる」つもりならそれが良さそうな気もするけどね。いや…というか、明日のワタシの行動範囲内だと、10～15時の間はずっと「傘がいらんのではないか」という程度の雨量でしかない、てことか、あ、それでいいならそれでいいわ。あんまし悩む必要はなかったか、GPV 可視化に基づく方は。ので問題は「天気予報」で提供されてる雨量予測との乖離だけな。6mm言うてたのよ、ほかの天気予報だと15時以降。GPV 予測はそうは言ってないんだよねぇ。まぁどっちにしても、明日は早めに行動を済ませよ、てことで FA か、ワタシのケースでは。
実際発表されてる気象情報って、この GPV 生予測にどんな補正をかけて伝えてるんだろうなぁ? その正確なことが知りたいよね。どっかに書かれてたりするんだろか??
2022-03-26追記(2)(「2022-03-26 AM2時頃追記」の「顛末」＋α):
結局 GPV が正しかった。11時過ぎに外に出て、まさに「傘がいらないくらいの降雨」でしばらく、して、12時ちょっと前くらいに、数分だけ通り雨的に「ちゃんと」降った。ちょっとだけ民家の駐車場の屋根を借りて雨宿りしつつ傘を出し、少しだけ傘をさして歩くことになった。
それはそうと、こういうのをごくごくたまーにやりたくなるんだけれど、それをするたびに毎度ちょっと気持ち悪いと思ってて、それでも調べてこなかったことを、やっとちゃんと調べようかなと。京都大学のアーカイブについてのメモね。ただしワタシは MSM-GPV しか使わないので、それについてだけ。
ちゃんと知っとかないといけないなと思ったのは「最新の予測」のありようについて。どういうことかと言うと、われわれ「京都大学のアーカイブをあてにする立場」からみて、たぶん気象庁から購入するのとは違うタイムラグがあるんだろうなぁ、と思うから。きっとこういうことだと思うの:
- 京都大学が気象庁より生データを受け取る
- 受け取ったデータを何かしら加工とかしてるかもしんない(してないかもしんない)
- データをアーカイブサイトに配置する
- 配置したデータを可視にする(要するにインデクスページを更新する)
この「1.」の前の話と、「1.」以降の話について、ちゃんと整理して理解しないといけないなと。そういうことね。
たとえば「Z__C_RJTD_20211208000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin」の「20211208000000」部分と「FH00-15」部分。後者はこれは「0時間後～15時間後予測」という意味。前者が「初期値」と呼ばれている値で、これはシミュレーションに与えた初期値の時刻なのだが、ワタシの理解だと「観測データを反映したもの」。で、ワタシが今回知らないといけないなと思ったのは、この「初期値～を、ワタシがいつ手に出来るのか」ね。
キャプチャした画像にみえている「Last Modified」は、まずは JST で、で、おそらく「気象庁から京都大学サーバがデータを受け取った時刻」を反映しているんだと思う。これを 2009 年のものから本日のものまでざっと見てみたところ、どれも「初期値の時刻 + 2時間15分くらい」らしい。日によって多少前後するが、どうやらそのくらい。で、この時間差というのが何を意味するのかだが、これはむろん「わからない」が正解なんだけれど、少なくとも「気象庁のスーパーコンピュータが観測データを初期値として入力して、それとモデルに基づいて予測演算する時間」があるはずで、そこから京都大学のサーバに至るまでの時間があるわけね、それが平均で2時間15分ということなのか、あるいは単に京都大学のサーバのアーカイブのためのスケジューラの都合なのかもしれない。これは本当に気象庁からデータを購入してみないとわからない。そこらへんの事情は書いてくれてないから。
ここまでは、過去のものだけ眺めてればわかることなのだが、これまでちゃんと見てこなかったのは、先の列挙の「4. 配置したデータを可視にする」部分なのだわ。実際ずっと貼りついて更新ボタンし続けてみるまでわからないと思うんだけれど、上にはりつけた画像の「2022-03-26 12:32」とあるデータは、これは「2022-03-26 12:32 JST」には我々のもとにはどうやら届かない。MSM GPV の Surface の場合 3時間ごとの配信があるんだけど、どうも本日の様子をみてると、「2022-03-26 12:32 JST の 1～3時間後」にようやくダウンロード可能になるらしい、つまり丁度2周回遅れくらいみたい。なるほど。(というかインデクスページの更新がそうなだけで、実体は既に存在していて、ファイル名がわかっているなら直接ダウンロード出来るのかもしれんけど、「絶賛生成中」とならない時間を選んでるんだろうから、まぁインデクスページが更新されるのを待つのが無難だろうね。)
つまり、「京都大学のアーカイブを利用する場合は、初期値時刻からみた 5～6 時間後までの予測は、ほぼ過去についての予測になる」てことかなと。FH00-15 が15時間後までの予測なので、それでも 06-15 の9時間分は「予測として役立つ」てことね。
まぁずっとそうだろうなと理解していたことを、ちゃんと一日観察してみた、てだけのことなの。直接購入するとこのタイムラグがどの程度なくなるんだろうねぇ? まぁ我々のような一般人ユーザにとっては、この6時間は許容範囲かなと思ってる。
2022-03-27追記:
昨日の二つ目の追記は「cartopy ネタ」として機能してなかったわね、すまん。ただ、今回の追記も、元の「tileたい人々への補足事項的な」という本題からは外れるので、きっと「またしてもすまん」。まぁでも cartopy ネタだよ。
昨晩の整理の後も少し cartopy 遊びを続けてたのね、ちょっと久しぶりだったから。その中で、昨年このシリーズを書いた際には書かなかったネタがあったことに気付いたので、ここに追記しとく。
「1mm以下の降水量」の件で、「0.1 で切るか 0.5 で切るか」みたいなのを「二つの絵」として別々に作るのか、subplots の rows, cols を駆使して並べるのか、の話があるでしょ。一連の cartopy ネタの中で、このパターンはそういえばやってなかったんだよね。個人的には複雑になるので「「二つの絵」として別々に作る」のを優先することをお勧めしたいけれど、でもどうしてもやりたいことはあるんだよね。大した話ではないんだけど、なにせ matplotlib だけで存分に複雑なので、例は多いほうが良い、と。たとえばこうね:
1 # -*- coding: utf-8 -*-
2 import os
3 import sys
4 import datetime
5
6 import pygrib
7 import numpy as np
8 import matplotlib.pyplot as plt
9 from matplotlib.colors import BoundaryNorm, ListedColormap
10 import cartopy.crs as ccrs
11 from cartopy.io.img_tiles import GoogleTiles
12
13
14 _today = datetime.datetime.utcnow()
15
16
17 def _read(srcfn):
18 src = pygrib.open(srcfn)
19 for s in src:
20 #if _today < s.validDate:
21 yield (s.name, s.parameterName), s.validDate, s.latlons(), s.values
22 src.close()
23
24
25 def main(gpv):
26 key = "Total precipitation"
27 _alldata = {}
28 (lats, lons) = None, None
29 for yk, vd, latlons, values in _read(gpv):
30 (lats, lons) = latlons
31 if key in yk:
32 print(vd)
33 _alldata[vd] = values
34 if lats is None:
35 sys.exit(0)
36 X, Y = lons[0,:], lats[:,0]
37 #
38 # NHK方式の [0, 1, 5, 10, 20, 30, 50, 80]
39 cmap = ListedColormap(
40 ['white', 'cyan', 'skyblue', 'blue', 'yellow', 'orange', 'red', 'magenta'])
41 norms = []
42 clip_min = (0.1, 0.5)
43 for clmin in clip_min:
44 norms.append(BoundaryNorm(
45 [clmin, 1, 5, 10, 20, 30, 50, 80], cmap.N, extend="max"))
46 #
47 tocho = (139.691717, 35.689568)
48 #
49 tiler = GoogleTiles(style='satellite')
50 for vd in _alldata.keys():
51 outname = "{}_{}_1_{}_bn.png".format(
52 os.path.basename(gpv), key,
53 str(vd).partition(":")[0])
54 if os.path.exists(outname):
55 continue
56 fig = plt.figure()
57 fig.set_size_inches(16.53 * 1.5, 11.69)
58 axes = []
59 axes.append(fig.add_subplot(1, 2, 1, projection=tiler.crs))
60 axes.append(fig.add_subplot(1, 2, 2, projection=tiler.crs))
61 #ext = [X.min(), X.max(), Y.min(), Y.max()]
62 ext = [tocho[0] - 6, tocho[0] + 1, tocho[1] - 3, tocho[1] + 3]
63 for ax in axes:
64 ax.add_image(tiler, 6)
65 ax.set_extent(ext, ccrs.PlateCarree())
66 CS = []
67 for i, clmin in enumerate(clip_min):
68 v = _alldata[vd].copy()
69 v[np.where(v < clmin)] = np.nan
70 CS.append(axes[i].pcolormesh(
71 X, Y, v,
72 norm=norms[i],
73 transform=ccrs.PlateCarree(),
74 shading="auto",
75 cmap=cmap, alpha=0.7))
76 for i, ax in enumerate(axes):
77 ax.coastlines(resolution="10m")
78 ax.gridlines()
79 plt.colorbar(CS[i], ax=axes[i])
80 axes[i].set_title("{} on {} JST\n(v[np.where(v < {})] = np.nan)".format(
81 key, str(vd + datetime.timedelta(hours=9)), clip_min[i]))
82 plt.savefig(outname, bbox_inches="tight")
83 plt.close(fig)
84
85
86 if __name__ == '__main__':
87 # 入力は "Z__C_RJTD_20210520030000_MSM_GPV_Rjp_Lsurf_FH34-39_grib2.bin" など。
88 main(sys.argv[1])
この例の場合は二つの Axes 両方を projection=tiler.crs で使うので Figure の add_subplot でなく pyplot.subplots も使えると思う。けど projection を軸ごとに変えたい場合、特に一方を cartopy とは関係ないグラフにしたい場合などもあるから Figure の add_subplot がいいと個人的には思う。こんな絵になる:
