軽石漂流シミュレーションのcsvファイルを読む

はじめに

お仕事で、軽石の漂流シミュレーションをしている。結果の数値データは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)

結果のanimation.gifは

以上のjupyter notebookは
https://gist.github.com/tmiyama/f6d02bcc3d894a915a2346fbc1ae06af