Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationEventAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
19
Help us understand the problem. What are the problem?

More than 1 year has passed since last update.

posted at

updated at

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

はじめに

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

pillow_imagedraw2.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形式のマップを開くと以下のようにプロット出来ていることがわかると思います。
無題.png

プロットした各日にちのマップを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について一切触れていないので、ご指摘や知見についてなにかあれば教えてください。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
19
Help us understand the problem. What are the problem?