mpl_toolkits.basemap 開発はストップ、cartopy 使えやヲラ、てことなので (4) – tileたい人々への補足事項的な

本質的にはタイルそのものとは関係ないんだけれど、「典型的に」タイルを使う場合に悩みがち、なある事柄について。

(3)では、alpha を使って重ねる方法だけをやってみせたのだけれど、まぁ気になる人はすぐに気になると思うのだけど、「読みにくい」絵になりがちなわけね。

今回やってみせるこれは、本質的にはもう全然 cartopy と関係ないんだけれど、ただ、タイルを使わず海岸線だけ描画してくれるのを期待するだけの場合はこの悩みは発生しにくく、cartopy でタイルを使おうとして初めて悩む人もいるかもしらんので、一応紹介しとこうかと思った次第。

以下のように何も工夫しない K-Index ビジュアライズ:

version 0
 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])

これは例えばこんな絵になる:
Z__C_RJTD_20210520000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_K-Index_2021-05-20 00_0

一つには「せっかくの GoogleTiles が台無しだよ参ったな」に目が行くと思うが、まぁ焦らない。「K-Index」のビジュアライズにおいて肝要なのは「値がどれだけなのか」という直値そのものではない。COVID-19 新規感染者数を、感染者数そのもので見るのではなく「ステージ」に分類するのに似た話。K-Index の場合は「20未満の値は「存在しないのと一緒」」「35超えるのがやっべーんだぜ」という色分け、これをして初めて K-Index のビジュアライズとして実用になる。ので、まずは例えばこうする:

version 1
 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])

こうなる:
Z__C_RJTD_20210520000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_K-Index_2021-05-20 00_1

MidpointNormalize は古い matplotlib であっても使えるが、新しい matplotlib を前提にするなら、matplotlib.colors.TwoSlopeNormmatplotlib.colors.CenteredNorm も使えるので、使えるならそちらを使えばいい。

GoogleTiles 台無し問題を度外視すれば、まぁこれでもいいことはいいのだが、「せっかくなので」とこうするのが…:

version 2 (不満)
 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])

これはやだなぁてことよ:
Z__C_RJTD_20210520000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_K-Index_2021-05-20 00_2
かっこわるい、ということだけでなくて、そもそも非常に読み取りにくい。

こうした場合に何を考えればいいかと言えば、問題を「全ての値域に対して色を使おうとしてることが悪いのだ」とすり替えてしまうと良い。事実今の例の場合、K-Index が 20 未満のどぎつい紫色が美しい GoogleTiles の衛星画像を覆い隠してしまっている、ということなのだから。つまり:

version 3
 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 は「値域外につき色なし」として扱った結果が、タイルとの重ね合わせ時には透ける:
Z__C_RJTD_20210520000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_K-Index_2021-05-20 00_3
Z__C_RJTD_20210520000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_K-Index_2021-05-20 03_3
Z__C_RJTD_20210520000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_K-Index_2021-05-20 06_3
Z__C_RJTD_20210520000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_K-Index_2021-05-20 09_3
Z__C_RJTD_20210520000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_K-Index_2021-05-20 12_3
Z__C_RJTD_20210520000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_K-Index_2021-05-20 15_3
これが唯一の解かどうかはわからないが、多くのケースでは、これだけの措置で十分と思ってる。(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 # ...

Z__C_RJTD_20210520000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_K-Index_2021-05-20 00_3_gist_heat

1 # ...
2         CS = ax1.pcolormesh(
3             X, Y, george,
4             norm=norm,
5             transform=ccrs.PlateCarree(),
6             shading="auto",
7             cmap="hsv_r")
8 # ...

Z__C_RJTD_20210520000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_K-Index_2021-05-20 00_3_hsv_r

1 # ...
2         CS = ax1.pcolormesh(
3             X, Y, george,
4             norm=norm,
5             transform=ccrs.PlateCarree(),
6             shading="auto",
7             cmap="jet")
8 # ...

Z__C_RJTD_20210520000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_K-Index_2021-05-20 00_3_jet
みたいに望みのものになる、はず。


ところで、「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])

Z__C_RJTD_20210520000000_MSM_GPV_Rjp_Lsurf_FH34-39_grib2.bin_Total precipitation_0_2021-05-21 09
色々試行錯誤してみたが、全然ダメ。なんだろうなぁこれ。


ところで、「ついでなので」。一つ前の「ダメでした」で雨量を例にしたのは実は意味があって。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])

Z__C_RJTD_20210520000000_MSM_GPV_Rjp_Lsurf_FH34-39_grib2.bin_Total precipitation_1_2021-05-21 09
「弱い雨の強弱にはあまり興味がない」とかあるいはその逆、みたいなことで vcenter を調整すればいい。たとえば同じ入力に対してTwoSlopeNorm(vcenter=1)とすればこうなる:
Z__C_RJTD_20210520000000_MSM_GPV_Rjp_Lsurf_FH34-39_grib2.bin_Total precipitation_1_2021-05-21 09vc1
TwoSlopeNorm(vcenter=10)とすればこうなる:
Z__C_RJTD_20210520000000_MSM_GPV_Rjp_Lsurf_FH34-39_grib2.bin_Total precipitation_1_2021-05-21 09vc10

完全にテレビの天気予報と同じにするのだと、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])

以下の絵は上のスクリプトそのものの結果ではなくて、九州付近をクリップしたもの:
Z__C_RJTD_20200704000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin_Total precipitation_1_2020-07-04 00_bn



Related Posts