Pythonで自分だけの世界地図を作ってみる ～オープンデータを自分の欲しい形にしてみよう～

本日は世界地図を作ってみようと思います。オープンデータ関連で検索すれば世界地図がダウンロードできるサイトがいくつか検索できます。ただし、すべて海外のサイトのため、表記が英語になっていて非常に見づらいです。なので、今回は日本人向けの世界地図を作ってみようと思います。

※今回使用したソースと作成したGDBはGitHubに公開しています。
GitHub - sanvarie/MinnanoArcGIS

ダウンロード

このサイトからShapeファイルをダウンロードします。
http://www.naturalearthdata.com/

ダウンロードページ＞Large scale data, 1:10m＞Cultural　からダウンロードします。

こんな感じの表記なので、うーん、といった感じでしょうか。ただ、グルジアがジョージアになっていたりと、データ自体は新しそうです。

実行環境

• Windows7 64bit
• ArcGIS10.2
• Python2.7

事前準備

• 「Map.gdb」を作成します。

ポリゴンを格納するフィーチャクラスを作成する

ダウンロードした「ne_10m_admin_0_countries.shp」を「Map.gdb」にコピーします。

サンプルコード①

ダウンロードした「ne_10m_admin_0_countries.shp」を「Map.gdb」にコピーして、不要なカラムを削除するサンプルです。また、正式名称がないフィーチャに関しては名称で置き換えます（クック諸島とか）。

■copyWorldMap.py

# -*- coding: utf-8 -*-
import arcpy

#対象GDB
arcpy.env.workspace = "C:\ArcPySample\Map.gdb"

#コピーするフィーチャクラス名
inFeature = "WorldMap"

#ShapeをMap.gdbにコピー
arcpy.CopyFeatures_management(r"D:\python\WorldMap\ne_10m_admin_0_countries.shp", inFeature)

#コードブロック
codeblock = """
def UpdateColumn(FORMAL_EN,NAME_LONG):
    if FORMAL_EN == " ":
        return NAME_LONG
    else:
        return FORMAL_EN
"""

# 条件式を設定
expression = "UpdateColumn(!FORMAL_EN!,!NAME_LONG!)"

#フィールド演算(正式名称が半角スペースになっているものをNAME_LONGに置き換える)
arcpy.CalculateField_management(inFeature, "FORMAL_EN", expression, "PYTHON_9.3", codeblock)

#カラム追加
arcpy.AddField_management(inFeature, "NAME_EN", "Text",field_length = 60,field_alias="英語名")

#フィールド演算(FORMAL_ENをNAME_ENに置き換える)
arcpy.CalculateField_management(inFeature, "NAME_EN", "!FORMAL_EN!", "PYTHON_9.3")

field_list = []
fields = arcpy.ListFields(inFeature)
for field in fields:
    if field.type != "Geometry":
        if field.name != "NAME_EN" and field.name != "OBJECTID" and field.name != "Shape_Length" and field.name != "Shape_Area":
            field_list.append(field.name)

#不要なフィールドを削除
arcpy.DeleteField_management(inFeature,
                             field_list)

こんな感じになれば成功です。

世界地図に属性を付与する

属性カラムがすっきりしたところで、新たに属性を付与させようと思います。

対象サイト

外務省のサイトがよさげだったので今回はこのサイトから情報を取得しようと思います。世界の国々 ｜ 外務省
ただし、世界と日本のデータを見る（世界の国の数，国連加盟国数，日本の大使館数など） ｜ 外務省で確認すると、「世界の国の数＝196か国。これは，現在，日本が承認している国の数である195か国に日本を加えた数です。」とあるようにフィーチャ数（255）と一致しません。残念ですが、一致するものだけに属性を付与したいと思います。

また、「独立年」がうまくとってこれないので、ここはカットしようと思います(泣)。Pandasでスクレイピングする限界ですかね・・・

サンプルコード②

必要なカラムを追加するサンプルです。

■addColumn.py

# -*- coding: utf-8 -*-
import arcpy

#対象GDB
arcpy.env.workspace = "C:\ArcPySample\Map.gdb"

#コピーするフィーチャクラス名
inFeature = "WorldMap"

#カラム追加
arcpy.AddField_management(inFeature, "NAME", "Text",field_length = 50,field_alias="国名")
arcpy.AddField_management(inFeature, "CAPITAL", "Text",field_length = 50,field_alias="首都")
arcpy.AddField_management(inFeature, "LANGUAGE", "Text",field_length = 50,field_alias="主要言語")
arcpy.AddField_management(inFeature, "AREASQUARE", "Double",field_length = 50,field_alias="面積（1,000平方キロ）")
arcpy.AddField_management(inFeature, "POPULATION", "Double",field_length = 50,field_alias="人口（100万人）")
arcpy.AddField_management(inFeature, "CURRENCY", "Text",field_length = 50,field_alias="通貨単位")

サンプルコード③

上述したサイトから各国の情報を取得して「WorldMap」フィーチャクラスに属性を付与するサンプルです。ものすごい力技になってしまいましたが。。。

■setCountryData.py

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

def updateAttribute(row,cur,conTable):
    row.setValue("NAME", conTable[0])
    row.setValue("CAPITAL", conTable[2])
    row.setValue("LANGUAGE", conTable[4])

    AreaSquere = conTable[5]
    #変な文字列がまじっているので
    if isinstance(AreaSquere, float):
        row.setValue("AREASQUARE", AreaSquere)
    elif isinstance(AreaSquere, long):
        row.setValue("AREASQUARE", AreaSquere)
    #なぜかdatetimeとして認識されるものがあるのでそれは無視
    elif isinstance(AreaSquere, datetime.date):
        pass
    else:
        if AreaSquere.find(u"平方キロ") > -1:
            row.setValue("AREASQUARE", AreaSquere[0:AreaSquere.find(u"平方キロ")])

    population = conTable[6]
    #変な文字列がまじっているので
    if isinstance(population, float):
        row.setValue("POPULATION", population)
    else:
        if population.find(u"約") > -1:
            if population.find(u"人") > -1:
                population = population[population.find(u"約")+1:]
                row.setValue("POPULATION", population[0:population.find(u"人")])
            else:
                row.setValue("POPULATION", population[population.find(u"約")+1:])
        elif population.find(u"人") > -1:
            row.setValue("POPULATION", population[0:population.find(u"人")])

    currency = conTable[7]
    if currency.find(u"（") > 1:
        row.setValue("CURRENCY", currency[0:currency.find(u"（")])
    else:
        row.setValue("CURRENCY", currency)

    cur.updateRow(row)

tableList = []
countryList = []

#対象のURL
countryList.append("http://www.mofa.go.jp/mofaj/kids/ichiran/i_europe.html")
countryList.append("http://www.mofa.go.jp/mofaj/kids/ichiran/i_africa.html")
countryList.append("http://www.mofa.go.jp/mofaj/kids/ichiran/i_chuto.html")
countryList.append("http://www.mofa.go.jp/mofaj/kids/ichiran/i_asia.html")
countryList.append("http://www.mofa.go.jp/mofaj/kids/ichiran/i_oceania.html")
countryList.append("http://www.mofa.go.jp/mofaj/kids/ichiran/i_n_america.html")
countryList.append("http://www.mofa.go.jp/mofaj/kids/ichiran/i_latinamerica.html")

#対象のフィーチャクラス
arcpy.env.workspace = "C:\ArcPySample\Map.gdb"
spatial_reference=arcpy.SpatialReference(4326)

for l in countryList:
    #HTMLを読込
    df = pd.io.html.read_html(l)
    tableList.append(df[0])

#データフレームを結合
conTable = pd.concat(tableList, ignore_index=True)
conTable.columns = ["NAME","ENGLISH_NAME","CAPITAL","INDE_YEAR","LANGUAGE","AREASQUARE","POPULATION","CURRENCY"]

for i in conTable.index:
    dfEnglishName = conTable.ix[i][1]
    cursor = arcpy.UpdateCursor("WorldMap")
    for row in cursor:
        if dfEnglishName == row.NAME_EN:
            updateAttribute(row,cursor,conTable.ix[i])
        #色々力技
        elif dfEnglishName.find("of") > -1:
            if dfEnglishName.find("Australia") > -1 and row.NAME_EN.find("Australia") > -1:
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName.find("Vatican") > -1 and row.NAME_EN.find("Vatican") > -1:
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName.find("Hungary") > -1 and row.NAME_EN.find("Hungary") > -1:
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName.find("Verde") > -1 and row.NAME_EN.find("Verde") > -1:
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName.find("Gambia") > -1 and row.NAME_EN.find("Gambia") > -1:
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName.find("Comoros") > -1 and row.NAME_EN.find("Comoros") > -1:
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName.find("Principe") > -1 and row.NAME_EN.find("Principe") > -1:
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName.find("Nepal") > -1 and row.NAME_EN.find("Nepal") > -1:
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName.find("Viet") > -1 and row.NAME_EN.find("Viet") > -1:
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName.find("Bahamas") > -1 and row.NAME_EN.find("Bahamas") > -1:
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName == "Republic of Guinea" and row.NAME_EN == "Republic of Guinea":
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName == "Independent State of Papua New Guinea" and row.NAME_EN == "Independent State of Papua New Guinea":
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName.find("Cote") > -1  and row.NAME_EN == "Republic of Ivory Coast":
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName == "Republic of Congo" and row.NAME_EN == "Republic of Congo":
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName == "Democratic Republic of the Congo" and row.NAME_EN == "Democratic Republic of the Congo":
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName.find("Togo") > -1 and row.NAME_EN == "Togolese Republic":
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName.find("Lebanon") > -1 and row.NAME_EN == "Lebanese Republic":
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName == "Republic of Korea" and row.NAME_EN == "Republic of Korea":
                updateAttribute(row,cursor,conTable.ix[i])
            elif dfEnglishName.find(row.NAME_EN[row.NAME_EN.find("of")+3:]) > -1 and row.NAME_EN.find("of") > -1 \
                 and dfEnglishName != "Republic of Korea" and dfEnglishName != "Republic of Congo" and dfEnglishName != "Independent State of Papua New Guinea" \
                 and dfEnglishName != "Republic of Equatorial Guinea" and dfEnglishName != "Republic of Guinea-Bissau" and dfEnglishName != "Democratic Republic of the Congo":
                updateAttribute(row,cursor,conTable.ix[i])
        elif dfEnglishName.find("Spain") > -1 and row.NAME_EN.find("Spain") > -1:
                updateAttribute(row,cursor,conTable.ix[i])
        elif dfEnglishName.find("Brunei") > -1 and row.NAME_EN.find("Brunei") > -1:
            updateAttribute(row,cursor,conTable.ix[i])
        elif dfEnglishName.find("Nevis") > -1 and row.NAME_EN.find("Nevis") > -1:
            updateAttribute(row,cursor,conTable.ix[i])
        elif row.NAME_EN.find(dfEnglishName) > -1 and row.NAME_EN.find("Kingdom") == -1  and dfEnglishName.find("of") == -1 \
             and row.NAME_EN != "South Georgia and South Sandwich Islands" and row.NAME_EN != "Indian Ocean Territories" and row.NAME_EN != "British Indian Ocean Territory":
            updateAttribute(row,cursor,conTable.ix[i])

こんな感じになりました。面積や人口がうまく入らないところもあったのですが、とりあえずよしとします。。。

サンプルコード④

Shape_Length、Shape_Areaの位置が気持ち悪いのでフィーチャクラスをコピーします。

■copyFeatures.py

# -*- coding: utf-8 -*-
import arcpy

#対象GDB
arcpy.env.workspace = "C:\ArcPySample\Map.gdb"

#ShapeをMap.gdbにコピー
arcpy.CopyFeatures_management("WorldMap", "WorldMap2")

OKですね。後は「WolrMap」フィーチャクラスを削除して「WolrMap2」を「WolrMap」にリネームしてもOKだと思います。

課題

1. 「独立年」がとってこれなかった
2. 「面積」「人口」が一部とってこれなかった
3. 承認していない国の情報がない
4. プログラムが力技すぎ

課題は多いです。ですが、とりあえず、日本が承認している国の日本語名を入れることができたので最低限のことはできたかと思います。台湾がないのは中国との政治的なあれですかね。今後、この世界地図を使って色々情報を追加していきたいですね。「サッカー日本代表と各国との対戦成績を分析」みたいに。

以上、今回は「Pythonで自分だけの世界地図を作ってみる　～オープンデータを自分の欲しい形にしてみよう～」でした。

【PythonとGIS】GIS的視点で高校サッカーを分析してみる

さて、とうとうサッカーネタに手を出してみようと思います。最初は高校サッカーにしてみました。「どの県が最強なのか」というのを分析・視覚化してみようと思います。

※ソースと作成したデータはGitHubで公開しています。https://github.com/sanvarie/MinnanoArcGIS

分析の詳細

• 各年の選手権の成績において、各校にそれぞれポイントを付与します(以下参照)。

優勝 ・・・10ポイント
準優勝 ・・・5ポイント
ベスト4・・・3ポイント

• 上記のポイントを県ごとに集計してポイント数で順位付を行います。

• 集計は戦後からとします。

主な技術的要素

1. スクレイピング　　　・・・成績表から高校名をスクレイピングします。今回はこのサイトを対象にしてみました。全国高校サッカー選手権 歴代優勝校(ベスト4)
2. ジオコーディング　　・・・スクレイピングした高校名からジオコーディングをしてみます。ジオコーディングについてはこちらをご参照ください。Pythonでジオコーディングをやってみる - GIS奮闘記
3. ポイントのプロット　・・・ジオコーディングした緯度経度でポイントをプロットします。
4. 集計　　　　　　　　・・・各県ごとにポイントの集計（上述した「分析の詳細」参照）を行います。

バージョン

ArcGIS10.2、Python2.7

インストール

Pandas、geocoderのインストールをお願いします。また、lxml、html5lib、BeautifulSoup4も必要ですのでインストールをお願いします。

事前準備

1. 「Soccer.gdb」を作成します。
2. 前回作成した日本地図を「Soccer.gdb」にコピーします（フィーチャクラス名は「Japan」）。
3. フィーチャクラス「Japan」に「Point」カラムを追加します（型はLong Integer）。
4. 「HighShcool」フィーチャクラスを作成します。

事前準備完了後のイメージです。「Japan」フィーチャクラスには日本地図が格納されており、「HighShcool」フィーチャクラスは空の状態です。

「HighShcool」フィーチャクラスにポイントをプロット

ArcPyを使用して「HighShcool」フィーチャクラスにポイントをプロットします。

サンプルコード

無理矢理感は否めませんが、なんとかポイントをプロットすることができました。

■highSchoolSoccer.py

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

def getCoordinate(location_name):
    try:
        #地名から座標を取得する
        ret = geocoder.google(location_name)
    except KeyError, e:
        print "KeyError(%s)" % str(e)
        return

    return ret

def createPoint(name,points,year):

    #ジオコーディング
    loc = getCoordinate(name)

    if loc.lat is not None:
        point = arcpy.Point()
        point.X = loc.lng
        point.Y = loc.lat
        pointGeometry = arcpy.PointGeometry(point,spatial_reference)

        cur = arcpy.da.InsertCursor(infc, ["SHAPE@","Name","Point","Year"])

        #校名のあとに県などがついている場合、分解する
        if name.find(" ") > 0:
            if len(name.split(" ")) == 2:
                name,ken = name.split(" ")
            if len(name.split(" ")) == 3:
                name,ken,nihon = name.split(" ")
        cur.insertRow((pointGeometry,name,points,year))
        del cur
    else:
        print name + u"の座標がとれない"

        #多々良学園などの座標がとれないのでこれで対応
        if name != "":
            name = name + " " + u"日本"
            #浦和南 埼玉 日本とかだとだめなのでこの場合は南浦和にする
            if len(name.split(" ")) == 4:
                name,ken,nihon,nihon2 = name.split(" ")
                createPoint(name,points,year)
            else:
                createPoint(name,points,year)

#対象のフィーチャクラス
infc = "C:\ArcPySample\Soccer.gdb\HighShcool"
spatial_reference=arcpy.SpatialReference(4612)

#対象サイト
html = 'http://www.tigerkaz.info/highschool/senshuken.html'

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

#表の部分を取得
table = pd.DataFrame(dataframes[0])

#カラム作成
table.columns = ['LINK','Time','Year','Champion','Prefecture','Time2','Finalist','Prefecture2','Best4','Prefecture3','Best4_2','Prefecture4']

schoolList = []

for key,row in table.iterrows():

    #年度を保持
    if row.Year.year != -1:
        year = row.Year.year

    #戦後の結果のみを対象とする
    if row.Year.year == datetime.datetime(1945,1,1).year:
        break

    if row.Champion == u"優勝":
        continue

    #両校優勝の場合、変なとこに優勝校の一つが入っているので
    if isinstance(row.Champion, float):
        yusho = row.LINK

        if yusho.find("(") > 0:

            yusho = yusho[0:yusho.find("(")]

        #帝京とかだと中国にジオコーディングされてしまうので
        yusho = yusho  + " " + u"日本"

        schoolList.append([yusho,"","","",year])
    else:
        #校名 + 県名でジオコーディングする。国見とかだと変なとこに飛ぶので。
        #ただし、これをやると高校から微妙に座標がずれる。が、とりあえず近くなのでよしとする。
        yusho = row.Champion   + " " + row.Prefecture

        if yusho.find("(") > 0:
            yusho = yusho[0:yusho.find("(")] + " " + row.Prefecture

        #両校優勝の場合列がずれているので
        if row.Finalist == u"(両校優勝)":
            best4 = row.Prefecture2 + " " + row.Best4
            best4_2 = row.Prefecture3 + " " + row.Best4_2
            schoolList.append([yusho,"",best4,best4_2,row.Year.year])
        else:
            junYusho = row.Finalist + " " + row.Prefecture2
            best4 = row.Best4 + " " + row.Prefecture3
            best4_2 = row.Best4_2 + " " + row.Prefecture4
            schoolList.append([yusho,junYusho,best4,best4_2,row.Year.year])

#リストをデータフレームに変換
schooDf = pd.DataFrame(schoolList)
schooDf.columns = ['Champion','Finalist','Best4','Best4_2','Year']

for key,rowS in schooDf.iterrows():
    #ポイントをプロット
    createPoint(rowS.Champion,10,rowS.Year)
    createPoint(rowS.Finalist,5,rowS.Year)
    createPoint(rowS.Best4,3,rowS.Year)
    createPoint(rowS.Best4_2,3,rowS.Year)

実行後、それっぽいところにポイントがプロットされているのがわかります。ただ、ジオコーディングの関係上（コメント参照）、本当の位置からは若干ずれている気がしますが、県からはずれていないのでよしとします。

「Japan」フィーチャクラスに各校の成績を反映

ArcPyを使用して「Japan」フィーチャクラスの「Point」カラムに各校のポイントを集計した値を付与します。

サンプルコード

■updatePoint.py

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

arcpy.env.workspace = "C:\ArcPySample\Soccer.gdb"

#インターセクト
highShcool = arcpy.Intersect_analysis(["HighShcool", "Japan"],"HighShcool_Japan" ,"", "", "")

field_list = []
fields = arcpy.ListFields(highShcool)
for field in fields:
    if field.type != "Geometry":
        field_list.append(field.name)

#属性テーブルをPandasに格納
column = arcpy.da.FeatureClassToNumPyArray(highShcool,field_list,null_value=-9999)
df = pd.DataFrame(column)

#ポイントの集計
df_group = df.groupby('KEN')['Point'].sum()

#集計値をJapanに格納
for i in df_group.index:
    cursor = arcpy.UpdateCursor("Japan")
    for row in cursor:
        if i == row.KEN:
            row.setValue("Point", df_group.ix[i])
            cursor.updateRow(row)
    del cursor


#コードブロック
codeblock = """
def UpdatePoint(Point):
    if Point is None:
        return 0
    else:
        return Point
"""

# 条件式を設定
expression = "UpdatePoint(!Point!)"

#ポイントが0の県に対してフィールド演算
arcpy.CalculateField_management("Japan", "Point", expression, "PYTHON_9.3", codeblock)

このような結果になれば成功です。各県のポイントが算出できたら何かテンションあがってきました！この後の順位付が楽しみですね♪

「Japan」フィーチャクラスの「Point」カラムによるレンダリング

まずは「Japan」フィーチャクラスのシンボル設定を以下のようにします。OKを押すとマップ上から日本地図が消えますが問題ありません。

サンプルコード

各県ごとに付与されたポイントを以下レンジでレンダリングを行います。数値が大きいほど、その県が強いということがわかりますね。
0
1～10
11～30
31～50
51～100
101～150
151～200

■rendering.py

# -*- coding: utf-8 -*-
import arcpy

mxd = arcpy.mapping.MapDocument("CURRENT")
for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.symbologyType == "GRADUATED_COLORS":
        if lyr.symbology.valueField == "":
            lyr.symbology.valueField = "Point"
        lyr.symbology.classBreakValues = [0,10,30,50,100,150,200]
        arcpy.RefreshActiveView()

このような感じになりました。

結果の分析

1. 都市圏が強い。有名私立が多いですし、順当な結果といったところでしょうか。
2. 雪国の成績は今一歩。やはり雪の影響で冬は外でトレーニングをするのが難しいからでしょうか。
3. 強豪校がある県が強い。当たり前ですが。国見がある長崎は94ポイントを獲得しています。また、雪国ですが、石川県は星陵だけで21ポイント、秋田県は秋田商業だけで34ポイントを獲得していますね。
4. 2校出場できる東京の順位が思ったより高くない。89ポイントを獲得して全体で7位ですが、有名高校の数、毎年2校出場の割には少ないような気もします。
5. 沖縄、香川、佐賀、山形、新潟、長野、鳥取は残念ながら0ポイント。確かにこの辺の高校はあまり印象に残っていない気がします。
6. 最強の県は静岡。やはりサッカー王国と言われるだけあります。しかし、近年の成績はあまりふるわなく、1995年に静岡学園が優勝して以来、優勝からは遠ざかっていますね。おそらく直近20年の分析結果では首位ではないような気がします。ちなみに各県の順位は以下のようになりました。

課題

1. 関東などの地域ごとのポイント集計も行えるようにした方がより分析が面白くなりそうですね。
2. 優勝、準優勝、ベスト4といったようなおおざっぱなポイント付与ではなく、一勝ごとにポイントを付与するなど、もう少し細かい分析が必要。ただ、そのデータをどこから手に入れるのかという問題がありますが・・・

☆番外編☆

せっかくなので、高校ごとの順位もつけてみようと思います。

テーブルの作成

「HighSchool_Record」テーブルを作成します。

作成後、GDBの中身は以下のようになります。「HighShcool_Japan」は「updatePoint.py」でインターセクトした結果です。

サンプルコード

高校ごとのポイントを集計した結果をテーブルに格納します。

■highSchoolRecord.py

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

arcpy.env.workspace = "C:\ArcPySample\Soccer.gdb"

field_list = []
fields = arcpy.ListFields("HighShcool")
for field in fields:
    if field.type != "Geometry":
        field_list.append(field.name)

#属性テーブルをPandasに格納
column = arcpy.da.FeatureClassToNumPyArray("HighShcool",field_list,null_value=-9999)
df = pd.DataFrame(column)

#ポイントの集計
df_group = df.groupby('Name')['Point'].sum()

#ポイントの降順でソート
df_group.sort("Point",ascending=False)

cursor = arcpy.InsertCursor("HighSchool_Record")

#集計値を「HighSchool_Record」に格納
for i in df_group.index:
    row = cursor.newRow()
    row.setValue("Name", i)
    row.setValue("Point", df_group.ix[i])
    cursor.insertRow(row)
del cursor

結果を見てみると、帝京が見事に一位を獲得しました！最近はあまり元気がありませんが、やはり名門ですね。

全体の順位は以下のようになっています。有名選手の出身高はこの中でも上位が多い気がします。帝京＝田中達也、国見＝大久保嘉人、藤枝東＝長谷部誠、市立船橋＝カレン・ロバート、浦和市立＝ごめんなさい、わからない・・・、鹿児島実業＝遠藤保仁、韮崎＝中田英寿、などなど上位は選手がぱっと浮かぶ高校が多いですね。

色々やってみて最後に気づいたことがあります。今年の結果がない(笑)。最新の結果を反映したら少し違った結果になるかと思います。今回は単純に選手権の成績だけでしたが、人口密度や高校数を加味した分析などもおもしろそうですね。また、Jリーグなどでも今回とは違った分析が行えそうだと思います。例えば「Jリーグの順位における対象地域への経済効果の差について」のような感じですかね。それをやるには自分自身もっとスキルを上げる必要があるので、頑張りたいと思います。

以上、今回は「【PythonとGIS】GIS的視点で高校サッカーを分析してみる」でした。

PythonでESRIの日本地図(全国市区町村界データ)を加工してみる ～オープンデータを自分の欲しい形にしてみよう～

本日は当ブログでもよく使用しているESRIさんの日本地図を加工してみようと思います。このデータは全国の市区町村ごとのポリゴンデータなのですが、県ごとのポリゴンが欲しいという方もいるかと思います。単純に県ごとにするだけなら簡単なのですが、人口や世帯数といったデータも一緒に集約することはArcMapではできない？！のではないでしょうか。←自分がやり方を知らないだけかもですが・・・

また、今回のソースと加工したデータはGitHubで公開しました。
https://github.com/sanvarie/MinnanoArcGIS

URL

以下サイトでダウンロード可能です。
全国市区町村界データ | 製品 | ESRIジャパン

仕様

ただのコピペですが。

データ仕様

・フォーマット シェープファイル
・図形タイプ ポリゴン
・座標系 地理座標系 日本測地系 2000（緯度経度 JGD 2000）

属性情報

フィールド名
JCODE・・・市区町村コード
KEN・・・都道府県名
SICHO ・・・支庁名・振興局名
GUN・・・郡名（町村部のみ）
SEIREI・・・政令指定都市の市名
SIKUCHOSON・・・市区町村名
CITY_ENG・・・市区町村名（英語）
P_NUM・・・人口
H_NUM・・・世帯数

県ごとのポリゴンを作成

ディゾルブを行います。以下スクリプトを実行してください。また、事前に「ArcPyJapan.gdb」というgdbを作成しておいてください。ちなみにここまでは【Pythonで分析】ArcpyとPandasを使用して将来推計人口を視覚化する - GIS奮闘記でやりましたね。よく考えたら今回はほとんどこの時と同じようなことをしているような・・・でもまぁ気にしません(笑)

サンプルコード

dissolve.py

# -*- coding: utf-8 -*-
import arcpy

#ディゾルブ
arcpy.Dissolve_management("D:\python\soccer\japan_ver80.shp", "C:\ArcPySample\ArcPyJapan.gdb\Japan",
                           "KEN", "", "MULTI_PART",
                           "DISSOLVE_LINES")

結果を確認します。県ごとのポリゴンが完成しました。ただ、このままだと人口、世帯数がありませんね。

カラム追加

上記で作成したポリゴンに人口、世帯数カラムを追加します。

サンプルコード

addcolumn.py

# -*- coding: utf-8 -*-
import arcpy

arcpy.env.workspace = "C:\ArcPySample\ArcPyJapan.gdb"
arcpy.AddField_management("Japan", "P_NUM", "Long")
arcpy.AddField_management("Japan", "H_NUM", "Long")

カラムが追加されましたね。

県ポリゴンに人口、世帯数を与える

上記で追加した人口、世帯数カラムに値を付与します。

インストール

Pandasのインストールをお願いします。

サンプルコード

updatejapanmap.py

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

#日本地図のShape
inFeatures = "D:\python\soccer\japan_ver80.shp"

#更新するフィーチャクラスがあるgdb
arcpy.env.workspace = "C:\ArcPySample\ArcPyJapan.gdb"

field_list = []
for field in arcpy.ListFields(inFeatures):
    if field.type != "Geometry":
        field_list.append(field.name)

df = pd.DataFrame(arcpy.da.FeatureClassToNumPyArray(inFeatures,field_list,null_value=-9999))

#グルーピング
df_group = df.groupby('KEN')['P_NUM','H_NUM'].sum()

for key,row in df_group.iterrows():
    cursorJ = arcpy.UpdateCursor("Japan")
    for rowJ in cursorJ:
        if key == rowJ.KEN:
            rowJ.setValue("P_NUM", row.P_NUM)
            rowJ.setValue("H_NUM", row.H_NUM)
            cursorJ.updateRow(rowJ)

人口、世帯数を持った県ごとのポリゴンの作成に成功しました。

最近、オープンデータの数が増加している気がしますが、自分の求めているものとちょっと違っていたりする経験がある方が多いのではないでしょうか。こういった技術を使えば、オープンデータを自分の欲しかったデータに加工することが可能です(今回はかなり簡単な例でしたが)。

以上、本日は「PythonでESRIの日本地図(全国市区町村界データ)を加工してみる」でした！

ArcPyを使用してフィーチャクラスの情報をCSV出力する

さて、本日は「ArcPyを使用してフィーチャクラスの情報をCSV出力する」です。実際に業務では数十個のフィーチャクラスがあり、それぞれにサブタイプなどを設定していたりします。例えば、サブタイプごとのアイテム数を抽出したい場合、一つ一つ確認していくのは現実的ではありません。そういった場合に対応するPythonスクリプトを書いてみました。恥ずかしながら最近Gistを使い始めたので、そこにもアップしてみました。
フィーチャクラスのサブタイプごとのアイテム数を取得するスクリプト（複数のGDBには対応していません） · GitHub

データ

以下のようなGDBとフィーチャクラスを作成しました。

各フィーチャクラスのサブタイプ設定です。PIPEフィーチャクラスはサブタイプ設定なしにしました。
■BOUNDARY

■CHARACTER

■LANDMARK

■PIPE

■ROAD

インストール

ArcMapがインストールされていればNumpyはインストール済だと思いますので、Pandasのインストールをお願いします。

サンプルコード

各フィーチャクラスもしくはフィーチャクラスのサブタイプごとのアイテム数をCSV出力するPythonスクリプトです。

以下の項目を出力します。
・gdb
・レイヤ名
・エイリアス
・ジオメトリタイプ
・サプタイプ
・アイテム数

# -*- coding: utf-8 -*-
import arcpy
from os.path import join
import pandas as pd
import numpy as np

def getInfo(dataSource,fList):

    #サブタイプを取得
    subtypefields = arcpy.da.ListSubtypes(dataSource)
    subDict = {}
    for name in subtypefields.keys():
        groupDict = subtypefields[name]
        for stkey in list(groupDict.keys()):
            if stkey == 'SubtypeField':

                fields = groupDict[stkey]

                #サブタイプが設定されていない場合
                if fields == "":
                    outputstring = fList[0] \
                                   + "," + fList[1] \
                                   + "," + fList[2] \
                                   + "," + fList[3] \
                                   + "," + "" \
                                   + "," + str(fList[4]) + "\n"

                    f.write(outputstring.encode("SHIFT-JIS"))

                #サブタイプが設定されている場合
                else:
                    for stkey in list(groupDict.keys()):
                        if stkey == 'Name':

                            #フィーチャをnumpyに変換してpandasに格納
                            column = arcpy.da.FeatureClassToNumPyArray(dataSource,("SUBTYPE_CD"),null_value=-9999)
                            dfWater = pd.DataFrame(column)

                            fields = groupDict[stkey]
                            subDict[fields] = name

                            #サブタイプごとのデータを抽出
                            dfWaterValue = pd.DataFrame(dfWater[dfWater.SUBTYPE_CD == subDict[fields]])

                            outputstring = fList[0] \
                                           + "," + fList[1] \
                                           + "," + fList[2] \
                                           + "," + fList[3] \
                                           + "," + fields \
                                           + "," + str(len(dfWaterValue)) + "\n"

                            f.write(outputstring.encode("SHIFT-JIS"))

# ファイルの出力先を指定
f = open(r"D:\python\featureclass/list.csv", "w")

# 対象のGDBが格納されているフォルダ
folder = 'C:\ArcPySample'
gdb = "ArcPyTest.gdb"
data = join(folder,gdb)

# ファイルにヘッダを出力
outputstring = u"GDB,レイヤ名,エイリアス,ジオメトリタイプ,サブタイプ,アイテム数\n"
f.write(outputstring.encode("SHIFT-JIS"))

arcpy.env.workspace = data

# フィーチャクラスから属性を取得
for fc in arcpy.ListFeatureClasses():

    fList = []

    dataSource = join(data,fc)

    #フィーチャクラスの名称を取得
    name = arcpy.Describe(dataSource).Name
    #フィーチャクラスのエイリアスを取得
    aliasName = arcpy.Describe(dataSource).AliasName
    #フィーチャの件数取得
    cnt = arcpy.GetCount_management(dataSource)
    #フィーチャタイプ取得
    fType = arcpy.Describe(dataSource).FeatureType
    #ジオメトリタイプ取得
    gType = arcpy.Describe(dataSource).shapeType
    #アノテーションの場合（アノテーションのジオメトリタイプはポリゴンになってしまうのでこの処理を追加）
    if fType == "Annotation":
        gType = "Annotation"

    fList.append(gdb)
    fList.append(name)
    fList.append(aliasName)
    fList.append(gType)
    fList.append(int(cnt.getOutput(0)))

    getInfo(dataSource,fList)

f.close()

結果を確認します。

各フィーチャクラスもしくはフィーチャクラスのサブタイプごとのアイテム数が出力されましたね。GDBが複数存在する場合はちょっとプログラムを変更する必要がありますが、このような感じで既存のデータのアイテム数を知ることが出来ます。アイテム数だけではなく、知りたい項目に合わせてプログラムを修正することもできますね。

簡単ではありますが、本日は以上です！

【PythonとSQL】Zmap-TOWNⅡをファイルジオデータベースに変換して住所データ(地番戸番)を抜き出す方法 ～やっぱり便利なPandas～

さて、本日はZmap-TOWNⅡから住所データ(地番戸番)を抜き出してみようと思います。大変便利なZmap-TOWNⅡですが、これの住所一覧が欲しいなぁと思ってもゼンリンさんはそんなデータは用意してくれていませんね。なら自分で作ってしまいましょう！

Zmap-TOWNⅡをファイルジオデータベースに変換

まず、Zmap-TOWNⅡは独自の形式なのでファイルジオデータベースに変換します。

変換ツールの使用

以下のようにESRIから変換ツールをダウンロードします。
Zmap-TOWNII、Zmap-AREAII 変換ツールが更新！洗練された地図表現に : ESRIジャパン ArcGISブログ

ArcMapにこのようなツールを追加することができました。

これで変換ですね。

無事変換を終えたのですが、なんとレイヤ構成がZmap-TOWNⅡと違ってるじゃん！どうなってるのさ？！こんな感じなのですが、明らかにレイヤ構成がおかしい。

ということでESRIに問い合わせた結果、変換したファイルジオデータベースは通常のレイヤ追加ではなく、変換ツールからレイヤを追加しなくてはいけないとのこと。

実行した結果、そうそうこれ。これですよ！ESRIさんありがとうございますm(_ _)m

データの内容

今回一番欲しい建物のデータ定義はこんな感じになっています。

中身はこんな感じですね。大字や字丁目などがコードで入っているので、別テーブルと結合して名称を取得する必要がありますね。

サンプルコード

各テーブルやフィーチャクラスをPandasのデータフレームに格納してCSVに吐き出します。そのCSVをDB（SQLSERVER）に突っ込み、住所データ(地番戸番)を取得したいと思います。本当はCSVに出力ではなく、そのままDBに格納したかったのですが、どうしてもうまくいかなかったです・・・なので、それは今後の課題として、今回はCSV出力を行いたいと思います。

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

#カラム名と属性を返す関数
def GetAttribute(fc):
    field_list = []
    for field in arcpy.ListFields(fc):
        field_list.append(field.name)

    cur = arcpy.SearchCursor(fc)
    vallist2 = []
    for row in cur:
        vallist = []
        for field in field_list:
            vallist.append(row.getValue(field))
        vallist2.append(vallist)
    arr = np.array(vallist2)

    del row,cur
    return arr,field_list

#GDBに変換したZmapを指定
arcpy.env.workspace = ur"D:\python\zenrin\ZmapTown.gdb"

table = arcpy.ListTables()
table.sort()
jCode = ""

#テーブルから属性を取得
for tb in table:
    if tb ==  "ZmapTOWN_ZMAP_PST":
        data,olist = GetAttribute(tb)
        pst = pd.DataFrame(data)
        pst.columns = olist
        #一行だけでいいのでフィルターかける(もじかしたらここは市町村によっては変える必要ありかも)
        pstOneRow = pd.DataFrame(pst[pst.OID == "1"])
        jCode = pstOneRow["JCODE"].values[0]

    if tb ==  "ZmapTOWN_ZMAP_TATEMONO_NAME":
        data,olist = GetAttribute(tb)
        tatemono = pd.DataFrame(data)
        tatemono.columns = olist

        #重複データとCHIBAN、NAMEが空のデータは削除
        tatemonoValue = pd.DataFrame(tatemono[tatemono.CHIBAN <> ""])
        tatemonoValue = pd.DataFrame(tatemonoValue[tatemonoValue.NAME <> ""])

        #ここはキャストしない
        tatemonoValue = pd.DataFrame(tatemonoValue[tatemonoValue.JCODE == jCode])
        tatemonoValue = pd.DataFrame(tatemonoValue.drop_duplicates(['ACODE','CCODE','GCODE','CHIBAN','TPOLYCD']))

#GDBに変換したZmapのZmapTOWNフィーチャデータセットを指定
arcpy.env.workspace = ur"D:\python\zenrin\ZmapTown.gdb\ZmapTOWN"

# フィーチャクラスから属性を取得
for fc in arcpy.ListFeatureClasses():
    if fc ==  "ZmapTOWN_ZMAP_OOAZA":
        data,olist = GetAttribute(fc)
        ooaza = pd.DataFrame(data)
        ooaza.columns = olist

        #重複データとNAMEが空のデータは削除
        ooazaValue = pd.DataFrame(ooaza[ooaza.NAME.notnull()])
        ooazaValue = pd.DataFrame(ooazaValue[ooazaValue.JCODE == int(jCode)])
        ooazaValue = pd.DataFrame(ooazaValue.drop_duplicates(['NAME']))

    elif fc ==  "ZmapTOWN_ZMAP_AZA":
        data,olist = GetAttribute(fc)
        aza = pd.DataFrame(data)
        aza.columns = olist

        #重複データとNAMEが空のデータは削除
        azaValue = pd.DataFrame(aza[aza.NAME.notnull()])
        azaValue = pd.DataFrame(azaValue[azaValue.JCODE == int(jCode)])
        azaValue = pd.DataFrame(azaValue.drop_duplicates(['ACODE','CCODE','NAME']))

    elif fc ==  "ZmapTOWN_ZMAP_GAIKU":
        data,olist = GetAttribute(fc)
        gaiku = pd.DataFrame(data)
        gaiku.columns = olist

        #重複データとNAMEが空のデータは削除
        gaikuValue = pd.DataFrame(gaiku[gaiku.NAME.notnull()])
        gaikuValue = pd.DataFrame(gaikuValue[gaikuValue.JCODE == int(jCode)])
        gaikuValue = pd.DataFrame(gaikuValue.drop_duplicates(['ACODE','CCODE','GCODE','NAME']))

#データフレームをCSV出力
gaikuValue.to_csv("D:\python\zenrin\gaiku.csv",encoding = "utf-8")
azaValue.to_csv(r"D:\python\zenrin\aza.csv",encoding = "utf-8")
ooazaValue.to_csv("D:\python\zenrin\ooaza.csv",encoding = "utf-8")
tatemonoValue.to_csv(r"D:\python\zenrin\tatemono.csv",encoding = "utf-8")

出力したCSVはBULK INSERTなりなんなりでDBに突っ込んでください。←雑ですみません
テーブル定義はZmap-TOWNⅡの各テーブル、フィーチャクラスと同じにしました。

■大字

■字丁目

■街区

■建物

どのテーブルもコードばかりでわけがわかりませんね。定義書を見た感じ以下を考慮したSQLを作ればいいのではないかと思います。
・JCODE（拡張市町村コード）
・ACODE（大字コード）
・CCODE（字丁目コード）
・GCODE（街区コード）

select d.NAME   大字,
       c.NAME   字,
       b.NAME   街区,
       a.CHIBAN 地番戸番,
       a.NAME   名称
from tatemono a left outer join gaiku_zenrin b 
                    on (a.ACODE = b.ACODE and a.CCODE = b.CCODE and a.GCODE = b.GCODE)
                left outer join aza c 
                    on (a.ACODE = c.ACODE  and a.CCODE = c.CCODE)
                left outer join ooaza d 
                    on a.ACODE = d.ACODE
order by d.NAME,
         c.NAME,
         b.NAME,
         a.CHIBAN,
         a.NAME

結果を見るとこんな感じですね。けっこうあっさりできました。あとはこれをテーブルに突っ込んだり、EXCELに出力したりして利用することができますね！

全然記事の趣旨からはそれるのですが、今回の変換したZmap-TOWNⅡやデータコレクションなどを見ていると、明らかに動きが早いと感じています←やはりこれがESRIさんの実力なのか？！
フィーチャの数は結構多いのに、動きがかなり軽くて驚きですね。今度はどうしてこんなに動きを早くすることができるのか、その辺の分析を行ってみたいと思います。

以上、本日は「【PythonとSQL】Zmap-TOWNⅡをファイルジオデータベースに変換して住所データ(地番戸番)を抜き出す方法　～やっぱり便利なPandas～」でした。

ArcGISコミュニティの発足 ～Q&Aや情報交換用のFacebookグループ～

このたびArcGISコミュニティを発足することになりました。

ArcGISは便利なエンジンですが、世に出回る情報量の少なさが玉に瑕ですよね（あったとしても英語ばかりといったことも多々）。ですので、困った時に質問や情報収集が出来るコミュニティの存在が必要だと感じていました。このFacebookグループが日本におけるArcGISコミュニティの先駆けになれればと思っています。

■対象者

ArcGISに興味のある方全員

■対象技術

・ArcMap

・ArcObjects

・ArcPy

 ・その他

ArcGISに関わるすべての技術を対象とします。またトピック的なものも含めて情報交換できればと思っています。

■今後の展望

・定期交流会の開催（業種を問わず情報交換できればと思っています。学生さんも大歓迎です）

・ArcGISハッカソンの開催

ブログ主はまだまだArcGISを勉強中の身ですので、もっとArcGISを覚えたいという方と一緒に技術を高められれば幸いです。お気軽にご参加ください。

https://www.facebook.com/groups/1669235376647851/

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

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

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

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

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

このページをいったん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')

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

データ準備

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

ArcMapで読み込みます。

ポイントのフィーチャークラスを作成します。

フィーチャークラスの属性を以下の地理座標系にします。

サンプルコード②

「サンプルコード①」で出力した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()

ポイントがプロットされました。位置もあっているっぽいですね。

ポイントに属性も入っていますね。

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

サンプルコード③

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

# -*- 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ではできない模様。うー貧弱すぎる・・・

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

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

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

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