プログラミング系科目オール「C」の俺がPythonに一目惚れして局所天気図の描写に挑む話 ～第1話～
あらすじ
大学の必修科目に設定されたプログラミング(C言語)の講義を機に、一気にアンチプログラミング人間になってしまった俺ことうえはーす。とある講義で利用したPythonの使い勝手の良さに一目惚れし、AMeDAS観測点データから局所天気図を描き、熱的低気圧と雷雨性高気圧を表せないかと試みる。一応は描けたように見えた天気図であったが、どうやら本人はその出来に納得いかないようで……
環境について
話を進める前に、万が一にもこの記事を読んで似たようなことをしてみたいという人が現れた時のために、Pythonの開発環境について述べておく。
実行するためのプロンプトとプログラムを構築する開発環境には、それぞれ大学の講義の一環でインストールしたAnacondaとVisual Code Studioを利用している。ダウンロード方法や設定方法については詳しいサイトがゴロゴロあるはずなのでここでは割愛。
Pythonそのもののバージョン及び天気図の描写に用いたパッケージと使用したバージョンは以下の通り。一応、それぞれの用途についてもコメントアウトの形で記してしておく。
Python 3.8.18
Package Version
------------------- ------------
numpy 1.24.4 # 数値計算
pandas 2.0.3 # CSVファイルの読み込み
scipy 1.10.1 # AMeDASや格子点の座標データ間を補間
matplotlib 3.7.2 # 図、グラフの描写
Cartopy 0.21.1 # 日本地図の描写CSVファイルの確認
AMeDAS観測点の観測データについては以下の気象庁のサイトからダウンロードできる。今回は2018年7月15日の、気圧を観測している観測点に限定してダウンロード。突風の発生時間は15~16時だが追々時間変化も確認できるよう時間は09時から21時までとしておくか。
ダウンロードしたデータには様々な文字情報も付随している。できるだけPythonで読み込みやすくできるようにしてやりたい。
まず、CSVファイルを扱ったことのある人ならわかるだろうが、とにかく文字化け防止や読み込み時の形式指定が面倒くさすぎるのでデータは英数字だけにしてしまう。
また、本来ダウンロードしたCSVファイルにはデータを扱ううえで欠かせない情報の精度に関するデータ(品質情報)もある。しかしこのデータまで逐一読み込ませて参照していたのでは本来の予定を1年超過して、なおかつもう1年使うことになったノートPCくんの負担なので、泣く泣く消すことにする。
あと何より緯度経度のデータが無い。これなくして天気図など書きようがないので10進法表記で打ち込んでやる。そうしてできたCSVファイルはこんな感じ。
所在地についても記録しているのは後々正確な位置にプロットできているかの確認のためだ。
日本地図の描写
天気図を描写するうえで日本地図は絶対に必要だ。まずはその描写からしていきたいが……ぶっちゃけてしまうと偉大なる先駆者の素晴らしいサイトに記載されていたコードをほとんどそのまま利用してしまっている。なので描写の仕方については当該サイトを参考されたし。
一応、描写範囲をあらかじめ指定しておくことで逐一書き換える必要の無いようにはした。少し先の話だが、格子点の作成にも使える。
# 描写範囲
W = 136 # 西端
E = 140 # 東端
S = 34 # 南端
N = 38 # 北端
-------------
# 一例としてこんな感じ
ax.set_extent([W, E, S, N]) # 図の範囲を指定するコード
# リスト[値, 値, …]にしておかないとエラー吐くぞ。何も知らなかった俺はいっぺんやらかした。できた地図はこれだ。白地図では味気ないのでせめてもの海に色を付けてみた。
地図に観測点を表示する
まずはpandasでCSVファイルを読み込んでもらう。
import pandas as pd # pandasで読み込んでもらうためのおまじない
# CSV読み込み
df = pd.read_csv('gust_20180715.csv') # gust = 突風このCSVファイルから緯度と経度のデータを抜き出したい。とりあえずこんな感じでやってみた。
# CSVから経度、緯度のデータを指定して抽出
lon = df['lon'] # 経度(x方向)。lonの列から抜き出してね、という意になるらしい。
lat = df['lat'] # 緯度(y方向)。latの(以下同文
print(lat)
print(lon) # 緯度から表示したほうが見やすいのでこの順にする。
---------------------------------------
0 34.628333
1 35.166667
2 35.400000
3 36.155000
4 34.940000
5 34.733333
6 34.068333
7 34.761667
8 35.275000
9 36.055000
10 36.588333
11 37.390000
12 36.708333
13 36.791667
14 37.106667
15 37.893333
16 36.661667
17 36.246667
18 36.045000
19 35.523333
20 34.753333
21 34.603333
22 34.975000
23 35.113333
24 34.603333
25 35.045000
26 36.405000
27 35.666667
28 35.438333
29 36.548333
30 35.990000
31 36.150000
32 35.691667
33 34.748333
34 34.123333
35 34.986667
36 37.488333
Name: lat, dtype: float64
0 137.093333
1 136.965000
2 136.761667
3 137.253333
4 136.580000
5 136.518333
6 136.193333
7 136.141667
8 136.243333
9 136.221667
10 136.633333
11 136.895000
12 137.201667
13 137.055000
14 138.246667
15 139.018333
16 138.191667
17 137.970000
18 138.108333
19 137.821667
20 137.711667
21 138.213333
22 138.403333
23 138.925000
24 138.841667
25 139.091667
26 139.060000
27 138.553333
28 139.651667
29 139.868333
30 139.073333
31 139.380000
32 139.750000
33 139.361667
34 139.520000
35 139.865000
36 139.910000
Name: lon, dtype: float64とりあえずデータは正しく抜き出せているようだ。
ではここから日本地図にマーカーをプロットしてもらおう。
# AMeDAS観測点プロット
ax.scatter(lat, lon, marker="s", color="magenta", s=25, edgecolors="none", transform=ccrs.PlateCarree(), zorder=6)
# markerで形状指定。sはsquare、つまり四角。
# マーカーの形状によってはedgecolorsを指定するとエラー吐くことも。"+"はダメだった。
# transform=ccrs.PlateCarree()は正距円筒図法で描いてねというおまじない的なやつ。これ抜くと座標がずれたり表示されないことも。
# zorderは図の描写の優先順位。マーカープロットまでの段階であれば、とりあえずこいつの値を一番デカくしてやればよい。
plt.show()
よし、ここまで順調に来られたぞ。この調子でマーカーに地名も重ねて表示してやろう。……と意気込んでいたうえはーすであったが、実は初回作成時にここで大きく躓いている。次の話はこの辺りの話からしていこう。
【次の話】