地震のデータを地図上にプロットして可視化する

はじめに

位置情報データの可視化にめっちゃ便利なPythonのライブラリfoliumを使って、このように地震の大きさを可視化してみた。
ついでにPILを使ったgifファイルの作り方についても簡単に紹介したいと思います。

データを集める

気象庁の震源データをスクレイピングしました。
気象庁では震源地を60進数で記載している為、データとして取得する際には10進数による座標変換を行います。

例
「39度 50分 9秒」を10進数に変換する。
「度」はそのままで、「分」を60、「秒」を3,600で割りそれらを足す。
　= 139 + ( 50 ÷ 60 ) + ( 9 ÷ 60 ÷ 60 )
　= 139.691749

以下取得したデータ例

日	時	緯度	経度	深さ	震度	地名
2019/2/1	0:09:29	39.8358333333333	142.351111111111	22	1.2	岩手県沖
2019/2/1	0:26:44	24.5	122.250833333333	27	4.2	台湾付近
2019/2/1	0:27:24	37.1177777777777	141.134444444444	49	2.4	福島県沖
2019/2/1	0:28:09	34.8833333333333	135.65	12	0.8	大阪府北部
2019/2/1	0:30:52	35.2847222222222	135.302222222222	9	0.6	京都府北部

地図上にプロットする

取得したデータをfoliumを使ってプロットしていきます。
foliumの使い方についてはこちらから確認することができます。

folium.Marker(["緯度","経度"],"地図の拡大サイズ","地名")

でマーカを地図上にプロットすることが可能です。

folium.CircleMarker(["緯度","経度"],"円の大きさ","地名","円の外枠の色","円内の色")

で円を地図上にプロットすることが可能です。

今回のソースコードはこんな感じ。
地震の大きさで円を色分けしようと思ってIF条件が多くなってしまっていますが、簡単にプロット出来ると思います。
ソースコードが大変見苦しいものになってしまっているのでGW中にきれいにします。

# import libraries
import folium
import pandas as pd
import csv
from datetime import datetime, date, timedelta

#各日にちごとにマップファイルを作成する
def plot_earthquack(data_set):

    save_html_path = "./plot_html/"

    today = datetime.today()
    first_day = today.replace(day=1)
    day = datetime.strftime(first_day, '%Y-%m-%d')
    for k in range(1,500):
        print(day)

        lat1 = []
        lot1 = []
        locate1 = []
        magnitude1 = []
        day1 = []
        for j in range(len(data_set)):
            if str(data_set[j][0]) == str(day):
                lat1.append(float(data_set[j][2]))
                lot1.append(float(data_set[j][3]))
                magnitude1.append(abs(float(data_set[j][5])))
                locate1.append(str(data_set[j][6]))
                day1.append(str(data_set[j][0]))

        # Make a data frame with dots to show on the map
        data = pd.DataFrame({
        'lat':lot1,
        'lon':lat1,
        'name':locate1,
        'mag':magnitude1,
        'day':day1
        })


        #print(data)

        # Make an empty map
        m = folium.Map(location=[35.681236, 139.767125], tiles="Mapbox Bright",  zoom_start=5)

        for i in range(0,len(data)):
            if data.iloc[i]['mag'] >= 1.0 and data.iloc[i]['mag'] < 3.0:
                #folium.Marker([data.iloc[i]['lon'], data.iloc[i]['lat']], popup=data.iloc[i]['name']).add_to(m)

                folium.CircleMarker([data.iloc[i]['lon'], data.iloc[i]['lat']],
                    radius= 5 * data.iloc[i]['mag'],
                    popup=data.iloc[i]['name'],
                    color='#b7fffd',
                    fill_color='#b7fffd'
                    ).add_to(m)

            if data.iloc[i]['mag'] >= 3.0 and data.iloc[i]['mag'] < 5.0:
            #folium.Marker([data.iloc[i]['lon'], data.iloc[i]['lat']], popup=data.iloc[i]['name']).add_to(m)

                folium.CircleMarker([data.iloc[i]['lon'], data.iloc[i]['lat']],
                    radius= 5 * data.iloc[i]['mag'],
                    popup=data.iloc[i]['name'],
                    color='#004cff',
                    fill_color='#004cff'
                    ).add_to(m)

            if data.iloc[i]['mag'] >= 5.0 and data.iloc[i]['mag'] < 6.0:
                #folium.Marker([data.iloc[i]['lon'], data.iloc[i]['lat']], popup=data.iloc[i]['name']).add_to(m)

                folium.CircleMarker([data.iloc[i]['lon'], data.iloc[i]['lat']],
                    radius= 5 * data.iloc[i]['mag'],
                    popup=data.iloc[i]['name'],
                    color='#e66bff',
                    fill_color='#e66bff'
                    ).add_to(m)

            if data.iloc[i]['mag'] >= 6.0 and data.iloc[i]['mag'] < 7.0:
                #folium.Marker([data.iloc[i]['lon'], data.iloc[i]['lat']], popup=data.iloc[i]['name']).add_to(m)

                folium.CircleMarker([data.iloc[i]['lon'], data.iloc[i]['lat']],
                    radius= 5 * data.iloc[i]['mag'],
                    popup=data.iloc[i]['name'],
                    color='#f90081',
                    fill_color='#f90081'
                    ).add_to(m)

            if data.iloc[i]['mag'] >= 7.0:
                #folium.Marker([data.iloc[i]['lon'], data.iloc[i]['lat']], popup=data.iloc[i]['name']).add_to(m)

                folium.CircleMarker([data.iloc[i]['lon'], data.iloc[i]['lat']],
                    radius= 7 * data.iloc[i]['mag'],
                    popup=data.iloc[i]['name'],
                    color='#f91800',
                    fill_color='#f91800'
                    ).add_to(m)

        first_day = first_day - timedelta(days=1)
        day = datetime.strftime(first_day, '%Y-%m-%d')

        # Save it as html
        m.save(save_html_path + str(day) + '.html')

def get_data():
    data = []
    with open('./csv/today.csv',  newline='') as f:
        dataReader = csv.reader(f)
        for row in dataReader:
            data.append(row)
    return data

def main():
    data = get_data()
    plot_earthquack(data)


if __name__ == '__main__':
    main()

以下のように作成したMapオブジェクトはhtml形式として保存することができます。

# Save it as html
m.save(save_html_path + str(day) + '.html')

保存したhtml形式のマップを開くと以下のようにプロット出来ていることがわかると思います。

プロットした各日にちのマップをgifファイルを作る

pngファイルを作成する。

各日にちごとにプロットした地震マップをgifにしたいと思います。
seleniumのwebdriverを使って、html形式のマップのスクリーンショットを撮ることによってpngファイルを作成しました。

driver.save_screenshot("保存先ファイル名")

を利用することによってwebdriverで開くサイトやhtmlのファイルをスクリーンショットを撮ることが出来ます。
今回のソースは以下のようになっています。

import os
import time
from selenium import webdriver
from datetime import datetime, date, timedelta


read_html_path = "./plot_html/"
save_img_path = "./png/"

today = datetime.today()
first_day = today.replace(day=1)
first_day = first_day - timedelta(days=1)
day = datetime.strftime(first_day, '%Y-%m-%d')

driver_path = "./chromedriver.exe"
delay=5

for i in range(1,500):
    fn = str(day)+".html"
    tmpurl='file://{path}/plot_html/{mapfile}'.format(path=os.getcwd(),mapfile=fn)
    print(tmpurl)
    driver = webdriver.Chrome(driver_path)
    driver.get(tmpurl)
    #Give the map tiles some time to load
    time.sleep(delay)

    save = save_img_path + str(day) + ".png"
    driver.save_screenshot(save)
    driver.quit()

    first_day = first_day - timedelta(days=1)
    day = datetime.strftime(first_day, '%Y-%m-%d')

PILを使って、先ほど生産したpngファイル達からgifファイルを生成します。

ImageFont.truetype("Arial Bold.ttf", "文字サイズ")

を利用することによって、gifファイルの中に文字などを入れることが出来ます。
.tffファイルはweb上から入手することが可能で、好きなフォントを使うことが出来ます。

以下今回のソースプログラム

from PIL import Image, ImageDraw, ImageFont
from datetime import datetime, date, timedelta

read_img_path = "./png/"
im = []

today = datetime.today()
first_day = today.replace(day=1)
first_day = first_day - timedelta(days=1)
day = datetime.strftime(first_day, '%Y-%m-%d')

for i in range(1,300):
    png = read_img_path + str(day)+".png"
    print(png)
    img = Image.open(png)
    draw = ImageDraw.Draw(img)
    font = ImageFont.truetype("Arial Bold.ttf", 100)
    draw.text((40, 40),str(day),font=font,fill=(0,0,0))
    im.append(img)
    first_day = first_day - timedelta(days=1)
    day = datetime.strftime(first_day, '%Y-%m-%d')

im[0].save('pillow_imagedraw.gif',
               save_all=True, append_images=im[1:], optimize=False, duration=400, loop=0)

まとめ

foliumを使った地図上のプロットとPILを使ったgifファイルの作り方について簡単に紹介しました。
私自身、業務ではデータやPythonについて一切触れていないので、ご指摘や知見についてなにかあれば教えてください。