【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