読者です 読者をやめる 読者になる 読者になる

GIS奮闘記

現役GISエンジニアの技術紹介ブログ。主にPythonを使用。

【GISとPython】pandasで気象庁の震度データベースをスクレイピングして震源をマップにプロットする

さて、本日はpandasでHTMLの表をスクレイピングしてみようと思います。スクレイピングだったらBeautifulSoupやlxmlなどが有名ですが、HTMLの表をスクレイピングするのは割と面倒くさい作業です。「気象庁の震度データベース」というちょうどいいサンプルがあったので、今回はpandasを使ってみようと思います。←これがとっても便利なんです

ここで過去10年間で起きた震度5強以上の地震を検索してみたいと思います。
f:id:sanvarie:20160125080801p:plain

52件もあることがわかりました。やはり日本は地震大国ですね。そして、地震の検索結果が表になっていることがわかります。ここから緯度経度を取得して地震のポイントをマップにプロットしようと思います。
f:id:sanvarie:20160125080957p:plain

このページでもポイントをプロットしています。地震の深さでポイントの大きさを変えていますね。
f:id:sanvarie:20160125082714p:plain

このページをいったんHTMLとしてローカルに保存します(動的に作られているページなので)。そのHTMLの表をスクレイピングしてみます。

インストール

pandasだけでなくlxml、html5lib、BeautifulSoup4もインストールする必要があります。

サンプルコード①

HTMLの表をスクレイピングして結果をCSVに出力するサンプルです。

# -*- coding: utf-8 -*-
import pandas as pd

#保存したHTML
html = 'D:\python\earthquake\earthquake.html'

#HTMLを読込
dataframes = pd.io.html.read_html(html)

#表の部分を取得
table = pd.DataFrame(dataframes[5])
table.columns = ['No','Time','Name','Latitude','Longitude','Depth','Magnitude','SeismicIntensity']

#CSV出力
table.to_csv('D:\python\earthquake\earthquake.csv' ,encoding = 'Shift-JIS')

見事に地震データのスクレイピングに成功しました!これだけでできるなんてびっくりですね!!
f:id:sanvarie:20160125083703p:plain

データ準備

まず、以前と同じようにESRIから全国地図をダウンロードします。
全国市区町村界データ | 製品 | ESRIジャパン

ArcMapで読み込みます。
f:id:sanvarie:20160126122232p:plain

ポイントのフィーチャークラスを作成します。
f:id:sanvarie:20160126122349p:plain

フィーチャークラスの属性を以下の地理座標系にします。
f:id:sanvarie:20160126122421p:plain

サンプルコード②

「サンプルコード①」で出力したCSVを元にマップ上にポイントをプロットするサンプルです。色々やっているので、要点だけ書いておきます。

1.CSVから読み込んだ緯度経度を60進数から10進数に変換
2.「1」の緯度経度からポイントをプロット
3.「earthquake」フィーチャクラスにカラム追加(Arcpy)
4.「3」のカラムにCSVから読み込んだ情報を付与(Arcpy)

# -*- coding: utf-8 -*-
import arcpy
import pandas as pd

pointGeometryList = []
def main():

    #CSV読込
    table = pd.read_csv(r"D:\python\earthquake\earthquake.csv",encoding="SHIFT-JIS")

    data_list = []
    for key,rowD in table.iterrows():
        if key > 0 :
            #緯度経度を度分秒に分解
            latDo = rowD.Latitude[0:rowD.Latitude.find(u"°")]
            lonDo = rowD.Longitude[0:rowD.Longitude.find(u"°")]
            latFun = rowD.Latitude[rowD.Latitude.find(u"°")+1:rowD.Latitude.find(u".")]
            lonFun = rowD.Longitude[rowD.Longitude.find(u"°")+1:rowD.Longitude.find(u".")]
            latByou = rowD.Latitude[rowD.Latitude.find(u".")+1:rowD.Latitude.find(u"′")]
            lonByou = rowD.Longitude[rowD.Longitude.find(u".")+1:rowD.Longitude.find(u"′")]

            #緯度経度を60進数から10進数に変換
            lat,lon = latLong(latDo,lonDo,latFun,lonFun,latByou,lonByou)
            data_list.append([rowD.No,rowD.Time,rowD.Name,lat,lon,rowD.Depth,rowD.Magnitude,rowD.SeismicIntensity])
            CreatePont(lat,lon)
    table2 = pd.DataFrame(data_list)
    table2.columns = ["No","Time","Name","Latitude","Longitude","Depth","Magnitude","SeismicIntensity"]
    arcpy.CopyFeatures_management(pointGeometryList, r"C:\ArcPySample\test.gdb\earthquake")

    arcpy.env.workspace = ur"C:\ArcPySample\test.gdb"

    #earthquakeレイヤにカラムを追加
    arcpy.AddField_management("earthquake", "Time", "Text")
    arcpy.AddField_management("earthquake", "Name", "Text")
    arcpy.AddField_management("earthquake", "lat", "Double")
    arcpy.AddField_management("earthquake", "lon", "Double")
    arcpy.AddField_management("earthquake", "Depth", "Text")
    arcpy.AddField_management("earthquake", "Magnitude", "Text")
    arcpy.AddField_management("earthquake", "SIntensity", "Text")

    for key,rowD in table2.iterrows():
        cursorJ = arcpy.UpdateCursor("earthquake")
        for rowJ in cursorJ:
            if rowD.No == rowJ.OBJECTID:
                rowJ.setValue("Time", rowD.Time)
                rowJ.setValue("Name", rowD.Name)
                rowJ.setValue("lat", rowD.Latitude)
                rowJ.setValue("lon", rowD.Longitude)
                rowJ.setValue("Depth", rowD.Depth)
                rowJ.setValue("Magnitude", rowD.Magnitude)
                rowJ.setValue("SIntensity", rowD.SeismicIntensity)
                cursorJ.updateRow(rowJ)

#ポイントを配置
def CreatePont(lat,lon):
    global pointGeometryList

    point = arcpy.Point()
    point.X = lon
    point.Y = lat

    #日本地図の空間参照を利用
    dataset = r"D:\python\earthquake\japan_ver80.shp"
    spatial_reference=arcpy.Describe(dataset).spatialReference
    pointGeometry = arcpy.PointGeometry(point,spatial_reference)
    pointGeometryList.append(pointGeometry)

#緯度経度を60進数から10進数に変換
def latLong(latDo,lonDo,latFun,lonFun,latByou,lonByou):
    lat1 = float(latDo)
    lat2 = float(latFun) / 60.0
    if type(latByou) == float:
        lat3 = float(latByou) / 60.0 / 60.0
    else:
        lat3 = float(0)
    lat1 = lat1 + lat2 + lat3

    long1 = float(lonDo)
    long2 = float(lonFun) / 60.0
    if type(lonByou) == float:
        long3 = float(lonByou) / 60.0 / 60.0
    else:
        long3 = float(0)
    long1 = long1 + long2 + long3
    return lat1,long1

if __name__ == '__main__':
    main()

ポイントがプロットされました。位置もあっているっぽいですね。
f:id:sanvarie:20160126132631p:plain

ポイントに属性も入っていますね。
f:id:sanvarie:20160126141215p:plain

次にシンボルの設定をします。個別値として設定します。←ここもArcpyでできるようにESRIさんお願いします。
f:id:sanvarie:20160126141510p:plain

サンプルコード③

震度でポイントの色分けを行うサンプルです。

# -*- coding: utf-8 -*-
import arcpy
mxd = arcpy.mapping.MapDocument("current")
lyr = arcpy.mapping.ListLayers(mxd, "earthquake")[0]
if lyr.symbologyType == "UNIQUE_VALUES":
  lyr.symbology.valueField = "SIntensity"
  lyr.symbology.addAllValues()
arcpy.RefreshActiveView()
arcpy.RefreshTOC()
del mxd

ポイントの色分けができました。次にシンボルの大きさを分けようと思ったのですが、どうやらこれもArcpyではできない模様。うー貧弱すぎる・・・
f:id:sanvarie:20160126141932p:plain

とりあえずArcMapで変更してみました。いい感じになりました。このようにデータを引っこ抜けば色々な角度から分析できるので便利ですね!
f:id:sanvarie:20160126142201p:plain

意外なことに都心部では過去10年間に震度5強以上の地震は起きてないようですね。また、東日本大震災震源はけっこう陸から離れているのにあの津波が起きたわけですから、規格外の大きさの地震ということがわかります。
f:id:sanvarie:20160126155433p:plain

結果を見ると全国的に大きな地震は起きているのですが、西日本は少ないことがわかりました。沖縄に関しては0です。この結果にプレートとの関係性を絡めたらもっと面白い結果が得られそうです。
とりあえず、今回はここまでとして、次回以降はさらに深く掘り下げた分析をしてみたいと思います。

以上、本日は「 【GISPython】pandasで気象庁の震度データベースをスクレイピングして震源をマップにプロットする」でした。