「誤魔化して書かずにいたこと」についての追記は一部前のにもしたが、もう一つ。
「MSM_GPV_Rjp_L-pall」なのに pressure (つまり「level」)を意識せんでどーすんねん、という「気象に関する基礎」を疎かにしてしまったことはまぁどうでもいいとして、本題の「cartopyの」で意図的に話から抜いたのが「resolution」の話。チラ見えてたし、Basemap の話をした際にはした話なので、その対応関係からはしといた方がいい話であることは自覚してた。ただ、正確に伝えられるならそうした方がいいだろうなぁ、と思ったので、後回しにしたの。
やってみたら拍子抜けするほど単純な話だったので、前回のに入れときゃよかったとも思った。「サンプル」としてはこんだけね:
1 # -*- coding: utf-8 -*-
2 import os
3 import sys
4 import pygrib
5 import numpy as np
6 import cartopy.crs as ccrs
7 import cartopy.feature as cfeature
8 import matplotlib.pyplot as plt
9
10
11 def _read(srcfn):
12 src = pygrib.open(srcfn)
13 for s in src:
14 yield s.level, (s.name, s.parameterName), s.validDate, s.latlons(), s.values
15 src.close()
16
17
18 def main(gpv):
19 key = "Temperature"
20 data = {}
21 inivd, inilev = None, None
22 for lev, yk, vd, latlons, values in _read(gpv):
23 if key in yk and key not in data:
24 data[key] = latlons, values
25 inivd, inilev = vd, lev
26 if len(data) == 1:
27 break
28 (lats, lons), values = data[key]
29 X, Y = lons[0,:], lats[:,0]
30 #
31 for resol in ("110m", "50m", "10m",):
32 fig = plt.figure()
33 fig.set_size_inches(16.53, 11.69)
34 crs = ccrs.Geostationary(X.mean())
35 ax1 = fig.add_subplot(projection=crs)
36 clat, clon = 35.224167, 139.025373
37 ax1.set_extent([clon - 5, clon + 5, clat - 4, clat + 4], ccrs.PlateCarree())
38 CS = ax1.pcolormesh(
39 X, Y, values - 273.15, transform=ccrs.PlateCarree(), shading="auto")
40 ax1.coastlines(resolution=resol)
41 ax1.gridlines()
42 plt.colorbar(CS, ax=ax1)
43 #ax1.add_feature(cfeature.OCEAN)
44 #ax1.add_feature(cfeature.LAND)
45
46 plt.savefig("{}_{}_{}_p{}_res{}.png".format(
47 os.path.basename(gpv), str(inivd).partition(":")[0],
48 key, inilev, resol), bbox_inches="tight")
49 plt.close(fig)
50 #plt.show()
51
52
53 if __name__ == '__main__':
54 # 入力は "Z__C_RJTD_20210514000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin" など。
55 main(sys.argv[1])
basemap では「resolution=’f’ (full) / resolution=’c’ (crude)」だので制御し、これが実行速度に大きく影響し、「試行錯誤の最初のうちから f を使うな」みたいな注意喚起が必要だったのよ、これが cartopy ではどうなったのか、てことをね、言わねばなぁ、と思ってたわけだけれど。
感覚的に、上の「110m、50m、10m」での体感的な実行速度の差はないように思う。というか basemap の「f」が極端にヘビィだったのと較べると、全然速い。ただ、まず解像度はこの3つしか選べないのだけれど、最高解像度の「10m」が果たして望む最高解像度なんだろうか、てのは、まだちょっと良くわかんない。basemap の full よりいいんだろうか、悪いんだろうか。上の結果はこうだよ:
この範囲での「110m」が論外の解像度であることだけはわかるだろう。
basemap の「f」ほどに極端に重いってことがないんであれば、割と無頓着に「常に10mで」って考えも悪くなさそうだな、とは思った。
2021-05-31追記:
cartopy.io.shapereader お試しがてら、国土地理院配布の shapefile と cartopy の海岸線描画比較が出来たので…:
Click here to download
青が国土地理院の「coastl_jpn」によるもの、黒が cartopy の coastlines の「resolution=”10m”」。全般には国土地理院のものの方が解像度が高いようだが、たとえば北海道のサロマ湖などのような表現は cartopy のものの方が細かいようで、一概にどちらか一方だけが優れているってもんでもなさそう。(一応国土地理院のものは目的別に、例えば基盤地図情報を使って自分で海岸線ポリゴンを作り出すのも不可能ではないことも含め、「最高精度の可能性」を目指す場合は国土地理院から情報取得する一択だけれど、今は「coastl_jpn」についての話だけをしている。)
