mpl_toolkits.basemap 開発はストップ、cartopy 使えやヲラ、てことなので

かつてのように pip install basemap 的なことをしても、もはや PyPI にも登録されてないわけよ、インストール出来ない、のな。

いつまで掲載され続けるかわからんけど、「deprecated になったぜのアナウンス」。deprecated どころか install 出来ないんで、お急ぎたい人にはお急ぎたい話になりうろう。

どういうやりとりがあっての流れなのかは良くわからないけれど、とにかく Basemap が 2016 年より「新規プロジェクトがこれを置き換える前提で」元の作者の Jeff Whitaker から「新規プロジェクト」の管理下に移されて、そして今まさに「置き換え完了した」てことなわけね。その「新規プロジェクト」てのが Cartopy project。ちょっとオモロイのは、Jeff Whitaker って、NOAA、つまり「アメリカ海洋気象庁」の人なんだけれども、Cartopy を起こした集団てこれ、イギリスの気象庁に強く関係してそうなんだよね。(プロジェクトそのものはイギリスの気象庁とは独立してるみたいだけども。)

こいつの導入は例によって「Windows でだけ鬱陶しい」。やってないけど Unix 系(Mac も)で苦労する要素はおそらくないと思う。細かく言うとこれの必須の依存物は:

  1. GEOS
  2. PROJ
    • PROJ 6.x 以降の場合は、これが sqlite3、curl、tiff に依存。
    • PROJ 8 は、説明されてないがまだ対応されてないので、7までを使うこと。

ワタシはこれまで PROJ 6.x 以前を使って何か苦労した記憶は、GEOS ともどもなかった。これは Windows 版ビルドですら同じで、何か問題があっても数秒で解決出来るような問題にしか巡り合ってきてない。PROJ 6.x 以上で sqlite3、curl、tiff が依存として増えたところで、Unix ならなんら苦労は増えないであろう。

PyPI に置かれてる公式の cartopy が「ビルド上等」、つまり、Windows ユーザが喜ぶタイプのバイナリ配布でないので、何も考えたくない人は、Christoph Gohlke の Unofficial Windows Binaries for Python Extension Packages のものを使うと良い:

1 [me@host: ~]$ py -3 -m pip install Cartopy-0.19.0.post1-cp39-cp39-win_amd64.whl

「非公式だとは何事だ」と気分悪さを感じる人や、あるいは何か事情があって正確にバージョン等々のコントロールが必要で、「むしろビルド上等上等!」という人のために、一応ワタシが成功したパターンを書いとく(当たり前だけど Visual Studio、Cython は自明の前提)。概要:

  1. GEOS は結構無頓着でもオッケーで、野良ビルドても、vcpkg のでもヨシ。
  2. 現時点最新の cartopy の setup.py の問題で、vcpkg の PROJ は使えない。proj.exe が必要なため。一応 PR を書いたが、受け容れられるかはわからん。
  3. PROJ について、共有ライブラリとしてのリンクがついに成功せず。
    • PROJ のビルド自体は出来るよ、そうじゃなくて。
    • cartopy の pyd (_crs.cp39-win_amd64.pyd 等)が proj.dll にリンクするが、このロードに失敗し、実行(インポート)出来ない。

そういうわけで、GEOS は vcpkg で取り寄せたものを、PROJ については、依存物のない PROJ 5系をダウンロードしてビルドすることにしよう、と。

c:/Program Files (x86)/Microsoft Visual Studio/2019/Community/VC/Auxiliary/Build/vcvars64.bat からスタートした(「管理者として実行」として起動した)コンソール(MSYS の bash)より:

CMake はインストールして、なおかつパスを通しておくのだぞ。
 1 [me@host: ~]$ tar zxvf proj-5.2.0.tar.gz
 2 [me@host: ~]$ cd proj-5.2.0
 3 [me@host: proj-5.2.0]$ mkdir build
 4 [me@host: proj-5.2.0]$ cd build
 5 [me@host: build]$ # 以下は「静的ライブラリとして作るぜ」。
 6 [me@host: build]$ cmake -DBUILD_SHARED_LIBS=OFF ..
 7    ... (snip) ...
 8 [me@host: build]$ # 以降は本来は cmake だけで実行できるはずのことなのだが、
 9 [me@host: build]$ # proj のバージョンごとに cmake のコントロールがマチマチで、
10 [me@host: build]$ # 情報通りにやっても意図したことにならないので、「cmake 箱庭
11 [me@host: build]$ # の外で遊ぶ」:
12 [me@host: build]$ export PATH="${PATH}":"/c/Program Files (x86)/Microsoft Visual Studio/2019/Community/MSBuild/Current/Bin/amd64/"
13 [me@host: build]$ MSBuild.exe -p:Configuration=Release -p:platform=x64 PROJ4.sln
14    ... (snip) ...
15 [me@host: build]$ cmake -P cmake_install.cmake
16    ... (c:/OSGeo4W にインストールされる、と思う) ...

バージョニングされた命名で c:/OSGeo4W/local/lib/proj_5_2.lib が出来上がるのだが、のちのリンクでこれが面倒で、幸いこれが static Library であるのだからその未来永劫名乗るべき名前には無頓着でいいわけで、なので、リンクするなりリネームするなりして c:/OSGeo4W/local/lib/proj.lib にしとくといい:

1 [me@host: build]$ ln c:/OSGeo4W/local/lib/proj_5_2.lib c:/OSGeo4W/local/lib/proj.lib

ワタシは vcpkg を「c:/develop」の下においたので、vcpkg 管理版の geos と今作った PROJ を参照してもらうようにする環境変数設定は例えばこうなる:

 1 [me@host: ~]$ # geos_c.dll が /bin の下にいるのであって、そこへパスを通す必要がある
 2 [me@host: ~]$ # のであって、geos-config.exe みたいなのがあってそれを起動できるように
 3 [me@host: ~]$ # するためではないよ。
 4 [me@host: ~]$ export PATH=/c/develop/vcpkg/packages/geos_x64-windows/bin
 5 [me@host: ~]$ # 対してこちら↓は proj.exe の存在を cartopy の setup.py が前提にするため。
 6 [me@host: ~]$ export PATH=/c/OSGeo4W/bin:$PATH
 7 [me@host: ~]$
 8 [me@host: ~]$ # LIB、INCLUDE は MSVC が参照してくれる環境変数ね。
 9 [me@host: ~]$ export LIB=c:/develop/vcpkg/packages/geos_x64-windows/lib";${LIB}"
10 [me@host: ~]$ export INCLUDE=c:/develop/vcpkg/packages/geos_x64-windows/include";${INCLUDE}"
11 [me@host: ~]$ export LIB=c:/OSGeo4W/local/lib";${LIB}"
12 [me@host: ~]$ export INCLUDE=c:/OSGeo4W/local/include";${INCLUDE}"

この状態にまで出来れば、きっと以下は問題なく成功する:

1 [me@host: ~]$ py -3 -m pip install cartopy

なお、「geos-config.exe なんかない」件に関して、スクリプト呼び出しの EXE 化を前提にして自分でまかなうことは出来るので、そちらが好みならそちらでどうぞ。その場合は geos-config.sh は例えばこんな感じ:

geos-config.sh
 1 #!/bin/sh
 2 
 3 # escape paths
 4 escape() {
 5   echo "$1" | sed 's/ /\\ /g'
 6 }
 7 
 8 prefix=`escape "c:/develop/vcpkg/packages/geos_x64-windows/"`
 9 # exec_prefix is interpolated in c:/develop/GEOS/lib
10 exec_prefix=`escape "c:/develop/vcpkg/packages/geos_x64-windows/bin"`
11 libdir=`escape "c:/develop/vcpkg/packages/geos_x64-windows/lib"`
12 
13 usage()
14 {
15   cat <<EOF
16 Usage: geos-config [OPTIONS]
17 Options:
18      [--prefix]
19      [--version]
20      [--libs]
21      [--clibs]
22      [--cclibs]
23      [--static-clibs]
24      [--static-cclibs]
25      [--cflags]
26      [--ldflags]
27      [--includes]
28      [--jtsport]
29 EOF
30   exit $1
31 }
32 
33 if test $# -eq 0; then
34   usage 1 1>&2
35 fi
36 
37 while test $# -gt 0; do
38   case "$1" in
39     -*=*) optarg=`echo "$1" | sed 's/[-_a-zA-Z0-9]*=//'` ;;
40     *) optarg= ;;
41   esac
42   case $1 in
43     --prefix)
44       echo ${prefix}
45       ;;
46     --version)
47       echo 3.9.1
48       ;;
49     ##--libs)
50     ##  # TODO: make an alias for --clibs
51     ##  # see http://trac.osgeo.org/geos/ticket/497
52     ##  echo -L${libdir} geos-3.9.1.lib
53     ##  ;;
54     --clibs)
55       echo -L${libdir} -lgeos_c
56       ;;
57     ##--cclibs)
58     ##  echo -L${libdir} geos.lib
59     ##  ;;
60     ##--static-clibs)
61     ##  echo -L${libdir} geos_c.lib geos.lib
62     ##  ;;
63     ##--static-cclibs)
64     ##  echo -L${libdir} geos.lib
65     ##  ;;
66     --cflags)
67       echo -I${prefix}/include
68       ;;
69     --ldflags)
70       echo -L${libdir}
71       ;;
72     --includes)
73       echo ${prefix}/include
74       ;;
75     ##--jtsport)
76     ##  echo 1.17.0
77     ##  ;;
78     *)
79       usage 1 1>&2
80       ;;
81   esac
82   shift
83 done

ともあれこれで無事導入が成功した、として。

2021-05-22追記:
すまん、色々無頓着で、雑に書いちゃダメな部分で雑にやっちゃった。ほかのものを検証中に二つ気づいたので、補足しておく。

一つは、「PROJ 5.x でビルド」は、公式にはこれは NG の可能性がある。cartopy ではないが pyproj の方ははっきりと最低レベルを PROJ 6.2 としていた。ビルドはできるし問題なく動いてはいるけれど、そういうのもあるので、Christoph Gohlke の方に頼っておいた方が無難かもしれない。彼のほうがワタシよりこうしたことへの措置には慣れているはずだから、同じ非公式でもワタシの決断よりはきっとあてになる、はず。

もう一つは Python 3.5 の話。2.7 と同じく、Python 3.5 も昨年既にそのライフサイクルを終えていて、なのでそもそもが「Python 3.5 をまだ使っている」ということ自体が不健全であることは念押ししつつも、個人が自由にやるんでなく組織として使う場合、乗り換えは簡単にはいかなかったりもすると思うので一応言っとく。ワタシがここに書いた情報だけで 3.5 版のビルドは出来ない。幸い Christoph Gohlke はまだ 3.5 版も保守し続けている。そちらを素直に使うといい。(3.4 も同じ話で、Christoph Gohlke 版は、本日時点で 2.7、3.4~3.9 を保守している。)どうしても 3.5 の自力ビルドを(今からわざわざ準備して)行う必要がある特殊な事情持ちの人は、こちらを参照のこと。この準備をすれば「たぶん可能」。

3.6 以降は、ワタシは実際に試して確認はしていないけれど、おそらく自力ビルドも大丈夫だろうと思う。(lru_dict では 3.5 が NG、3.6 以降大丈夫なのを実際に確認した、ので、それから考えればおそらく。)

さっそくお試したい、かんたんなへろうわルドはないものか、てのは、はっきりいってソースコードを持ってきて、その中にぶち込まれてる examples がまぁまぁ量があるので、それをお試せばよいかと思う:

1 [me@host: ~]$ py -3 cartopy/examples/gridlines_and_labels/gridliner.py

とか。そうなんだけれど、一応ワタシが最初に書いたヤツも見せとく:

 1 # -*- coding: utf-8 -*-
 2 import cartopy.crs as ccrs
 3 import cartopy.feature as cfeature
 4 import matplotlib.pyplot as plt
 5 
 6 
 7 def main():
 8     lat, lon = 35.224167, 139.025373  # 箱根駒ケ岳
 9     fig = plt.figure()
10     ax1 = fig.add_subplot(projection=ccrs.Geostationary(lon))
11     ax1.set_extent([lon - 2, lon + 2, lat - 2, lat + 2], ccrs.PlateCarree())
12     ax1.coastlines()
13     ax1.add_feature(cfeature.LAND)
14     #ax1.add_feature(cfeature.OCEAN)
15     ax1.gridlines()
16     # あとは普通に pclor だの countour だの普段どおりお好きに出来る。
17     plt.savefig("my1.png", bbox_inches="tight")
18     #plt.show()
19 
20 
21 if __name__ == '__main__':
22     main()

こんなん出ますよな: my1

なお、一応何もしなくても問題なく動いてるんだけれど、本来 PROJ_LIB をセットしてあげなきゃいかんのじゃないかな、と言うのだけは個人的に謎。設定しなければならんのだとしたら、たぶんこうなんじゃないかと思う:

1 [me@host: ~]$ export PROJ_LIB=c:/OSGeo4W/share

けどやってもやらなくても現状動いちゃってるんで、これの正しさがわからん。


まだほんとに軽く眺めただけなのだけれど、おそらく Basemap よりも綺麗に matplotlib に統合されてる。Basemap はちょっと癖が強いところがあったんだけど、cartopy の方が「matplotlib な日常」そのままに使える感がより強いんではないかしら、という印象を受けた。まぁ使い込みだすと印象変わってくる可能性もないではないけれどもね。

あと、これは想像なんだけれど、「pyproj」も、cartopy に完全に取り込まれるんじゃないかと思う。pyproj も Jeff Whitaker だからね。ついでにいうと pygrib も近い運命を辿るかもなぁて気がする。さてどうなるかね?


2021-05-15追記:
「あとは普通に pclor だの countour だの普段どおりお好き」が「pcolor, contour」の間違いなのはご愛嬌としても、最低でも「そうしたことが確かに出来るんだろうな?」てことを言わない情報って不安になると思う。それをしなかったのはワタシの準備が整ってなかったから。昨年頭に Windows 7 を退役させて Windows 10 に乗り換えたあとの環境整備を、色々ずっと怠ってきてたツケがね、今頃。

元の Basemap、それに今度のこの cartopy が威力を発揮する一番手始めなものは、やっぱり気象分野にある。例によって pygrib ね。Windows でのそれについては、ここに追記を書いておいた。として準備をしておけば、めでたく「お湯が湧く温度を予言出来る」、てことだけど、まぁもちっとシンプルな例から始めようかなと:

vis_msmgpv_temp_example.py
 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.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 = None
22     for yk, vd, latlons, values in _read(gpv):
23         if key in yk and key not in data:
24             data[key] = latlons, values
25             inivd = vd
26         if len(data) == 1:
27             break
28     (lats, lons), values = data[key]
29     X, Y = lons[0,:], lats[:,0]
30     #
31     fig = plt.figure()
32     fig.set_size_inches(16.53, 11.69)
33     crs = ccrs.Geostationary(X.mean())
34     ax1 = fig.add_subplot(projection=crs)
35     ax1.set_extent([X.min(), X.max(), Y.min(), Y.max()], ccrs.PlateCarree())
36     ax1.pcolormesh(X, Y, values, transform=ccrs.PlateCarree(), shading="auto")
37     ax1.coastlines()
38     ax1.gridlines()
39     #ax1.add_feature(cfeature.LAND)
40     #ax1.add_feature(cfeature.OCEAN)
41 
42     plt.savefig("{}_{}_{}.png".format(
43         os.path.basename(gpv), str(inivd).partition(":")[0], key), bbox_inches="tight")
44     #plt.show()
45 
46 
47 if __name__ == '__main__':
48     # 入力は "Z__C_RJTD_20210514000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin" など。
49     main(sys.argv[1])

Z__C_RJTD_20210514000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_2021-05-14 00_Temperature
気温予測のデータってもともと地形が綺麗に見えちゃうんで「coastlinesしないと位置とデータのマッピングがワカラン」てもんではなくて、そういう意味で「cartopy、ステキっ」がやや霞んでしまうんだけれど、逆にいえば最初におためすものとしては「あ、確かに正しいね」がすぐにわかるので、例としては結構ふさわしいもんではあったりする。

colorbar:

vis_msmgpv_temp_example2.py
 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.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 = None
22     for yk, vd, latlons, values in _read(gpv):
23         if key in yk and key not in data:
24             data[key] = latlons, values
25             inivd = vd
26         if len(data) == 1:
27             break
28     (lats, lons), values = data[key]
29     X, Y = lons[0,:], lats[:,0]
30     #
31     fig = plt.figure()
32     fig.set_size_inches(16.53, 11.69)
33     crs = ccrs.Geostationary(X.mean())
34     ax1 = fig.add_subplot(projection=crs)
35     ax1.set_extent([X.min(), X.max(), Y.min(), Y.max()], ccrs.PlateCarree())
36     CS = ax1.pcolormesh(X, Y, values, transform=ccrs.PlateCarree(), shading="auto")
37     ax1.coastlines()
38     ax1.gridlines()
39     plt.colorbar(CS, ax=ax1)
40     #ax1.add_feature(cfeature.LAND)
41     #ax1.add_feature(cfeature.OCEAN)
42 
43     plt.savefig("{}_{}_{}_cb.png".format(
44         os.path.basename(gpv), str(inivd).partition(":")[0], key), bbox_inches="tight")
45     #plt.show()
46 
47 
48 if __name__ == '__main__':
49     # 入力は "Z__C_RJTD_20210514000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin" など。
50     main(sys.argv[1])

Z__C_RJTD_20210514000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_2021-05-14 00_Temperature_cb

contourf (or contour)のも一応:

vis_msmgpv_temp_example3.py
 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.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 = None
22     for yk, vd, latlons, values in _read(gpv):
23         if key in yk and key not in data:
24             data[key] = latlons, values
25             inivd = vd
26         if len(data) == 1:
27             break
28     (lats, lons), values = data[key]
29     X, Y = lons[0,:], lats[:,0]
30     #
31     fig = plt.figure()
32     fig.set_size_inches(16.53, 11.69)
33     crs = ccrs.Geostationary(X.mean())
34     ax1 = fig.add_subplot(projection=crs)
35     ax1.set_extent([X.min(), X.max(), Y.min(), Y.max()], ccrs.PlateCarree())
36     #CS = ax1.contour(X, Y, values, transform=ccrs.PlateCarree())
37     CS = ax1.contourf(X, Y, values, transform=ccrs.PlateCarree())
38     ax1.coastlines()
39     ax1.gridlines()
40     plt.colorbar(CS, ax=ax1)
41 
42     #contourfでは以下は邪魔だが contour ではどちらかはあるといいか?
43     #ax1.add_feature(cfeature.OCEAN)
44     #ax1.add_feature(cfeature.LAND)
45 
46     #plt.savefig("{}_{}_{}_cont.png".format(
47     #    os.path.basename(gpv), str(inivd).partition(":")[0], key), bbox_inches="tight")
48     plt.savefig("{}_{}_{}_contf.png".format(
49         os.path.basename(gpv), str(inivd).partition(":")[0], key), bbox_inches="tight")
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])

Z__C_RJTD_20210514000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_2021-05-14 00_Temperature_cont
Z__C_RJTD_20210514000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_2021-05-14 00_Temperature_contf


2021-05-15 18:30追記:
出来てないことが理由で誤魔化してることって、すぐにバレてたりするかいね? tick についてワタシは何も言ってない。うまくいってなかったからだよ。

わかってみると、ちょっとイヤな事態になってるなぁと思う。一応解決の方向に向かおうとはしてるらしいのではあるけれど、リリースされてるやつに取り込まれてない気がする。

問題と思う事象が2つあって、一つは「明示的に set_xticks, set_yticks」しないと絶対に ticklabel を打ってくんない」点。もう一つがリンク先が問題にしてる「Cannot handle non-rectangular coordinate」問題。要は「PlateCarree ならいいが Geostationary だとダメ」てこと、例えば。

set_xticks, set_yticks を必ず行わねばならぬ、が鬱陶しいのは、以下のようにサボると:

vis_msmgpv_temp_example4.py
 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 from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
 9 import matplotlib.pyplot as plt
10 
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_u, key_v = "U component of wind", "V component of wind"
21     data = {}
22     inivd = None
23     for yk, vd, latlons, values in _read(gpv):
24         if key_u in yk and key_u not in data:
25             data[key_u] = latlons, values
26             if inivd is None:
27                 inivd = vd
28         if key_v in yk and key_v not in data:
29             data[key_v] = latlons, values
30             if inivd is None:
31                 inivd = vd
32         if len(data) == 2:
33             break
34     (lats, lons), values_u = data[key_u]
35     _, values_v = data[key_v]
36     X, Y = lons[0,:], lats[:,0]
37     #
38     fig = plt.figure()
39     fig.set_size_inches(16.53, 11.69)
40     crs = ccrs.PlateCarree(central_longitude=X.mean())  #Geostationary(X.mean())はダメ
41     ax1 = fig.add_subplot(projection=crs)
42     ax1.set_xticks(np.arange(X.min(), X.max() + 1), crs=ccrs.PlateCarree())
43     ax1.set_yticks(np.arange(Y.min(), Y.max() + 1), crs=ccrs.PlateCarree())
44     lon_formatter = LongitudeFormatter(zero_direction_label=True)
45     lat_formatter = LatitudeFormatter()
46     ax1.xaxis.set_major_formatter(lon_formatter)
47     ax1.yaxis.set_major_formatter(lat_formatter)
48     ax1.set_extent([X.min(), X.max(), Y.min(), Y.max()], ccrs.PlateCarree())
49     CS = ax1.pcolormesh(X, Y, np.sqrt(values_u**2 + values_v**2), transform=ccrs.PlateCarree(), shading="auto")
50     ax1.coastlines()
51     ax1.gridlines()
52     plt.colorbar(CS, ax=ax1)
53     #ax1.add_feature(cfeature.OCEAN)
54     #ax1.add_feature(cfeature.LAND)
55 
56     plt.savefig("{}_{}_{}.png".format(
57         os.path.basename(gpv), str(inivd).partition(":")[0], key_u[:-2]), bbox_inches="tight")
58     #plt.show()
59 
60 
61 if __name__ == '__main__':
62     # 入力は "Z__C_RJTD_20210514000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin" など。
63     main(sys.argv[1])

こんなん出てまう:
Z__C_RJTD_20210514000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_2021-05-14 00_U component of wi
「47.4°N」という tick を打ちたいわけないんであって、なおかつせっかく gridlines が引いてるとこにこそ打ちたい、てことであろう、その tick る位置をチマチマ計算で頑張って求める必要がある、ってこった。

まぁまだ新しいプロジェクトなのだから、てことで、多少の不足は受け容れつつ、てことかなぁ。


2021-05-16追記:
「その tick る位置をチマチマ計算で頑張って求める必要がある」に過剰に反応されると弱りそうなので、「困難なことに不平を言ってるわけじゃない」ことを一応示しとく:

「たとえば」。唯一の解だと思わないこと。
1     ax1.set_xticks(
2         [v for v in np.arange(np.ceil(X.min()), np.floor(X.max()) + 1)
3          if int(v) % 5 == 0],
4         crs=ccrs.PlateCarree())
5     ax1.set_yticks(
6         [v for v in np.arange(np.ceil(Y.min()), np.floor(Y.max()) + 1)
7          if int(v) % 5 == 0],
8         crs=ccrs.PlateCarree())

Z__C_RJTD_20210514000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin_2021-05-14 00_U component of wind
よく出来たインフラほど、こうしたことは「不満がある場合にだけカスタマイズ」としたくて書かねばならぬコード、という扱いとなり、「必ず書かなければならないもの」とはならない。てことね。


2021-06-09追記:
ticks の件なのだが、gallaly でキッチリ説明されてた:

vis_msmgpv_temp_example4.py
 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 #####from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
 9 import matplotlib.pyplot as plt
10 
11 
12 def _read(srcfn):
13     src = pygrib.open(srcfn)
14     for s in src:
15         yield s.level, (s.name, s.parameterName), s.validDate, s.latlons(), s.values
16     src.close()
17 
18 
19 def main(gpv):
20     key_u, key_v = "U component of wind", "V component of wind"
21     data = {}
22     inivd = None
23     for level, yk, vd, latlons, values in _read(gpv):
24         if key_u in yk and key_u not in data:
25             data[key_u] = latlons, values
26             if inivd is None:
27                 inivd = vd
28         if key_v in yk and key_v not in data:
29             data[key_v] = latlons, values
30             if inivd is None:
31                 inivd = vd
32         if len(data) == 2:
33             break
34     (lats, lons), values_u = data[key_u]
35     _, values_v = data[key_v]
36     X, Y = lons[0,:], lats[:,0]
37     #
38     fig = plt.figure()
39     fig.set_size_inches(16.53, 11.69)
40     crs = ccrs.Geostationary(X.mean())
41     ax1 = fig.add_subplot(projection=crs)
42     ######↓なんと、これらの苦労がいらんてことなのであった…
43     ######ax1.set_xticks(np.arange(X.min(), X.max() + 1), crs=ccrs.PlateCarree())
44     ######ax1.set_yticks(np.arange(Y.min(), Y.max() + 1), crs=ccrs.PlateCarree())
45     ######lon_formatter = LongitudeFormatter(zero_direction_label=True)
46     ######lat_formatter = LatitudeFormatter()
47     ######ax1.xaxis.set_major_formatter(lon_formatter)
48     ######ax1.yaxis.set_major_formatter(lat_formatter)
49     ext = [X.min(), X.max(), Y.min(), Y.max()]
50     ax1.set_extent(ext, ccrs.PlateCarree())
51     CS = ax1.pcolormesh(
52         X, Y, np.sqrt(values_u**2 + values_v**2),
53         transform=ccrs.PlateCarree(), shading="auto")
54     ax1.coastlines()
55     ax1.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)
56     plt.colorbar(CS, ax=ax1)
57     #ax1.add_feature(cfeature.OCEAN)
58     #ax1.add_feature(cfeature.LAND)
59 
60     plt.savefig("{}_{}_{}_1000hPa_ticks.png".format(
61         os.path.basename(gpv),
62         str(inivd).partition(":")[0], "wind"), bbox_inches="tight")
63     #plt.show()
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_20210609000000_MSM_GPV_Rjp_L-pall_FH36-39_grib2.bin_2021-06-10 12_wind_1000hPa_ticks

一つ前のは dms の価値があまりみえてなかったので、ちょっとだけ変えたヤツ
 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 matplotlib.pyplot as plt
 8 
 9 
10 def _read(srcfn):
11     src = pygrib.open(srcfn)
12     for s in src:
13         yield s.level, (s.name, s.parameterName), s.validDate, s.latlons(), s.values
14     src.close()
15 
16 
17 def main(gpv):
18     key_u, key_v = "U component of wind", "V component of wind"
19     data = {}
20     inivd = None
21     for level, yk, vd, latlons, values in _read(gpv):
22         if key_u in yk and key_u not in data:
23             data[key_u] = latlons, values
24             if inivd is None:
25                 inivd = vd
26         if key_v in yk and key_v not in data:
27             data[key_v] = latlons, values
28             if inivd is None:
29                 inivd = vd
30         if len(data) == 2:
31             break
32     (lats, lons), values_u = data[key_u]
33     _, values_v = data[key_v]
34     X, Y = lons[0,:], lats[:,0]
35     #
36     fig = plt.figure()
37     fig.set_size_inches(16.53, 11.69)
38     crs = ccrs.Stereographic(X.mean())
39     ax1 = fig.add_subplot(projection=crs)
40     ext = [X.min() + 8, X.max() - 8, Y.min() + 8, Y.max() - 8]
41     ax1.set_extent(ext, ccrs.PlateCarree())
42     CS = ax1.pcolormesh(
43         X, Y, np.sqrt(values_u**2 + values_v**2),
44         transform=ccrs.PlateCarree(), shading="auto")
45     ax1.coastlines()
46     ##### 「Thanks to the dms keyword, minutes are used when appropriate to display
47     ##### fractions of degree.」。dms の制御さえあればカスタムフォーマッタは
48     ##### あんましいらんであろう、というノリなんだと思うわ。
49     ax1.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)
50     plt.colorbar(CS, ax=ax1)
51 
52     plt.savefig("{}_{}_{}_1000hPa_ticks2.png".format(
53         os.path.basename(gpv),
54         str(inivd).partition(":")[0], "wind"), bbox_inches="tight")
55     #plt.show()
56 
57 
58 if __name__ == '__main__':
59     # 入力は "Z__C_RJTD_20210514000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin" など。
60     main(sys.argv[1])

Z__C_RJTD_20210609000000_MSM_GPV_Rjp_L-pall_FH36-39_grib2.bin_2021-06-10 12_wind_1000hPa_ticks2
上でリンクした issue はこの措置をしたということか…?



Related Posts