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

# はじめに

お仕事で、[軽石の漂流シミュレーション](https://www.youtube.com/playlist?list=PLKT1Tlr-tdGG85epmm6U9Gg7woBV4gd6p)をしている。結果の数値データはYouTubeの記述欄から入手可能にしている。

元々のデータは我々の業界でよく使うnetcdfで取り扱っているのだが、公開しているデータは誰でも読めるということでcsv形式にしている。とは言え、取り扱いやすいかは別問題なので、自分でcsvデータを読むプログラムをPythonで書いてみた。

# データをPandasで読む

粒子計算はJAXAが人工衛星で軽石を発見したところから計算を始めているので、軽石が発見された日付からスタートした粒子毎ファイルになっている。それを個別に[Pandas](https://pandas.pydata.org/)の[read_csv](https://pandas.pydata.org/docs/reference/api/pandas.read_csv.html)で読んで、一つのDataFrameにまとめることにする。

プログラムは以下の通り。

- (1) start,end で必要な時刻の始まりと終わりと指定する。
- (2) 海流のみで計算(windage=0)、海流+風速x0.5%で計算(windage=0.5)、海流+風速x1.0%で計算(windage=1.0)の3つのシナリオの計算しており、ここではwindage=0.5のケースを読む。
- (3) ファイルはディレクトリcsvの下にあるとする。[glob](https://docs.python.org/ja/3/library/glob.html)でwindageが共通するファイルを全て列挙する。forで回して、個別にファイルを読んで、合成する。
- (4) read_csvでcsvファイルを読む。行が粒子番号、列が日付になっているので、転置し、行が日付、列が粒子番号になるようにする。
- (5) 列の粒子番号はNで通しでカウントして、全体の通し番号にする。
- (6) 行はdatatime形式の日付時間にする。
- (7) [join](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.join.html)を使って日付(index)を軸に、一つのDataFrameに合成する。



```python
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の中身は以下の通り。



```python
lonall
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>P00001</th>
      <th>P00002</th>
      <th>P00003</th>
      <th>P00004</th>
      <th>P00005</th>
      <th>P00006</th>
      <th>P00007</th>
      <th>P00008</th>
      <th>P00009</th>
      <th>P00010</th>
      <th>...</th>
      <th>P23991</th>
      <th>P23992</th>
      <th>P23993</th>
      <th>P23994</th>
      <th>P23995</th>
      <th>P23996</th>
      <th>P23997</th>
      <th>P23998</th>
      <th>P23999</th>
      <th>P24000</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>2021-12-31 09:00:00</th>
      <td>123.896560</td>
      <td>127.66187</td>
      <td>123.704790</td>
      <td>126.180305</td>
      <td>124.029335</td>
      <td>123.290360</td>
      <td>125.401794</td>
      <td>128.02122</td>
      <td>125.44872</td>
      <td>123.742170</td>
      <td>...</td>
      <td>128.00153</td>
      <td>128.00153</td>
      <td>128.00452</td>
      <td>128.00186</td>
      <td>128.0107</td>
      <td>128.01974</td>
      <td>128.00093</td>
      <td>122.738075</td>
      <td>122.726560</td>
      <td>122.738400</td>
    </tr>
    <tr>
      <th>2021-12-31 10:00:00</th>
      <td>123.895645</td>
      <td>127.66187</td>
      <td>123.708250</td>
      <td>126.184425</td>
      <td>124.022410</td>
      <td>123.282050</td>
      <td>125.401794</td>
      <td>128.02122</td>
      <td>125.44872</td>
      <td>123.744540</td>
      <td>...</td>
      <td>128.00153</td>
      <td>128.00153</td>
      <td>128.00452</td>
      <td>128.00186</td>
      <td>128.0107</td>
      <td>128.01974</td>
      <td>128.00093</td>
      <td>122.737110</td>
      <td>122.725494</td>
      <td>122.737434</td>
    </tr>
    <tr>
      <th>2021-12-31 11:00:00</th>
      <td>123.894730</td>
      <td>127.66187</td>
      <td>123.712660</td>
      <td>126.189220</td>
      <td>124.017136</td>
      <td>123.276200</td>
      <td>125.401794</td>
      <td>128.02122</td>
      <td>125.44872</td>
      <td>123.748350</td>
      <td>...</td>
      <td>128.00153</td>
      <td>128.00153</td>
      <td>128.00452</td>
      <td>128.00186</td>
      <td>128.0107</td>
      <td>128.01974</td>
      <td>128.00093</td>
      <td>122.739790</td>
      <td>122.728690</td>
      <td>122.740160</td>
    </tr>
    <tr>
      <th>2021-12-31 12:00:00</th>
      <td>123.893814</td>
      <td>127.66187</td>
      <td>123.717240</td>
      <td>126.195435</td>
      <td>124.013220</td>
      <td>123.272250</td>
      <td>125.401794</td>
      <td>128.02122</td>
      <td>125.44872</td>
      <td>123.753105</td>
      <td>...</td>
      <td>128.00153</td>
      <td>128.00153</td>
      <td>128.00452</td>
      <td>128.00186</td>
      <td>128.0107</td>
      <td>128.01974</td>
      <td>128.00093</td>
      <td>122.745220</td>
      <td>122.735290</td>
      <td>122.745605</td>
    </tr>
    <tr>
      <th>2021-12-31 13:00:00</th>
      <td>123.892550</td>
      <td>127.66187</td>
      <td>123.721820</td>
      <td>126.203260</td>
      <td>124.010086</td>
      <td>123.269104</td>
      <td>125.401794</td>
      <td>128.02122</td>
      <td>125.44872</td>
      <td>123.758130</td>
      <td>...</td>
      <td>128.00153</td>
      <td>128.00153</td>
      <td>128.00452</td>
      <td>128.00186</td>
      <td>128.0107</td>
      <td>128.01974</td>
      <td>128.00093</td>
      <td>122.752490</td>
      <td>122.743660</td>
      <td>122.752880</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th>2022-01-13 05:00:00</th>
      <td>123.563030</td>
      <td>127.66187</td>
      <td>123.022540</td>
      <td>133.489500</td>
      <td>123.614210</td>
      <td>122.686890</td>
      <td>125.401794</td>
      <td>128.02122</td>
      <td>125.44872</td>
      <td>122.957820</td>
      <td>...</td>
      <td>128.00153</td>
      <td>128.00153</td>
      <td>128.00452</td>
      <td>128.00186</td>
      <td>128.0107</td>
      <td>128.01974</td>
      <td>128.00093</td>
      <td>124.752440</td>
      <td>124.500084</td>
      <td>124.724220</td>
    </tr>
    <tr>
      <th>2022-01-13 06:00:00</th>
      <td>123.565506</td>
      <td>127.66187</td>
      <td>123.018654</td>
      <td>133.518370</td>
      <td>123.621200</td>
      <td>122.695340</td>
      <td>125.401794</td>
      <td>128.02122</td>
      <td>125.44872</td>
      <td>122.952740</td>
      <td>...</td>
      <td>128.00153</td>
      <td>128.00153</td>
      <td>128.00452</td>
      <td>128.00186</td>
      <td>128.0107</td>
      <td>128.01974</td>
      <td>128.00093</td>
      <td>124.763460</td>
      <td>124.506220</td>
      <td>124.734490</td>
    </tr>
    <tr>
      <th>2022-01-13 07:00:00</th>
      <td>123.568214</td>
      <td>127.66187</td>
      <td>123.015495</td>
      <td>133.546750</td>
      <td>123.627556</td>
      <td>122.703840</td>
      <td>125.401794</td>
      <td>128.02122</td>
      <td>125.44872</td>
      <td>122.949585</td>
      <td>...</td>
      <td>128.00153</td>
      <td>128.00153</td>
      <td>128.00452</td>
      <td>128.00186</td>
      <td>128.0107</td>
      <td>128.01974</td>
      <td>128.00093</td>
      <td>124.773750</td>
      <td>124.512910</td>
      <td>124.744316</td>
    </tr>
    <tr>
      <th>2022-01-13 08:00:00</th>
      <td>123.571590</td>
      <td>127.66187</td>
      <td>123.013520</td>
      <td>133.575130</td>
      <td>123.633510</td>
      <td>122.711900</td>
      <td>125.401794</td>
      <td>128.02122</td>
      <td>125.44872</td>
      <td>122.948440</td>
      <td>...</td>
      <td>128.00153</td>
      <td>128.00153</td>
      <td>128.00452</td>
      <td>128.00186</td>
      <td>128.0107</td>
      <td>128.01974</td>
      <td>128.00093</td>
      <td>124.783646</td>
      <td>124.520460</td>
      <td>124.753930</td>
    </tr>
    <tr>
      <th>2022-01-13 09:00:00</th>
      <td>123.575870</td>
      <td>127.66187</td>
      <td>123.012794</td>
      <td>133.603520</td>
      <td>123.639460</td>
      <td>122.719170</td>
      <td>125.401794</td>
      <td>128.02122</td>
      <td>125.44872</td>
      <td>122.948654</td>
      <td>...</td>
      <td>128.00153</td>
      <td>128.00153</td>
      <td>128.00452</td>
      <td>128.00186</td>
      <td>128.0107</td>
      <td>128.01974</td>
      <td>128.00093</td>
      <td>124.793680</td>
      <td>124.528910</td>
      <td>124.763890</td>
    </tr>
  </tbody>
</table>
<p>313 rows × 24000 columns</p>
</div>



# アニメーションを作成

DataFrameを作ってしまえば、いかようにも解析できる。ある緯度経度範囲にある粒子数の数を数えてもいいだろう(例:[海域毎に全体の内どれだけ分布しているかの割合(%)とその時間変化を見た図](https://drive.google.com/file/d/19criUOXRsr7lTy-0sBCblOEZ1ZWoZfWh/view))。

ここでは[matplotlib](https://matplotlib.org/)と[cartopy](https://scitools.org.uk/cartopy/docs/latest/)で沖縄周辺のアニメーションを作成してみた。アニメーションにするにはmatplotlibの[FuncAnimation](https://matplotlib.org/stable/api/_as_gen/matplotlib.animation.FuncAnimation.html)を使っている。

```python
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は
![animation.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/77672/5212bd12-5e54-af27-31bb-9bd61f3ccb8b.gif)

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