はじめに
お仕事で、軽石の漂流シミュレーションをしている。結果の数値データはYouTubeの記述欄から入手可能にしている。
元々のデータは我々の業界でよく使うnetcdfで取り扱っているのだが、公開しているデータは誰でも読めるということでcsv形式にしている。とは言え、取り扱いやすいかは別問題なので、自分でcsvデータを読むプログラムをPythonで書いてみた。
データをPandasで読む
粒子計算はJAXAが人工衛星で軽石を発見したところから計算を始めているので、軽石が発見された日付からスタートした粒子毎ファイルになっている。それを個別にPandasのread_csvで読んで、一つのDataFrameにまとめることにする。
プログラムは以下の通り。
- (1) start,end で必要な時刻の始まりと終わりと指定する。
- (2) 海流のみで計算(windage=0)、海流+風速x0.5%で計算(windage=0.5)、海流+風速x1.0%で計算(windage=1.0)の3つのシナリオの計算しており、ここではwindage=0.5のケースを読む。
- (3) ファイルはディレクトリcsvの下にあるとする。globでwindageが共通するファイルを全て列挙する。forで回して、個別にファイルを読んで、合成する。
- (4) read_csvでcsvファイルを読む。行が粒子番号、列が日付になっているので、転置し、行が日付、列が粒子番号になるようにする。
- (5) 列の粒子番号はNで通しでカウントして、全体の通し番号にする。
- (6) 行はdatatime形式の日付時間にする。
- (7) joinを使って日付(index)を軸に、一つのDataFrameに合成する。
import pandas as pd
import glob
start="2021-12-31 09:00:00" # (1)
end="2022-01-13 09:00:00" # (1)
windage=0.5 # (2)
files=sorted(glob.glob(f"csv/latitude_from*_windage{windage:.1f}.csv")) # (3) 同じwindageのファイルを列挙する。
N=0 # (5)
for f in files: # (3)
lat=pd.read_csv(f,index_col=0).T # (4) 緯度データを読み込み、列と行を入れ替える
lon=pd.read_csv(f.replace("latitude","longitude"),index_col=0).T # (4) 同じく経度
lat.columns=[f"P{nn:05d}" for nn in range(N+1,N+len(lat.columns)+1)] # (5) 粒子に通しで番号を振る (列)
lon.columns=lat.columns # (5) 同じく緯度
lat.index=pd.to_datetime(lat.index) # (6) indexをdatatime formatにする
lon.index=pd.to_datetime(lon.index) # (6) 同じく緯度
if N==0:
latall=lat.loc[start:end]
lonall=lon.loc[start:end]
else:
latall=latall.join(lat.loc[start:end]) # (7) start から endまで選んで、indexを軸に、dataframeを合成する
lonall=lonall.join(lon.loc[start:end]) # (7)
N += len(lat.columns) # 粒子通し番号
lonallは緯度の、latallは経度のDataFrameになる。例えばlatallの中身は以下の通り。
lonall
P00001 | P00002 | P00003 | P00004 | P00005 | P00006 | P00007 | P00008 | P00009 | P00010 | ... | P23991 | P23992 | P23993 | P23994 | P23995 | P23996 | P23997 | P23998 | P23999 | P24000 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
2021-12-31 09:00:00 | 123.896560 | 127.66187 | 123.704790 | 126.180305 | 124.029335 | 123.290360 | 125.401794 | 128.02122 | 125.44872 | 123.742170 | ... | 128.00153 | 128.00153 | 128.00452 | 128.00186 | 128.0107 | 128.01974 | 128.00093 | 122.738075 | 122.726560 | 122.738400 |
2021-12-31 10:00:00 | 123.895645 | 127.66187 | 123.708250 | 126.184425 | 124.022410 | 123.282050 | 125.401794 | 128.02122 | 125.44872 | 123.744540 | ... | 128.00153 | 128.00153 | 128.00452 | 128.00186 | 128.0107 | 128.01974 | 128.00093 | 122.737110 | 122.725494 | 122.737434 |
2021-12-31 11:00:00 | 123.894730 | 127.66187 | 123.712660 | 126.189220 | 124.017136 | 123.276200 | 125.401794 | 128.02122 | 125.44872 | 123.748350 | ... | 128.00153 | 128.00153 | 128.00452 | 128.00186 | 128.0107 | 128.01974 | 128.00093 | 122.739790 | 122.728690 | 122.740160 |
2021-12-31 12:00:00 | 123.893814 | 127.66187 | 123.717240 | 126.195435 | 124.013220 | 123.272250 | 125.401794 | 128.02122 | 125.44872 | 123.753105 | ... | 128.00153 | 128.00153 | 128.00452 | 128.00186 | 128.0107 | 128.01974 | 128.00093 | 122.745220 | 122.735290 | 122.745605 |
2021-12-31 13:00:00 | 123.892550 | 127.66187 | 123.721820 | 126.203260 | 124.010086 | 123.269104 | 125.401794 | 128.02122 | 125.44872 | 123.758130 | ... | 128.00153 | 128.00153 | 128.00452 | 128.00186 | 128.0107 | 128.01974 | 128.00093 | 122.752490 | 122.743660 | 122.752880 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
2022-01-13 05:00:00 | 123.563030 | 127.66187 | 123.022540 | 133.489500 | 123.614210 | 122.686890 | 125.401794 | 128.02122 | 125.44872 | 122.957820 | ... | 128.00153 | 128.00153 | 128.00452 | 128.00186 | 128.0107 | 128.01974 | 128.00093 | 124.752440 | 124.500084 | 124.724220 |
2022-01-13 06:00:00 | 123.565506 | 127.66187 | 123.018654 | 133.518370 | 123.621200 | 122.695340 | 125.401794 | 128.02122 | 125.44872 | 122.952740 | ... | 128.00153 | 128.00153 | 128.00452 | 128.00186 | 128.0107 | 128.01974 | 128.00093 | 124.763460 | 124.506220 | 124.734490 |
2022-01-13 07:00:00 | 123.568214 | 127.66187 | 123.015495 | 133.546750 | 123.627556 | 122.703840 | 125.401794 | 128.02122 | 125.44872 | 122.949585 | ... | 128.00153 | 128.00153 | 128.00452 | 128.00186 | 128.0107 | 128.01974 | 128.00093 | 124.773750 | 124.512910 | 124.744316 |
2022-01-13 08:00:00 | 123.571590 | 127.66187 | 123.013520 | 133.575130 | 123.633510 | 122.711900 | 125.401794 | 128.02122 | 125.44872 | 122.948440 | ... | 128.00153 | 128.00153 | 128.00452 | 128.00186 | 128.0107 | 128.01974 | 128.00093 | 124.783646 | 124.520460 | 124.753930 |
2022-01-13 09:00:00 | 123.575870 | 127.66187 | 123.012794 | 133.603520 | 123.639460 | 122.719170 | 125.401794 | 128.02122 | 125.44872 | 122.948654 | ... | 128.00153 | 128.00153 | 128.00452 | 128.00186 | 128.0107 | 128.01974 | 128.00093 | 124.793680 | 124.528910 | 124.763890 |
313 rows × 24000 columns
アニメーションを作成
DataFrameを作ってしまえば、いかようにも解析できる。ある緯度経度範囲にある粒子数の数を数えてもいいだろう(例:海域毎に全体の内どれだけ分布しているかの割合(%)とその時間変化を見た図)。
ここではmatplotlibとcartopyで沖縄周辺のアニメーションを作成してみた。アニメーションにするにはmatplotlibのFuncAnimationを使っている。
import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib.animation as animation
import cartopy.crs as ccrs
import cartopy.feature as cfeature
DATES=pd.date_range(start=start,end=end,freq="H")
def update(i,ax):
if i != 0:
plt.cla() #現在描写されているグラフを消去
plon=lonall.iloc[i]
plat=latall.iloc[i]
cl=ax.plot(plon,plat,"b.")
ax.coastlines()
ax.add_feature(cfeature.LAND,color="gray")
ax.set_title(DATES[i].strftime("%Y/%m/%d %H"),fontsize=18)
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=2, color='gray', alpha=0.5, linestyle='--')
gl.top_lables = False
gl.right_labels = False
ax.set_extent([127,129,26,27.5])
fig=plt.figure(figsize=(12,8))
proj=ccrs.PlateCarree()
ax = fig.add_subplot(111, projection=proj)
ani = animation.FuncAnimation(fig, update, fargs = (ax,),interval = 100, frames = len(latall))
writergif = animation.PillowWriter(fps=30)
ani.save("animation.gif", writer = writergif)
以上のjupyter notebookは
https://gist.github.com/tmiyama/f6d02bcc3d894a915a2346fbc1ae06af