さて、本日はpandasでHTMLの表をスクレイピングしてみようと思います。スクレイピングだったらBeautifulSoupやlxmlなどが有名ですが、HTMLの表をスクレイピングするのは割と面倒くさい作業です。「気象庁の震度データベース」というちょうどいいサンプルがあったので、今回はpandasを使ってみようと思います。←これがとっても便利なんです
ここで過去10年間で起きた震度5強以上の地震を検索してみたいと思います。
52件もあることがわかりました。やはり日本は地震大国ですね。そして、地震の検索結果が表になっていることがわかります。ここから緯度経度を取得して地震のポイントをマップにプロットしようと思います。
このページでもポイントをプロットしています。地震の深さでポイントの大きさを変えていますね。
このページをいったんHTMLとしてローカルに保存します（動的に作られているページなので）。そのHTMLの表をスクレイピングしてみます。
インストール
pandasだけでなくlxml、html5lib、BeautifulSoup4もインストールする必要があります。
サンプルコード①
HTMLの表をスクレイピングして結果をCSVに出力するサンプルです。
import pandas as pd
html = 'D:\python\earthquake\earthquake.html'
dataframes = pd.io.html.read_html(html)
table = pd.DataFrame(dataframes[5])
table.columns = ['No','Time','Name','Latitude','Longitude','Depth','Magnitude','SeismicIntensity']
table.to_csv('D:\python\earthquake\earthquake.csv' ,encoding = 'Shift-JIS')
見事に地震データのスクレイピングに成功しました！これだけでできるなんてびっくりですね！！
サンプルコード②
「サンプルコード①」で出力したCSVを元にマップ上にポイントをプロットするサンプルです。色々やっているので、要点だけ書いておきます。
1.CSVから読み込んだ緯度経度を60進数から10進数に変換
2.「1」の緯度経度からポイントをプロット
3.「earthquake」フィーチャクラスにカラム追加(Arcpy)
4.「3」のカラムにCSVから読み込んだ情報を付与(Arcpy)
import arcpy
import pandas as pd
pointGeometryList = []
def main():
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"′")]
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"
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)
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()
ポイントがプロットされました。位置もあっているっぽいですね。
ポイントに属性も入っていますね。
次にシンボルの設定をします。個別値として設定します。←ここもArcpyでできるようにESRIさんお願いします。
サンプルコード③
震度でポイントの色分けを行うサンプルです。
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ではできない模様。うー貧弱すぎる・・・
とりあえずArcMapで変更してみました。いい感じになりました。このようにデータを引っこ抜けば色々な角度から分析できるので便利ですね！
意外なことに都心部では過去10年間に震度5強以上の地震は起きてないようですね。また、東日本大震災の震源はけっこう陸から離れているのにあの津波が起きたわけですから、規格外の大きさの地震ということがわかります。
結果を見ると全国的に大きな地震は起きているのですが、西日本は少ないことがわかりました。沖縄に関しては0です。この結果にプレートとの関係性を絡めたらもっと面白い結果が得られそうです。
とりあえず、今回はここまでとして、次回以降はさらに深く掘り下げた分析をしてみたいと思います。
以上、本日は「 【GISとPython】pandasで気象庁の震度データベースをスクレイピングして震源をマップにプロットする」でした。
