この記事は、学生LT Advent Calendar 2019 - Qiita11日目の記事です。
こんにちは。
私は大気科学を専攻している大学院生です。
今回は先行記事でもいくつか見られますが、Pythonを用いた気象データの描画に焦点を当てて投稿しました。
紹介する描画コードはGitHubにて公開しておりますので、気になっていただけた方はぜひダウンロード下さい。
はじめに
気象データの描画や取り扱いに関連して、以下のような記事が見受けられました。
非常に分かりやすくまとまっておりますが、どちらも一つの要素の描画までの紹介に留まっています。
私は下記のようなある程度の完成度を持つ図のデータの整理から描画までのコードを紹介することで上記の記事との差別化を図りたいと思います。
今回、作成する完成予定図はこちら。
ターゲットは2016年台風10号の事例を捉えています。
| ・図１ | ・図２ |
|---|---|
 |
 |
各々の図の紹介
図の自分なりの解釈を紹介します。(...需要は低いかもしれませんが)
- 図１
2016年８月21日06時から2016年８月31日00時までの気象要素の組み合わせ描画をしております。
コンターが高度場を、青いベクトルが10m/s以上の風速場を、斜線+黄色部分は指定した高度場における特に強い低気圧部分を表しています。
※強い低気圧部分の抽出方法は今回は自分の裁量で判断しましたが、あくまで学術的ではなく描画方法の紹介なので深く凝ってはおりません。
→上陸するまでに徐々に勢力が拡大している様子が確認できます。
- 図2
ここ20年間(2000-2019年)の台風がどの地点をより通ってきていたのかをヘキサグリッドで表現しました。またプロットでは特定の台風経路を表しており、今回は２０１６年台風10号の経路を図１にならいプロットしてしています。
→これまでの大方の台風が台湾の東側付近を通過して北進し、日本に接近してきたことがなんとなくヘキサグリッドから読み取れますが、2016年の台風10号は特殊なルートであったことが確認できます。
用意するデータ
- 全球数値予報モデルGPV
京都大学生存圏研究所（RISH: Research Institute for Sustainable Humanosphere）から取得することができます。詳しい格子情報や配信時間等は気象行支援センターをご覧ください。
※上記のサイトご利用時は、下記の点にご注意お願いいたします。
教育研究機関向けにデータを提供しています。企業活動等のためにデータを頻繁に必要とされる方は、気象業務支援センターからデータを直接購入し、データ提供スキーム全体の維持発展にご協力ください。（サイトより引用）
今回使用するデータは上記のサイトで入手したデータをwgirb2を用いて東西風・南北風・鉛直風・気温・高度場のデータをbinary形式で切り出しています。
切り出したデータフォーマットをFortran90で記述しました。...参考程度に。
format.f90
integer, parameter :: nx = 720, ny = 361, np = 12
integer, parameter :: elem = 5
real(4) :: gpv(elem, hgt, nx, ny)
real(4) :: pr(hgt)
character(4) :: elem_list(elem)
data pr / 1000., 925., 850., 700., 600., 500., 400., 300., 250., 200., 150., 100. /
data elem_list / "U ", "V ", "W ", "T ", "Z " /
open(*, fille=*, form='unformatted', access='direct', recl=4*nx*ny*np*elem)
close(*)
- 台風の位置表(csv形式)
気象庁から年度別で取得することができます。
台風発生時から温帯低気圧になるまでの情報が6時間間隔で記載されています。
実行コード
どちらのコードもmain_driverを基準に読んでいただければと思います。
実行をされる方はPATHの設定を変更お願いいたします。
気象要素描画コード(図１)
はじめに気象要素を描画したコードを紹介します。
weather_bin2plot.py
# -*- coding: utf-8 -*-
"""
Created on 2019.12.11
@author: Toyo_Daichi
"""
import numpy as np
import pandas as pd
import glob
import itertools
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap
class mapping:
def __init__(self):
pass
def gpv_data_coef(self):
nx, ny, hgt = 720, 361, 12
elem = 5
return nx, ny, hgt, elem
def prj_coef(self, nx, ny):
dx, dy = 0.5, 0.5
lon, lat = [], []
for ix in range(nx):
lon += [ float('{:.2f}'.format(dx*ix)) ]
for iy in range(ny):
lat += [ float('{:.2f}'.format(-90. + dy*iy)) ]
X, Y = np.meshgrid(lon, lat)
return X, Y
def preparating_data(self, yyyy, mm, dd, hh):
time_list = []
num_time_list = len(time_list)
month_thirtyone = [ 1, 3, 5, 7, 8, 10, 12 ]
month_thirty = [ 4, 6, 9, 11 ]
month_twntynine = [ 2 ]
while num_time_list < 40:
time_list.append(str(yyyy) + str('%02d' % mm) + str('%02d' % dd) + str('%02d' % hh) + '00')
hh = hh - 6
if hh < 0 and dd == 1:
mm, hh = mm - 1, 18
if mm in month_thirty:
dd = 30
elif mm in month_thirtyone:
dd = 31
elif mm in month_twntynine:
if yyyy % 4 == 0:
dd = 28
else:
dd =29
elif hh < 0:
dd, hh = dd - 1, 18
num_time_list += 1
return time_list
def setup_gpv_filelist(self, path, time_list):
filelist = []
for i_file in time_list:
filelist.append(path + '/data')
return filelist
def open_gpv_filelist(self, gpv_filelist, nx, ny, hgt, elem):
data = [ []*i for i in range(len(gpv_filelist)) ]
for num_gpv, gpv_file in enumerate(gpv_filelist):
with open(gpv_file, 'rb') as ifile:
"""
elem: 1-u, 2-v, 3-w, 4-tmp, 5-height
"""
data[num_gpv] = np.fromfile(ifile, dtype='>f', sep = '').reshape(elem, hgt, ny, nx)
print('..... Preparating data for ' + gpv_file, ', shapes :: ', data[num_gpv].shape)
return data
def main_mapping_tool(self, mode, path, time_list, nx, ny, *, gpv_datalist='None', snap_step=0, level='1000'):
fig, ax = plt.subplots()
outpath = path + '/fig/'
lat_min, lat_max = 17, 50
lon_min, lon_max = 120, 155
mapping = Basemap( projection='lcc', resolution="i", lat_0=35, lon_0=140, fix_aspect=(1,1), llcrnrlat=lat_min, urcrnrlat=lat_max, llcrnrlon=lon_min, urcrnrlon=lon_max )
mapping.drawcoastlines(color='black', linewidth=0.5)
mapping.drawmeridians(np.arange(0, 360, 5), labels=[False, True, False, True], fontsize='small', color='gray', linewidth=0.5)
mapping.drawparallels(np.arange(-90, 90, 5), labels=[True, False, False, True], fontsize='small', color='gray', linewidth=0.5)
lon_list, lat_list = self.prj_coef(nx, ny)
x, y = mapping(lon_list, lat_list)
if mode == 1 or mode == 2:
if level == "1000":
hPa, min_list = 0, [0, 200]
elif level == "925":
hPa, min_list = 1, [0, 600]
elif level == "850":
hPa, min_list = 2, [0, 750]
elif level == "700":
hPa, min_list = 3, [0, 1500]
elif level == "600":
hPa, min_list = 4, "None"
elif level == "500":
hPa, min_list = 5, "None"
elif level == "400":
hPa, min_list = 6, "None"
elif level == "300":
hPa, min_list = 7, "None"
elif level == "250":
hPa, min_list = 8, "None"
elif level == "200":
hPa, min_list = 9, "None"
elif level == "150":
hPa, min_list = 10, "None"
elif level == "100":
hPa, min_list = 11, "None"
gpv_u_data = gpv_datalist[snap_step][0][hPa]
gpv_v_data = gpv_datalist[snap_step][1][hPa]
speed = np.sqrt(gpv_u_data*gpv_u_data + gpv_v_data*gpv_v_data)
speed_list = list(range(0, 50, 25))
gpv_height_data = gpv_datalist[snap_step][4][hPa]
num_list = list(range(0, 7500, 10))
contour = mapping.contour(x, y, gpv_height_data, linewidths=0.25, linestyles='-', levels=num_list, colors='m')
contour.clabel(fmt='%1.1f', fontsize=6.5)
if not min_list == "None":
lines = mapping.contourf(x, y, gpv_height_data, min_list, alpha=0.5, hatches=['///'], lw=1., zorder=0, edgecolor='r', colors="y")
for i_nx, i_ny in itertools.product(range(nx), range(ny)):
if speed[i_ny][i_nx] > 10 and lon_min <= lon_list[i_ny][i_nx] <= lon_max and lat_min <= lat_list[i_ny][i_nx] <= lat_max:
print('...... Halfway step, ', i_nx, i_ny, speed[i_ny][i_nx])
vector = mapping.quiver(x[i_ny][i_nx], y[i_ny][i_nx], gpv_u_data[i_ny][i_nx], gpv_v_data[i_ny][i_nx], color='c', units='dots', scale=2.0, alpha=0.6)
plt.title(time_list[snap_step] + ' @GSM ' + level + 'hPa' , loc='left', fontsize=10)
plt.quiverkey(vector, 0.75, 0.9, 10, '10 m/s', labelpos='W', coordinates='figure')
plt.savefig(outpath + 'GPV_elem_' + time_list[snap_step] + '.png')
print('...... Saving fig :: ', outpath + 'GPV_elem_' + time_list[snap_step] + '.png')
#plt.show()
def main_driver(self, mode, indir, time_list, level):
nx, ny, hgt, elem = self.gpv_data_coef()
gpv_filelist = self.setup_gpv_filelist(indir, time_list)
gpv_datalist = self.open_gpv_filelist(gpv_filelist, nx, ny, hgt, elem)
if mode == 1:
snap_step = 0
self.main_mapping_tool(mode, indir, time_list, nx, ny, gpv_datalist=gpv_datalist, level=level, snap_step=snap_step)
elif mode == 2:
for snap_step in range(len(gpv_datalist)):
self.main_mapping_tool(mode, indir, time_list, nx, ny, gpv_datalist=gpv_datalist, level=level, snap_step=snap_step)
if __name__ == "__main__":
mapp = mapping()
#input dir
indir = '/your_directory/'
#target time
yyyy, mm, dd, hh = 2016, 8, 31, 0
time_list = mapp.preparating_data(yyyy, mm, dd, hh)
"""
Target altitude list
1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100 hPa
"""
level = '925'
#choose mode, if you append new func. more anynum.
"""
2019.12.11
mode 1: Normal weather element info at GPV GSM snap shot.
mode 2: Normal weather element info at GPV GSM gif.
"""
mode = 1
#main_driver
mapp.main_driver(mode, indir, time_list, level)
このコードで工夫した点は下記の点でした。
- 風速ベクトルを10m/s以上のみ描画する。
- 各高度の大きい値を斜線で描画する。
- gif形式の画像作成とスナップショットの画像作成をモードで選択できるようにする。
台風進路の描画コード(図２)
次に過去20年の台風進路度合いと特定の台風の進路を作成したコードを紹介します。
typhoon_csv2plot.py
# -*- coding: utf-8 -*-
"""
Created on 2019.12.11
@author: Toyo_Daichi
"""
import numpy as np
import pandas as pd
import glob
import itertools
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap
class mapping:
def __init__(self):
pass
def setup_csv_filelist(self, path, *, year='*'):
filelist = []
filelist = glob.glob(path + 'table' + str(year) + '.csv')
return filelist
def open_csv_filelist(self, csv_filelist, *, typhoon_number='None'):
csv_datalist = [ []*i for i in range(len(csv_filelist)) ]
list_num = list(range(11))
list_option = ( 'year', 'month', 'day', 'hour(UTC)', 'typhoon_number', 'typhoon_name', 'rank','latitude', 'longitude', 'central_pressure', 'max_wind')
for num_file, infile in enumerate(csv_filelist):
print('..... Preparating data for ' + str(num_file) + ' ' + str(infile))
tmp_data = pd.read_csv(infile, usecols=list_num, skiprows=1, names=list_option, sep=',')
if typhoon_number == "None":
tmp_lat_list = tmp_data['latitude'].values.tolist()
csv_datalist[num_file].append(tmp_lat_list)
tmp_lon_list = tmp_data['longitude'].values.tolist()
csv_datalist[num_file].append(tmp_lon_list)
tmp_centpre_list = tmp_data['central_pressure'].values.tolist()
csv_datalist[num_file].append(tmp_centpre_list)
else:
csv_datalist = []
specific_data = tmp_data.query('typhoon_number == %s' % typhoon_number)
tmp_lat_list = specific_data['latitude'].values.tolist()
csv_datalist.append(tmp_lat_list)
tmp_lon_list = specific_data['longitude'].values.tolist()
csv_datalist.append(tmp_lon_list)
tmp_centpre_list = specific_data['central_pressure'].values.tolist()
csv_datalist.append(tmp_centpre_list)
return csv_datalist
def main_mapping_tool(self, path, csv_datalist, csv_specific_datalist='None', typhoon_info='None'):
fig, ax = plt.subplots()
outpath = path + '/fig/'
lat_min, lat_max = 17, 50
lon_min, lon_max = 120, 155
mapping = Basemap( projection='lcc', resolution="i", lat_0=35, lon_0=140, fix_aspect=(1,1), llcrnrlat=lat_min, urcrnrlat=lat_max, llcrnrlon=lon_min, urcrnrlon=lon_max )
mapping.drawcoastlines(color='black', linewidth=0.5)
mapping.drawmeridians(np.arange(0, 360, 5), labels=[False, True, False, True], fontsize='small', color='gray', linewidth=0.5)
mapping.drawparallels(np.arange(-90, 90, 5), labels=[True, False, False, True], fontsize='small', color='gray', linewidth=0.5)
full_lat_list = list(map(lambda x: x[0], csv_datalist))
full_lon_list = list(map(lambda x: x[1], csv_datalist))
full_lat_list = np.sum(full_lat_list, axis=0)
full_lon_list = np.sum(full_lon_list, axis=0)
lat_list, lon_list = [], []
for i_num in range(len(full_lat_list)):
if(lat_min <= full_lat_list[i_num] <= lat_max and lon_min <= full_lon_list[i_num] <= lon_max):
lat_list.append(full_lat_list[i_num])
lon_list.append(full_lon_list[i_num])
x, y = mapping(lon_list, lat_list)
hexbin = mapping.hexbin(np.array(x), np.array(y), gridsize=[30, 30], cmap='Blues', edgecolors='gray', mincnt=8)
if not csv_specific_datalist == "None":
specific_lat_list, specific_lon_list = csv_specific_datalist[0], csv_specific_datalist[1]
specific_centpre_list = csv_specific_datalist[2]
lat_list, lon_list, centpre_list = [], [], []
for i_num in range(len(specific_lon_list)):
if(lat_min <= specific_lat_list[i_num] <= lat_max and lon_min <= specific_lon_list[i_num] <= lon_max):
lat_list.append(specific_lat_list[i_num])
lon_list.append(specific_lon_list[i_num])
centpre_list.append(specific_centpre_list[i_num])
case_x, case_y = mapping(lon_list, lat_list)
mapping.plot(case_x, case_y, linewidth=0.5, color='red', ls='--', marker='o', ms=2)
cbar = plt.colorbar(hexbin)
cbar.set_label('number', fontsize=8)
if not csv_specific_datalist == "None":
plt.title('Course of typhoon 2000-2019' + ', ' + 'Typhoon track: ' + str(typhoon_info[1]), loc='left', fontsize=9)
else:
plt.title('Course of typhoon 2000-2019', loc='left', fontsize=10)
plt.savefig(outpath + 'typhoon_strength2000-2019' + '_track_' + str(typhoon_info[1]) + '.png')
print('..... savefig ' + outpath + 'typhoon_strength2000-2019' + '_track_' + str(typhoon_info[1]) + '.png')
#plt.show()
def main_driver(self, indir, *, typhoon_info='None'):
csv_filelist = self.setup_csv_filelist(indir)
csv_datalist = self.open_csv_filelist(csv_filelist)
"""
csv_datalist.shape = [year][latitude][longitude][central_pressure]
"""
if not typhoon_info == "None":
print('..... Check specific case filelist')
csv_specific_filelist = self.setup_csv_filelist(indir, year=typhoon_info[0])
csv_specific_datalist = self.open_csv_filelist(csv_specific_filelist, typhoon_number=typhoon_info[1])
self.main_mapping_tool(indir, csv_datalist, csv_specific_datalist=csv_specific_datalist, typhoon_info=typhoon_info)
else:
self.main_mapping_tool(indir, csv_datalist)
if __name__ == "__main__":
mapp = mapping()
#input dir
indir = '/your_directory/'
"""
If you want to write a specific typhoon route, enter the typhoon information.
typhoon_info = [year, typhoon_number]
"""
typhoon_info = [2016, 1610]
#main_driver
mapp.main_driver(indir, typhoon_info=typhoon_info)
print('Normal END')
このコードで工夫した点は下記の点でした。
- ２０年間の台風経路の頻度をヘキサグリッドで描画する。
- 特定の台風番号を記述することで、その台風の経路を描画する。
今回使用した使用環境
申し訳程度の環境を記述しました。
- Python 3.6.7
requirement.txt
basemap==1.2.0
matplotlib==3.1.1
最後に
作成したコードの中には細々としていますが、いくつかの自分なりの工夫をしているので今後紹介していきたいです。
より良い書き方などお気づきになった方はコメント等でご教授いただければ幸いです。
最後まで読んでくださり、ありがとうございました。
参考にしたサイト
- 気象×Python 〜衛星降水データ〜 https://qiita.com/OSAKO/items/d54da432fff34d4be1d8