動機
YouTube で見かけた以下の動画を作ってみたくなった．ここまで綺麗なのは作れませんが．．
完成したもの
コードの概要
libcomcat というpythonのライブラリを通じてアメリカUSGSの地震データベース(ANSS Comprehensive Earthquake Catalog, ComCat)にアクセスできるようなのでそれを使ってデータを取得します．
取得したデータには日時，場所，マグニチュードなどの情報があるので，それを使って世界地図上に震源地をプロット，アニメーションを作成していきます．
必要なライブラリ
- libcomcat-python: 地震データへのアクセス
- cartopy: 世界地図データと描画用
- pandas: データハンドリング
- matplotlib: 画像描画用
- numpy: 主に配列操作
- tqdm: 進捗の表示用，無くてもOK
- opencv: 画像のリストから動画にする用
工夫した点
マグニチュードに応じて円の大きさと円が消えていくスピードを変化させました．円の濃さはマグニチュードと発生日からの経過時間に応じて薄くなるようにし，減衰していく様子を表しました．
コード要点解説
コード23-28行目
ここで地震データを取得しています．今回は比較的大きめの地震に絞るため，マグニチュード5以上のデータだけを取得するようにしています．
### retrieve the earthquake data
time_events = search(
starttime=datetime.datetime(2011, 1, 1, 0, 0),
endtime=datetime.datetime(2011, 12, 31, 23, 59),
minmagnitude = 5.0
)
コード30-43行目
ここで基本的なデータの作成を行なっています．個人的にPandasのデータフレームが扱いやすいのでそれに変換．
### create dataframe of earthqueakes
time_list = [e.time for e in time_events]
mag_list = [e.magnitude for e in time_events]
lat_list = [e.latitude for e in time_events]
long_list = [e.longitude for e in time_events]
dep_list = [e.depth for e in time_events]
df = pd.DataFrame({
'Date time': time_list,
'Magnitude': mag_list,
'Latitude' : lat_list,
'Longitude': long_list,
'Depth' : dep_list
})
あとはひたすらグラフを描いていくだけです．
発生日時からの経過日数をグラフの円の透過率に使いたいので，プロットする日付のリストを用意（コード46行目），それを使ってループを回します．
### count days of data collection
days = (df.iloc[-1,:]['Date time'] - df.iloc[0,:]['Date time']).days
なお，世界地図はこの1行で下準備完了です（72行目）．
### prepare axis
ax1 = plt.axes(projection=ccrs.PlateCarree(central_longitude=150))
ループ内でのポイントは，まず，発生日からの経過日数を評価するところです（83行目）．これを使って円の濃度（透過率）をコントロールして地震が減衰してく様子を表します．
### difference between current plot and past date
days_passed = np.array([(date - d.date()).days for d in plot_df['Date time']])
次に上で用意した経過日数とマグニチュードから円の透過率を0-1で求めます．減衰と言えば指数関数なので経過日数に応じて指数関数的に薄くなるようにしています．マグニチュードで割ることで，大きな地震の時は減衰が弱くなるようにして，あとに残る影響が長い様子を表しています．
### determine transparency depending on days passed and magnitude
alpha = np.exp(- alpha_const * days_passed / plot_df['Magnitude'])
個人的に情報収集に時間がかかったのが，太平洋を分割したくなかったので日本が真ん中あたりにくるマップにすることです．それ自体は72行目のcentral_longitudeの値を調整するだけでいいのですが，肝心の震源地が調整前の場所から動かずにいました．これは描写時のオプションtransferにcartopyのオブジェクトを渡すことで解決するようです（91-98行目）．
### plot eathquake data on the date
ax1.scatter(
plot_df['Longitude'],
plot_df['Latitude'],
s=np.exp(size_const * plot_df['Magnitude']),
facecolor='None', edgecolors=lc,
transform=ccrs.PlateCarree()
)
ちなみに付属しているmakemovie.pyは連番の画像から動画を作る関数です．よく使います．
import cv2
import glob
import tqdm
def make_movie(image_file_path, frame_rate=30.0, video_name='video', image_format='png'):
filelist = glob.glob(image_file_path + '/*.' + image_format)
filelist.sort()
im0 = cv2.imread(filelist[0])
height, width, color = im0.shape
fourcc = cv2.VideoWriter_fourcc('m', 'p', '4', 'v')
video = cv2.VideoWriter('{0}.mp4'.format(video_name), fourcc, frame_rate, (width, height))
for file in tqdm.tqdm(filelist):
img = cv2.imread(file)
video.write(img)
if __name__ == '__main__':
make_movie('data/test_api_002', video_name='EQ_movie')
最後に
冒頭の動画と本質的には似たようなものが作れました．見た目とか音とかまで含めるとまだまだですが，個人的には満足したのでこの話はここで終わろうと思います．
参考文献
座標の変換にtransformが有効であることを発見
usgsの地震データベース
USGS地震データベースのapi