0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

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

Last updated at Posted at 2021-12-31

はじめに

お仕事で、軽石の漂流シミュレーションをしている。結果の数値データはYouTubeの記述欄から入手可能にしている。

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

データをPandasで読む

粒子計算はJAXAが人工衛星で軽石を発見したところから計算を始めているので、軽石が発見された日付からスタートした粒子毎ファイルになっている。それを個別にPandasread_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を作ってしまえば、いかようにも解析できる。ある緯度経度範囲にある粒子数の数を数えてもいいだろう(例:海域毎に全体の内どれだけ分布しているかの割合(%)とその時間変化を見た図)。

ここではmatplotlibcartopyで沖縄周辺のアニメーションを作成してみた。アニメーションにするには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は
animation.gif

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

0
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?