GIS奮闘記

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

ArcPyで壊れたフィーチャクラスを抽出する方法

さて、本日は「ArcPyで壊れたフィーチャクラスを抽出する方法」です。フィーチャクラスが壊れるって何?!と思う方もいらっしゃるかと思いますが、フィーチャクラスは壊れるんです。本当に貧弱です。

確認方法

属性テーブルのヘッダカラムをダブルクリックしてください。

例えば、こういったフィーチャクラスがある場合
f:id:sanvarie:20160401155648p:plain

フィーチャの件数は861件です。フィーチャクラスが壊れている状態でヘッダカラムをダブルクリックすると
f:id:sanvarie:20160401155850p:plain

フィーチャの件数が737件に減ってしまいました。

もうこうなってしまったら最後、新しくフィーチャクラスを作成するしかありません。また、確認方法ですが、上記の方法でも問題ないのですが、フィーチャクラスが数多くある場合、あまりにも面倒です。なので、Pythonでの抽出方法を考えてみました。

確認方法

arcpy.da.FeatureClassToNumPyArray を使用します。これはFeatureClassをNumPyArrayに変換してくれる便利な関数です。

①この関数の第二引数にOBJECTIDのみを渡す
→正常なフィーチャ数を取得(861件)

②この関数の第二引数にOBJECTID以外も渡す
→異常なフィーチャ数を取得(737件)

この仕組みを利用すれば、壊れたフィーチャクラスを抽出できますね。

環境

Windows7 64bit
ArcGIS10.2.0
Python2.7.3

サンプルソース

import arcpy

arcpy.env.workspace = "" #GDBが存在する一個上のフォルダを指定
broken_list = []

for w in arcpy.ListWorkspaces(workspace_type="FileGDB"):

    arcpy.env.workspace = w

    for f in arcpy.ListFeatureClasses():

        array1 = arcpy.da.FeatureClassToNumPyArray(f,["OID@"],null_value=-9999)
        array2 = arcpy.da.FeatureClassToNumPyArray(f,["OID@","SHAPE@XY"],null_value=-9999)

        if len(array1) != len(array2):
            #GDB、フィーチャクラス、正常な件数、異常な件数を格納
            broken_list.append([w,f,len(array1),len(array2)])

print broken_list

結果はこのようになります。壊れたフィーチャクラスの結果のみ抽出することができました。
f:id:sanvarie:20160401161336p:plain

そもそもこんなしょっちゅうフィーチャクラスが壊れるってどういうことですかESRIさん!と言いたいところですが、言ってもどうしようもありません。とりあえず、運悪く壊れてしまった場合はこのように対象フィーチャクラスを抽出すればいいですね。

簡単ですが、本日は以上です。

【PythonとGIS】ArcPyで頂点の間引き処理を進化させてみた

さて、本日は頂点の間引き処理を行ってみたいと思います。ArcMapの頂点の間引きというジオプロセシングツールを使えばそれを実現することが可能です。頂点の間引き (Generalize)—ヘルプ | ArcGIS for Desktop
ただ、以下のような問題もあります。

  • Standardライセンスからのみ使用可能で、Basicまでしか持っていないという場合は使用不可

今回はBasicライセンスしか持っていないという場合でもそれに対応できるPythonスクリプトを作成してみました。

ソースはGistにも公開しています。https://gist.github.com/sanvarie/51bcdd9e1668081bde75

環境

windows7 64bit
ArcGIS10.2.0
Python2.7.3

事前準備

  1. GDBを作成する(今回は「Generalize.gdb」を作成しました)
  2. フィーチャクラスを作成する(今回は「Line」というラインのフィーチャクラスを作成しました)

処理概要

  1. フィーチャクラスを作成する(間引いたラインを作成するフィーチャクラス)
  2. 始点終点は間引き処理対象外とする
  3. 頂点間の角度と距離を計測して間引き対象/対象外を選別する
  4. 頂点間の角度の差異が1度以下、かつ、頂点間の長さが1m以上を間引き対象とする

対象データ

こういったデータを作成しました。
f:id:sanvarie:20160325092748p:plain

各頂点はこのようになっています。
f:id:sanvarie:20160325093200p:plain

サンプルソース

■generalize.py

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

class Generalize(object):

    def __init__(self,workspace,in_feature,out_feature):
        arcpy.env.workspace = workspace
        self.workspace = workspace
        self.in_feature = in_feature
        self.out_feature = out_feature
        self.generalize_list = list()
        self.check_feature_class()

    def check_feature_class(self):
        """out_featureに指定したフィーチャクラスが無い場合、作成します"""

        exist_out_feature = arcpy.ListFeatureClasses(self.out_feature)

        if len(exist_out_feature) == 0:
            try:
                print u"フィーチャクラスを作成します"
                arcpy.CreateFeatureclass_management(self.workspace, self.out_feature, "POLYLINE")
            except Exception:
                u"フィーチャクラスの作成に失敗しました"

    def set_generalize_list(self):
        """角度と距離から削除対象の頂点か判定してリストに格納する

        self.generalize_list --ラインの始点終点などの情報
                         [0]:OBJECTID
                         [1]:X
                         [2]:Y
                         [3]:始点X
                         [4]:終点X
                         [5]:削除対象区分 0:削除対象,1:始点で削除非対称,2:終点で削除非対称
                         [6]:ポイントジオメトリ
                         [7]:角度
                         [8]:距離

        """

        for i,pt in enumerate(self.generalize_list):

            #削除対象のポイントの場合
            if pt[5] == 0:

                #角度計算用に対象の頂点とその前後の頂点の角度を取得する
                angle_first = self.generalize_list[i-1][7]
                angle_next = pt[7]
                angle_last = self.generalize_list[i+1][7]

                #対象の頂点とその前後のポイントの角度の差が1以下、かつ、対象の頂点間の長さが1以上の場合
                if math.fabs(math.fabs(angle_first) - math.fabs(angle_next)) <= 1 \
                and math.fabs(math.fabs(angle_next) - math.fabs(angle_last)) <= 1 and pt[8] >= 1:

                    self.generalize_list[i].append(1)

                else:

                    self.generalize_list[i].append(0)
            else:

                self.generalize_list[i].append(0)

    def calc_angle_length(self):
        """各頂点の角度や長さなどを計算する関数し、計算結果をリストに格納する関数

        self.generalize_list --ラインの始点終点などの情報
                         [0]:OBJECTID
                         [1]:X
                         [2]:Y
                         [3]:始点X
                         [4]:終点X
                         [5]:削除対象区分 0:削除対象,1:始点で削除非対称,2:終点で削除非対称
                         [6]:ポイントジオメトリ
                         [7]:角度
                         [8]:距離
                         [9]:削除対象区分2  0:削除対象外,1:削除 ←[5]はこれを抽出する前段階のもの
        """

        x2,y2,pg = float(),float(),object()
        generalize_list2 = list()

        for i,pt in enumerate(self.generalize_list):

            #始点、終点でもない場合
            if pt[5] == 0:

                #角度を計算
                result = math.atan2(pt[2]-y2,pt[1]-x2)

                #角度をセット
                self.generalize_list[i].append(math.degrees(result))

                #距離をセット
                self.generalize_list[i].append(pt[6].distanceTo(pg))

            #始点の場合
            elif pt[5] == 1:

                #角度を計算
                result = math.atan2(pt[2]-self.generalize_list[i+1][2],pt[1]-self.generalize_list[i+1][1])

                #角度をセット
                if math.degrees(result) < 0:
                    self.generalize_list[i].append(math.degrees(result) + 180)
                else:
                    self.generalize_list[i].append(math.degrees(result) - 180)

                #距離は0でOK
                self.generalize_list[i].append(0)

            #終点の場合
            elif pt[5] == 2:
                #角度を計算
                result = math.atan2(pt[2]-generalize_list2[i-1][2],pt[1]-generalize_list2[i-1][1])

                #角度をセット
                if math.degrees(result) < 0:
                    self.generalize_list[i].append(math.degrees(result))
                else:
                    self.generalize_list[i].append(math.degrees(result))

                #距離は0でOK
                self.generalize_list[i].append(0)

            generalize_list2.append(self.generalize_list[i])
            x2,y2,pg = pt[1],pt[2],pt[6]

        #角度と距離から削除対象の頂点か判定
        self.set_generalize_list()

    def get_point_coord(self):
        """ラインの頂点ごとの情報をリストに格納する"""

        lines = [row for row in arcpy.da.SearchCursor(self.in_feature, ["OID@","SHAPE@"])]

        #頂点ごとにOID、X、Y、始点、終点の情報を格納する
        for line in lines:
            for parts in line[1]:
                for part in parts:
                    self.generalize_list.append([line[0],part.X,part.Y,line[1].firstPoint.X,line[1].lastPoint.X])

    def calc_first_last(self):
        """各頂点が始点なのか終点なのかそれ以外なのか(削除対象)を判定する関数 """

        for i,pt in enumerate(self.generalize_list):

            #対象のポイントのXが始点の場合
            if pt[1] == pt[3]:

                #始点で削除非対象
                self.generalize_list[i].append(1)

            elif pt[1] == pt[4]:

                #終点で削除非対象
                self.generalize_list[i].append(2)

            else:

                #削除対象
                self.generalize_list[i].append(0)

    def calc_point_distance(self):
        """各種関数を呼び出し、self.generalize_listに要素を追加していく"""

        #ラインの頂点ごとの情報を取得
        self.get_point_coord()

        #各頂点の始点終点それ以外を判定
        self.calc_first_last()

        #頂点のジオメトリ情報を取得
        self.create_point_geometry()

        #頂点の角度と長さの計算
        self.calc_angle_length()

        #削除対象の頂点を間引いてラインを作図する
        self.generalize_line()

    def generalize_line(self):
        """頂点を間引いてラインを作図する関数

        self.generalize_list --ラインの始点終点などの情報
                         [0]:OBJECTID
                         [1]:X
                         [2]:Y
                         [3]:始点X
                         [4]:終点X
                         [5]:削除対象区分 0:削除対象,1:始点で削除非対称,2:終点で削除非対称
                         [6]:ポイントジオメトリ
                         [7]:角度
                         [8]:距離
                         [9]:削除対象区分2  0:削除対象外,1:削除 ←[5]はこれを抽出する前段階のもの
        """

        oid = int()
        array = arcpy.Array()

        with arcpy.da.InsertCursor(self.out_feature, ["OID@", "SHAPE@"]) as cur:

            for i,pt in enumerate(self.generalize_list):

                #OBJECTIDが次に行ったタイミングでラインを作図する。
                if pt[0] != oid and i != 0:
                    cur.insertRow((oid,arcpy.Polyline(array)))
                    array = arcpy.Array()
                    array.add(arcpy.Point(pt[1], pt[2]))
                else:
                    #削除対象の頂点ではない場合のみ頂点の情報を格納していく
                    if pt[9] != 1:
                        array.add(arcpy.Point(pt[1], pt[2]))

                oid = pt[0]

            cur.insertRow((oid,arcpy.Polyline(array)))

    def create_point_geometry(self):
        """ポイントジオメトリを作成してリストに格納"""
        point = arcpy.Point()

        for i,pt in enumerate(self.generalize_list):
            point.X = pt[1]
            point.Y = pt[2]
            self.generalize_list[i].append(arcpy.PointGeometry(point))

if __name__ == '__main__':

    #対象のgdb
    workspace = r"C:\ArcPySample\Generalize.gdb"

    #間引き元のフィーチャクラス名
    in_feature = "Line"

    #間引き先のフィーチャクラス名
    out_feature = "LineGeneralize"

    #インスタンス化
    generalize = Generalize(workspace,in_feature,out_feature)

    #頂点を間引く
    generalize.calc_point_distance()

結果を確認すると問題なく、頂点が間引きされていることがわかりました。
f:id:sanvarie:20160325100013p:plain

課題

  1. データ量が多い場合は少し時間がかかってしまう
  2. 属性が存在する場合は、間引いた先のフィーチャクラスにはその属性は引き継がれない(generalize_lineメソッドを修正すればいける)

課題は少しありますが、角度や距離をうまく使用すれば、ArcMapの機能よりも柔軟な使用ができそうですね。簡単ではありますが、本日は以上です。

【PythonとGIS】ArcPyでラインの端点の角度を取得する

さて、本日はラインの端点の角度を取得しようと思います。ラインの端点にはシンボルを置くことが多いので、端点の角度を取得して、そこにのっかるポイントに角度を与える場合などに使えるかと思います。

環境

Windows7 64bit
ArcGIS10.2
Python2.7.3

事前準備

  1. 対象GDBの作成(今回はGeneralize.gdbを作成しました)
  2. フィーチャクラスの作成(今回はLineを作成しました)

プログラムで取得できる値

対象のフィーチャクラスにおいて、端点の角度を含む以下情報をListで返します。端点以外の頂点の情報は格納しません。

[0]:OBJECTID
[1]:X
[2]:Y
[3]:頂点区分 1:始点 2:終点
[4]:角度

対象データ

以下のようなサンプルデータを作成しました。
f:id:sanvarie:20160322151530p:plain

左から順に作図しました。
f:id:sanvarie:20160322151650p:plain

サンプル

■calcangle.py

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

class CalcAngle(object):

    def __init__(self,workspace,in_feature):
        arcpy.env.workspace = workspace
        self.in_feature = in_feature
        self.generalize_list = list()

    def calc_angle(self):
        """各頂点の角度を計算する関数

        self.generalize_list --ラインの始点終点などの情報
                         [0]:OBJECTID
                         [1]:X
                         [2]:Y
                         [3]:頂点区分 1:始点,2:終点
                         [4]:角度
        """

        for i,pt in enumerate(self.generalize_list):

            #始点の場合
            if pt[3] == 1:

                #角度を計算
                result = math.atan2(pt[2]-self.generalize_list[i+1][2],pt[1]-self.generalize_list[i+1][1])

                #角度をセット
                if math.degrees(result) < 0:

                    self.generalize_list[i].append(abs(math.degrees(result)))

                else:

                    self.generalize_list[i].append(abs(math.degrees(result) - 360))

            #終点の場合
            elif pt[3] == 2:

                #角度を計算
                result = math.atan2(pt[2]-self.generalize_list[i-1][2],pt[1]-self.generalize_list[i-1][1])

                #角度をセット
                if math.degrees(result) < 0:

                    self.generalize_list[i].append(abs(math.degrees(result)))

                else:

                    self.generalize_list[i].append(abs(math.degrees(result) -360))

            else:
                self.generalize_list[i].append(0)

    def get_point_coord(self):
        """ラインの頂点ごとの情報をリストに格納する"""

        lines = [row for row in arcpy.da.SearchCursor(self.in_feature, ["OID@","SHAPE@"])]

        #頂点ごとにOID、X、Y、始点、終点の情報を格納する
        for line in lines:
            for parts in line[1]:
                for part in parts:
                    #始点の場合
                    if part.X == line[1].firstPoint.X and part.Y == line[1].firstPoint.Y:
                        self.generalize_list.append([line[0],part.X,part.Y,1])
                    #終点の場合
                    elif part.X == line[1].lastPoint.X and part.Y == line[1].lastPoint.Y:
                        self.generalize_list.append([line[0],part.X,part.Y,2])

    def calc_point_distance(self):
        """各種関数を呼び出し、self.generalize_listに要素を追加していく"""

        #始点終点の情報を取得
        self.get_point_coord()

        #頂点の角度の計算
        self.calc_angle()

        return self.generalize_list

if __name__ == '__main__':

    #対象のgdb
    workspace = r"C:\ArcPySample\Generalize.gdb"

    #対象のフィーチャクラス名
    in_feature = "Line"

    #インスタンス化
    generalize = CalcAngle(workspace,in_feature)

    #角度計算
    angle_list = generalize.calc_point_distance()

    print angle_list

実行結果はこのようになりました。
f:id:sanvarie:20160322152039p:plain

OBJECTIDが1のデータで確認すると始点(以下画像の緑丸)の角度=156.1363677663898、終点(以下画像の赤丸)の角度=336.1363677663898 で問題なさそうですね。
f:id:sanvarie:20160322152247p:plain

これを使用すれば端点の角度を利用してそこに接するポイントに角度を付与することができるかと思います。また、これを応用して、各頂点ごとの角度などを取得してもおもしろそうですね。

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

【PythonとGIS】ArcPyでライン接合を行ってみた ~Advancedライセンス?!だが断る!~

さて、本日はArcPyでライン接合を行ってみようと思います。背景などでけっこうラインがぶちぶち切れてしまっていることが多いので、この機能があればラインのフィーチャ数を減らすことができますね。Advancedライセンスがあれば、ArcMapで実行できるのですが(うまくいくかは知らない)、いかんせんAdvancedはお高いですよね(けっこういい車が買えるお値段です)。そんなものを購入できるはずもなく、「使えないなら作ろう!」と思い立った次第であります。ArcMapでのライン接合についてはこのサイトに記載してあります。ラインの接合 (Unsplit Line)—データ管理ツールボックス | ArcGIS for Desktop

今回もソースをGistに公開しました。https://gist.github.com/sanvarie/410ea0c1bab366f85a68

環境

Windows7 64bit
ArcGIS10.2
Python2.7.3

事前準備

  1. ファイルジオデータベースの作成(今回はArcPyAdvanced.gdbという名前にしました)
  2. 接合対象のフィーチャクラスの作成(今回はLineSplitという名前にしました)
  3. 接合先のフィーチャクラスの作成(今回はLineUnsplitという名前にしました)←接合先のフィーチャクラスが存在しなかったら作成するようにしたので、作らなくても大丈夫です。
  4. Pandasのインストール

接合対象/対象外データ

以下参照
f:id:sanvarie:20160319170244p:plain

サンプルソース

ライン接合処理のサンプルです。以下に処理の概要を記載します。

  1. 接合対象ラインを抽出(つながりのあるラインを抽出)
  2. 接合対象ラインの中でループしているラインを判定
  3. 接合対象ラインの中で端のラインとそのラインの端点(他のラインと接していない端点)を抽出(ループしているラインは除く)
  4. 接合対象ラインの端点から作図順を決める
  5. 接合対象ラインを作図
  6. 接合対象外ラインを作図

■unsplitline.py

# -*- coding: utf-8 -*-
import arcpy
from collections import Counter
import os
import traceback
import pandas as pd

class UnsplitLine():

    column_unsplit_length = "UNSPLIT_LENGTH" #作成するカラム名

    def __init__(self,workspace,in_feature,out_feature):
        arcpy.env.workspace = workspace
        self.workspace = workspace
        self.in_feature = in_feature
        self.out_feature = out_feature
        self.line_list = []
        self.match_point_list = []
        self.finish_line_list = []
        self.divide_line_list = []
        self.list_draw_point_sort = []
        self.finish_oid = []
        self.count = UnsplitLine.counter(self)
        self.check_feature_class()
        self.check_columns()

    def check_feature_class(self):
        """out_featureに指定したフィーチャクラスが無い場合、作成します"""

        exist_out_feature = arcpy.ListFeatureClasses(self.out_feature)

        if len(exist_out_feature) == 0:
            print u"フィーチャクラスを作成します"
            arcpy.CreateFeatureclass_management(self.workspace, self.out_feature, "POLYLINE")

    def check_columns(self):
        """out_featureに指定したフィーチャクラスに「UNSPLIT_LENGTH」カラムが存在しない場合、カラムを作成します"""
        field_list = []
        fields = arcpy.ListFields(self.out_feature)
        for field in fields:
            if field.type != "Geometry":
                field_list.append(field.name)

        if not UnsplitLine.column_unsplit_length in field_list:
            print u"カラムを作成します"
            arcpy.AddField_management(self.out_feature, UnsplitLine.column_unsplit_length, "Double")

    def export_upsplit_lines(self):
        """ライン接合後、Shape_lengthやアイテム数をテキスト出力します"""
        try:
            #フィーチャ数のLengthの集計
            moto_line_length,saki_line_length = self.count_length()

            #フィーチャ数のカウント
            result_in_feature,result_out_feature = self.count_feature()

            file_dir = self.workspace[:self.workspace.rfind("\\")]

            with open(os.path.join(file_dir,"export_upsplit_lines.txt"), "w") as f:

                message = u"接合元フィーチャの総延長=" + str(moto_line_length)
                f.write(message.encode("Shift-JIS"))
                f.write("\n")

                message = u"接合先フィーチャの総延長=" + str(saki_line_length)
                f.write(message.encode("Shift-JIS"))
                f.write("\n")

                message = u"接合元フィーチャ数=" + str(result_in_feature)
                f.write(message.encode("Shift-JIS"))
                f.write("\n")

                message = u"接合先フィーチャ数=" + str(result_out_feature)
                f.write(message.encode("Shift-JIS"))
                f.write("\n")

                message = u"接合したラインのOID"
                f.write(message.encode("Shift-JIS"))
                f.write("\n")
                for line in self.line_list:
                    for l in line:
                        f.write(str(l))
                        f.write(" ")
                    f.write("\n")

        except IOError:
            print u"ファイルの出力に失敗しました。"
        except TypeError:
            print u"ファイルに出力する内容が不正です。"

    def count_feature(self):
        """接合元と接合先のフィーチャクラスのフィーチャ数をカウントします"""
        result_in_feature = arcpy.GetCount_management(self.in_feature)
        result_out_feature = arcpy.GetCount_management(self.out_feature)

        return result_in_feature,result_out_feature

    def count_length(self):
        """接合元と接合先のフィーチャクラスのShape_Lengthを集計します"""
        moto_line_length = float()
        saki_line_length = float()

        for line in [row for row in arcpy.da.SearchCursor(self.in_feature, ["OID@","SHAPE@"])]:
            moto_line_length += line[1].length

        for line in [row for row in arcpy.da.SearchCursor(self.out_feature, ["OID@","SHAPE@"])]:
            saki_line_length += line[1].length

        return moto_line_length,saki_line_length

    def set_join_oid(self,oid_list):
        """oid_listに格納されているOIDを分解します。[1,2,3]の場合、1,2,3に分解します。"""
        oid_list = set(oid_list)
        oid_list = list(oid_list)
        oid = ",".join([str(x) for x in oid_list])
        return oid

    def copy_feature(self):
        """接合対象ではないラインの作図を行います"""

        try:
            oid = self.set_join_oid(self.finish_line_list)

            with  arcpy.da.InsertCursor(self.out_feature, ["OID@","SHAPE@",UnsplitLine.column_unsplit_length]) as cur:
                for line in arcpy.da.SearchCursor(self.in_feature, ["OID@","SHAPE@"],"OBJECTID not in (%s)" % oid) :
                    length = self.get_shape_length(line[0])
                    for part in line[1]:
                        array = arcpy.Array()
                        for pnt in part:
                            if pnt:
                                array.add(arcpy.Point(pnt.X, pnt.Y))

                    cur.insertRow((line[0],arcpy.Polyline(array),length))

        except TypeError:
            print traceback.format_exc()

    def get_lines(self,oid=0):
        """接合元のラインを取得する関数です

        returns:
        self.moto_lines -- 接合元のラインの全フィーチャ
        line            -- ラインオブジェクト
        """

        if oid == 0:
            self.moto_lines = [row for row in arcpy.da.SearchCursor(self.in_feature, ["OID@","SHAPE@"])]
            return self.moto_lines
        else:
            line = [row for row in arcpy.da.SearchCursor(in_feature, ["OID@","SHAPE@"],"OBJECTID = %s" % oid)]
            line = [x for x in line for x in x]
            return line

    def search_line_combination(self,line2):
        """接合するラインの組み合わせを取得する関数です

            例:OBJECTID=1,2と4,5,6と10,11,12,13が接続対象の場合
               [[1,2],[4,5,6],[10,11,12,13]]というリストを作成します。

        arguments:
        line2 -- search_lineでのループから渡されるLineオブジェクト

        """

        point_list = []  #端点格納用リスト
        oid = line2[0]

        try:
            for line in self.moto_lines:
                if line[0] != oid:
                    #組み合わせ登録が完了したリストと端点が3か所重なっているライン(接合対象外)のリスト以外
                    if not line[0] in self.finish_line_list:
                        if (line[1].firstPoint.X == line2[1].lastPoint.X and line[1].firstPoint.Y == line2[1].lastPoint.Y) or \
                           (line[1].lastPoint.X == line2[1].lastPoint.X and line[1].lastPoint.Y == line2[1].lastPoint.Y) or \
                           (line[1].lastPoint.X == line2[1].firstPoint.X and line[1].lastPoint.Y == line2[1].firstPoint.Y) or \
                           (line[1].firstPoint.X == line2[1].firstPoint.X and line[1].firstPoint.Y == line2[1].firstPoint.Y):

                            point_list.append(line[0])
                            line_copy = line

            #隣り合っているラインが無い、もしくは、接合対象ラインの組み合わせリストができた場合のみこの分岐に入る
            if len(point_list) == 0:

                #接合対象ラインの組み合わせリストが存在する場合
                if len(self.match_point_list) != 0:
                    self.match_point_list = set(self.match_point_list) #上記処理だとOIDが重複してしまう可能性があるので、重複削除
                    self.line_list.append(list(self.match_point_list))
                    self.match_point_list = []

            #一番端のラインの場合はこの分岐に入る
            elif len(point_list) == 1:
                point_list.append(oid)
                self.finish_line_list.extend(point_list)
                self.match_point_list.extend(point_list)

                #再検索
                self.search_line_combination(line_copy)

            #両端にラインが存在する場合はこの分岐に入る(この分岐では組み合わせ登録は行わない)
            elif len(point_list) == 2:

                oid_list = []

                #端点が三つ以上重なっているか確認
                ret,point_list = self.check_line_combination(oid,point_list)

                #端点が三つ以上重なっている場合
                if ret == -1:
                    self.divide_line_list.extend(point_list)
                else:
                    oid_list.append(oid)
                    self.match_point_list.extend(oid_list)

                    #ループしている場合はこの分岐に入る
                    if len(self.match_point_list) > len(point_list):
                        self.finish_line_list.extend(self.match_point_list)

                    #重複していない要素を抽出
                    l = list(x for x in point_list if x not in self.match_point_list)

                    if len(l) == 1:
                        line = self.get_lines(l[0])
                        self.search_line_combination(line)
                    else:
                        self.search_line_combination(line_copy)

        except Exception:
            print traceback.format_exc()

    def check_line_combination(self,oid,point_list):
        """接合対象のラインにおいて、端点が三つ以上重なっているか確認する関数です

        arguments:
        oid -- search_lineでのループから渡されるOBJECTID
        point_list -- search_line_combinationで抽出した端点のリスト

        returns:
        1 or -1 -- 1:端点が三つ重なっていない -1:端点が三つ重なっている
        point_list -- point_list

        """

        point_list.append(oid)
        point_list = set(point_list)
        point_list = list(point_list)
        oid = ",".join([str(x) for x in point_list])

        devide_list = []

        #対象のOID(複数)で検索
        for row in arcpy.da.SearchCursor(in_feature, ["OID@", "SHAPE@"],"OBJECTID in (%s)" % oid):

            devide_list.append([row[0],row[1].firstPoint.X,row[1].firstPoint.Y,row[1].lastPoint.X,row[1].lastPoint.Y])

        #データフレームに格納
        df_devide = pd.DataFrame(devide_list)
        df_devide.columns = ["OID","FirstX","FirstY","LastX","LastY"]

        #始点と終点でグルーピング
        dic_first = df_devide.groupby("FirstX").groups
        dic_last = df_devide.groupby("LastX").groups
        list_first = [x for _,x in dic_first.iteritems()]
        list_last = [x for _,x in dic_last.iteritems()]

        #始点と終点で同じ座標が三つ以上ある場合は処理対象外とする(すべて始点、すべて終点の場合)
        for f,l in zip(list_first,list_last):
            if len(f) >= 3 or len(l) >= 3:
                return -1,point_list

        #始点と終点で同じ座標が三つ以上ある場合は処理対象外とする(始点と終点がそれぞれ3点に集中している場合)
        divide_line_list = []
        for key,row in df_devide.iterrows():
            for key2,row2 in df_devide.iterrows():
                if key !=  key2:
                    if row.FirstX == row2.LastX or row.FirstX == row2.FirstX:
                        divide_line_list.append(row.OID)
            if len(divide_line_list) >= 2:
                return -1,point_list

            divide_line_list = []

        return 1,point_list

    def search_line(self):
        """search_line_combinationでの結果を返す関数です

        returns:
        self.line_list -- 接合対象ラインのリスト(多次元)

        例:OBJECTID=1,2と4,5,6と10,11,12,13が接続対象の場合
        [[1,2],[4,5,6],[10,11,12,13]]というリストが返されます。

        """

        lines = self.get_lines()

        for line in lines:

            self.search_line_combination(line)

        return self.line_list

    def pop_list(self,pop_list):
        """リストから要素をポップする関数です

        arguments:
        pop_list -- ポップ対象のリスト

        """
        try:
            for p in pop_list:
                i = pop_list.index(p)
                self.moto_lines.pop(i)
        except ValueError:
            print traceback.format_exc()
            raise Exception

    def get_coord(self,oid):
        """接合対象ラインにおいて、各頂点の座標を取得する関数です

        arguments:
        oid -- OBJECTID → 接合対象のOBJECTIDが1,2,3 の場合、oid = 1,2,3 という形になります

        returns:
        coord_list -- 接合対象ラインの、各頂点の座標(多次元)
                     [0] = OBJECTID
                     [1] = ポイントごとに振ったID
                     [2] = X
                     [3] = Y
                     [4] = 始点X
                     [5] = 始点Y
                     [6] = 終点X
                     [7] = 終点Y

        loop_list -- 接合対象ラインがループしているのかどうかの判定用リスト(多次元)
                    [0] = OBJECTID
                    [1] = 始点X
                    [2] = 始点Y
                    [3] = 終点X
                    [4] = 終点Y

        """

        coord_list = []
        loop_list = []

        #対象のOID(複数)で検索
        for row in arcpy.da.SearchCursor(self.in_feature, ["OID@", "SHAPE@"],"OBJECTID in (%s)" % oid):

            loop_list.append([row[0],row[1].firstPoint.X,row[1].firstPoint.Y,row[1].lastPoint.X,row[1].lastPoint.Y]) #ループしたラインの判定用に使用する

            for part in row[1]:

                for i,pnt in enumerate(part):

                    if pnt:
                        #ポイントにIDを振る
                        pnt.ID = i

                        #ラインのOIDと座標を格納
                        coord_list.append([row[0],pnt.ID,pnt.X,pnt.Y,row[1].firstPoint.X,row[1].firstPoint.Y,row[1].lastPoint.X,row[1].lastPoint.Y])

        return coord_list,loop_list

    def check_geometry(self,df):
        """接合対象ラインにおいて、各頂点が始点なのか終点なのか判定する関数
           判定した結果を、データフレームに結合します。

        arguments:
        df -- Pandasのデータフレーム
              columns = ["OID","ID","X","Y","FirstX","FirstY","LastX","LastY"]


        returns:
        df_merge -- Pandasのデータフレーム
                    columns = ["OID","ID","X","Y","FirstX","FirstY","LastX","LastY","PointType"]
                    PointType 1:始点 2:終点 None:それ以外

        """

        try:
            list_point_type = []

            #各頂点が始点なのか終点なのかを確認
            for _ ,row in df.iterrows():
                for _ ,row2 in df.iterrows():
                    if row.OID == row2.OID and row.ID == row2.ID:
                        if row.X == row2.FirstX and row.Y == row2.FirstY:
                            #始点として登録
                            list_point_type.append([row2.OID,row2.ID,"1"])
                        elif row.X == row2.LastX and row.Y == row2.LastY:
                            #終点として登録
                            list_point_type.append([row2.OID,row2.ID,"2"])

            #リストをデータフレームに格納
            df_point_type = pd.DataFrame(list_point_type)
            df_point_type.columns = ["OID","ID","PointType"]
            df_point_type = pd.DataFrame(df_point_type.drop_duplicates(["OID","PointType"]))

            #データフレームを結合
            df_merge = pd.merge(df, df_point_type, on=["OID","ID"],how='left')

            return df_merge
        except Exception:
            print traceback.format_exc()

    def check_first_last(self,df):
        """接合対象ラインにおいて、端のラインと端が始点なのか終点なのかを判定する関数です

        arguments:
        df -- Pandasのデータフレーム。
              columns = ["OID","ID","X","Y","FirstX","FirstY","LastX","LastY","PointType"]
              PointType 1:始点 2:終点 None:それ以外

        returns:
        list_first_last -- 端のラインとPointTypeを格納したリスト
                           [0]:OBJECTID
                           [1]:頂点ごとに振ったID
                           [2]:PointType 1:始点 2:終点
        """

        list_first_last = []
        list_edge = []

        #ラインのつながりが二本の場合
        if len(df.groupby("OID").groups.keys()) == 2:
            for _ , row in df.iterrows():
                for _ , row2 in df.iterrows():
                    if row.OID !=  row2.OID:
                        if row.PointType == "1":
                            if (row.LastX == row2.FirstX and row.LastY == row2.FirstY) or \
                               (row.LastX == row2.LastX and row.LastY == row2.LastY):

                                list_first_last.append([row.OID,row.ID,"1"])

                        elif row.PointType == "2":

                            if (row.FirstX == row2.LastX and row.FirstY == row2.LastY) or \
                               (row.FirstX == row2.FirstX and row.FirstY == row2.FirstY):

                                list_first_last.append([row.OID,row.ID,"2"])

                        else:

                            continue

            return list_first_last

        else:
            for _ , row in df.iterrows():
                for _ , row2 in df.iterrows():
                    if row.OID != row2.OID:
                        if row.PointType == "1":
                            if (row.LastX == row2.FirstX and row.LastY == row2.FirstY) or \
                               (row.LastX == row2.LastX and row.LastY == row2.LastY):

                                list_edge.append([row.OID,row.ID,"1"])

                            if (row.FirstX == row2.FirstX and row.FirstY == row2.FirstY) or \
                               (row.FirstX == row2.LastX and row.FirstY == row2.LastY):

                                list_edge.append([row.OID,row.ID,"2"])

                        elif row.PointType == "2":
                            if (row.FirstX == row2.LastX and row.FirstY == row2.LastY) or \
                               (row.FirstX == row2.FirstX and row.FirstY == row2.FirstY):

                                list_edge.append([row.OID,row.ID,"2"])

                            if (row.LastX == row2.FirstX and row.LastY == row2.FirstY) or \
                               (row.LastX == row2.LastX and row.LastY == row2.LastY):

                                list_edge.append([row.OID,row.ID,"1"])

                #重複を削除([[1,2,3],[1,2,3]]があったら[[1,2,3]]にする)
                c = Counter(tuple(x) for x in list_edge)
                list_edge = [list(k) for k,v in c.items()]

                if len(list_edge) >= 2:
                    list_edge = []

                elif len(list_edge) == 1:
                    list_first_last.append([list_edge[0][0],list_edge[0][1],list_edge[0][2]])

                    return list_first_last

    def counter(self):
        """頂点の作図順用のカウンタ"""
        n = [0]
        def inc():
            n[0] += 1
            return n[0]
        return inc

    def search_dataframe(self,df,oid,line_id):
        """データフレームのレコードを抽出する用の関数"""
        df = pd.DataFrame(df[df.OID == oid])
        df = pd.DataFrame(df[df.ID == line_id])

        return df

    def draw_point_sort(self,df,oid,line_id,point_type):
        """接合対象ラインにおいて、各頂点の作図順を設定する関数です
      self.list_draw_point_sort に作図順を格納して後続処理で使用します

        arguments:
        df -- Pandasのデータフレーム。
              columns = ["OID","ID","X","Y","FirstX","FirstY","LastX","LastY","PointType"]
              PointType 1:始点 2:終点 None:それ以外
        oid -- OBJECTID(端のライン)
        line_id -- 各頂点に割り当てたID(端のライン)
        point_type -- 1:始点 2:終点 ←左記以外の値(None)は受け取らない想定(端のラインのどちらかの端点が入っているため)

        """

        #データフレームから対象のデータを抽出する(必ず一件になるはず)
        df_draw_point_sort = self.search_dataframe(df,oid,line_id)

        #対象のデータのX,Y,ポイントタイプ(始点or終点)
        first_x = df_draw_point_sort["FirstX"].values[0]
        last_x = df_draw_point_sort["LastX"].values[0]
        first_y = df_draw_point_sort["FirstY"].values[0]
        last_y = df_draw_point_sort["LastY"].values[0]

        if point_type == "":
            point_type = df_draw_point_sort["PointType"].values[0]

        #始点の場合
        if point_type == "1":
            #データフレームから対象のデータを抽出する(必ず複数件になるはず)
            df_asc = pd.DataFrame(df[df.OID == oid])

            #ポイントIDの昇順にソート(つまり、始点側から順にソートする)
            df_asc = df_asc.sort(['ID'], ascending=[True])

            #対象のラインの作図順を決定する
            for _ ,row in df_asc.iterrows():
                self.list_draw_point_sort.append([row.OID,row.ID,self.count()])

        #終点の場合
        elif point_type == "2":
            #データフレームから対象のデータを抽出する(必ず複数件になるはず)
            df_desc = pd.DataFrame(df[df.OID == oid])

            #ポイントIDの降順にソート(つまり、終点側から順にソートする)
            df_desc = df_desc.sort(['ID'], ascending=[False])

            #対象のラインの作図順を決定する
            for _ ,row in df_desc.iterrows():
                self.list_draw_point_sort.append([row.OID,row.ID,self.count()])

        #対象のラインに接続するラインで再度処理を行う
        for _ ,row in df.iterrows():
            if row.OID != oid:
                #対象のラインを始点から作図する場合は隣のラインは必ず対象のラインの終点に一致する
                if point_type == "1":

                    if row.X == last_x and row.Y == last_y:

                        if not row.OID in self.finish_oid:

                            self.finish_oid.append(oid)
                            self.draw_point_sort(df,row.OID,row.ID,"") #再検索(更に隣のラインを検索)

                #対象のラインを終点から作図する場合は隣のラインは必ず対象のラインの始点に一致する
                elif point_type == "2":
                    if row.X == first_x and row.Y == first_y:

                        if not row.OID in self.finish_oid:
                            self.finish_oid.append(oid)
                            self.draw_point_sort(df,row.OID,row.ID,"") #再検索(更に隣のラインを検索)

    def search_loop_line(self,df_loop):
        """接合対象ラインがループしているラインなのか判定する関数

        arguments:
        df_loop -- Pandasのデータフレーム。
                  columns = ["OID","FirstX","FirstY","LastX","LastY"]

        returns:
        1 or -1-- 1:ループしたラインではない -1:ループしたライン

        """

        #self.finish_line_listに格納されていないものはここで格納
        for _,row in df_loop.iterrows():
            if not row.OID in self.finish_line_list:
                self.finish_line_list.append(row.OID)

        #接合対象のラインが二本の場合
        if len(df_loop) == 2:

            for _,row in df_loop.iterrows():
                for _,row2 in df_loop.iterrows():
                    if row.OID !=  row2.OID:
                        if ((row.FirstX == row2.LastX and row.FirstY == row2.LastY) or \
                           (row.FirstX == row2.FirstX and row.FirstY == row2.FirstY)) and \
                           ((row.LastX == row2.FirstX and row.LastY == row2.FirstY) or \
                           (row.LastX == row2.LastX and row.LastY == row2.LastY)):

                            break

                #ループしたラインではない場合
                else:
                    return 1

            #ループしたラインの場合(breakしたらこっちに入る)
            return -1

        elif len(df_loop) >= 3:
            oid_list = []

            #ループしているラインを判定(接合対象ラインにおいて、全始点終点が他ラインの始点終点と一致した場合はループしている)
            for _,row in df_loop.iterrows():
                for _,row2 in df_loop.iterrows():
                    if row.OID !=  row2.OID:
                        if (row.FirstX == row2.FirstX and row.FirstY == row2.FirstY) or \
                           (row.FirstX == row2.LastX and row.FirstY == row2.LastY):

                            oid_list.append(1)

                        if (row.LastX == row2.FirstX and row.LastY == row2.FirstY) or \
                           (row.LastX == row2.LastX and row.LastY == row2.LastY):

                            oid_list.append(2)

                #ループしてないラインの場合はこの分岐に入る
                if len(oid_list) != 2:
                    return 1

                oid_list = []

            #ループしたラインの場合
            return -1

    def get_shape_length(self,oid):
        """接合対象ラインのShape_lengthの合計を取得する関数"""

        total_length = float()
        for line in arcpy.da.SearchCursor(self.in_feature, ["SHAPE@"],"OBJECTID in (%s)" % oid) :
             total_length += line[0].length

        return total_length

    def upsplit_line(self,cur,df,oid):
        """ラインを接合する関数

        arguments:
        cur -- カーソル
        df -- Pandasのデータフレーム
              columns = ["OID","ID","X","Y","FirstX","FirstY","LastX","LastY","PointSort"]
        oid -- OBJECTID(1,2,3が接合対象の場合oidには1,2,3 が格納されている)

        """

        try:
            array = arcpy.Array()

            length = self.get_shape_length(oid)

            for d in df.index:
                array.add(arcpy.Point(df.ix[d][2], df.ix[d][3]))
            cur.insertRow((arcpy.Polyline(array),length))
        except Exception:
            print traceback.format_exc()

    def create_line(self,lines):
        """各種関数を呼び出す処理。ここで、ラインの接合や接合対象ではないラインのコピー関数を呼び出す。

        arguments:
        lines -- search_lineで格納したリスト(多次元)

                 例:OBJECTID=1,2,3 また、4,5 というラインが接合対象の場合、
                    [[1,2,3],[4,5]] という要素が格納されている

        """

        try:
            if len(lines) == 0:
                print u"対象データがありません"
                raise Exception

            with  arcpy.da.InsertCursor(self.out_feature, ["SHAPE@",UnsplitLine.column_unsplit_length]) as cur:

                #接合対象ラインのリストでループ
                for line in lines:

                    #接合対象のラインのリストを分解(例 [1,2,3]→1,2,3)
                    oid = self.set_join_oid(line)

                    #接合対象ラインの座標をリスト化する
                    coord_list,loop_list = self.get_coord(oid)

                    df = pd.DataFrame(loop_list)
                    df.columns = ["OID","FirstX","FirstY","LastX","LastY"]

                    #対象のラインがループしているか判定
                    result = self.search_loop_line(df)

                    #ループしたラインがない場合
                    if result == 1:

                        #データフレーム作成
                        df_coord = pd.DataFrame(coord_list)
                        df_coord.columns = ["OID","ID","X","Y","FirstX","FirstY","LastX","LastY"]

                        #各頂点を始点、終点、それ以外に仕分け
                        df_checked = self.check_geometry(df_coord)

                        #各ラインのつながりにおいて、端のラインと端点を抽出する
                        list_checked = self.check_first_last(df_checked)

                        #上記で端のラインが判定できなかった場合、エラーとする
                        if not list_checked:
                            print u"端点が判断できないラインが存在しました" + " " + "OID=" + oid
                            raise Exception

                        #接合対象のラインにおいて、端のライン
                        oid_one,line_id,point_type = list_checked[0][0],list_checked[0][1],list_checked[0][2]

                        #隣りあっているラインを抽出して作図順を決めていく
                        self.draw_point_sort(df_checked,oid_one,line_id,point_type)

                        #リストをデータフレーム(各頂点の作図順を格納している)に格納
                        df_point_sort = pd.DataFrame(self.list_draw_point_sort)
                        df_point_sort.columns = ["OID","ID","PointSort"]

                        #データフレームを結合
                        df_checked = pd.merge(df_checked, df_point_sort, on=["OID","ID"],how="left")

                        #重複を削除
                        df_checked = pd.DataFrame(df_checked.drop_duplicates(["X","Y"]))

                        #各頂点の作図順でソート
                        df_checked = df_checked.sort(["PointSort"], ascending=[True])

                        #ラインを作図
                        self.upsplit_line(cur,df_checked,oid)

                        self.list_draw_point_sort = []

                    #ループしたラインがある場合
                    else:

                        #データフレーム作成
                        df_coord = pd.DataFrame(coord_list)
                        df_coord.columns = ["OID","ID","X","Y","FirstX","FirstY","LastX","LastY"]

                        #各頂点を始点、終点、それ以外に仕分け
                        df_checked = self.check_geometry(df_coord)

                        #隣りあっているラインを抽出して作図順を決めていく
                        self.draw_point_sort(df_checked,coord_list[0][0],coord_list[0][1],"")

                        #リストをデータフレーム(各頂点の作図順を格納している)に格納
                        df_point_sort = pd.DataFrame(self.list_draw_point_sort)
                        df_point_sort.columns = ["OID","ID","PointSort"]

                        #データフレームを結合
                        df_checked = pd.merge(df_checked, df_point_sort, on=["OID","ID"],how="left")

                        #各頂点の作図順でソート
                        df_checked = df_checked.sort(["PointSort"], ascending=[True])

                        #ラインを作図
                        self.upsplit_line(cur,df_checked,oid)

                        self.list_draw_point_sort = []

            #非接合対象のラインをコピーする
            self.copy_feature()

            #接合したラインの情報をテキストで出力
            self.export_upsplit_lines()

        except Exception:
            print traceback.format_exc()

if __name__ == '__main__':

    #対象のgdb
    workspace = r"C:\ArcPySample\ArcPyAdvanced.gdb"

    #接合元のフィーチャクラス名
    in_feature = "LineSplit"

    #接合先のフィーチャクラス名
    out_feature = "LineUnsplit"

    #インスタンス化
    unsplitLine = UnsplitLine(workspace,in_feature,out_feature)

    #接合対象のラインを抽出
    lines = unsplitLine.search_line()

    #対象のラインを接合
    unsplitLine.create_line(lines)


このソースをArcMap上で実行した結果が以下になります。
f:id:sanvarie:20160319171524p:plain


ちゃんと作られました。このように一本ずつだったラインが
f:id:sanvarie:20160319171549p:plain


このように接合されました!
f:id:sanvarie:20160319171605p:plain


プログラムの中で「UNSPLIT_LENGTH」というカラムを作成し、接合前のShape_Lengthを格納するようにしました(接合対象のラインが三本あり、それぞれの長さが10,20,30だった場合、60が格納される)。ただ、ちょっと誤差がでてしまっているので、接合結果を比較する場合は、小数点以下をはずすなどしないとだめそうですね・・・
f:id:sanvarie:20160319171709p:plain


出力前後の結果確認用にテキストも出力するようにしています(今回はC:\ArcPySampleの直下に作成)。こちらも微妙に誤差が出ていますね・・・
f:id:sanvarie:20160319171839p:plain

課題

  1. 接合前後でShape_Lengthに微妙に誤差がでてしまっている
  2. 大量のデータを処理する場合、パフォーマンスに難あり
  3. ソースが煩雑になりすぎた(ちょっと色々無理矢理やってしまったところがあります・・・)。共通化など、色々整理が必要ですね。
  4. バグが大量にありそうです・・・

課題は多いものの、とりあえずAdvancedの機能であるライン接合処理ができるようになったことは素直にうれしいですね。自分も色々勉強しながらだったので、今度作るときはもっとうまくプログラムを書くことができそうです。

以上、本日は「【PythonGIS】ArcPyでライン接合を行ってみた ~Advancedライセンス?!だが断る!~」でした。

ArcPyを使用してラインの頂点数を取得する

さて、本日はArcPyを使用してラインの頂点数を取得してみようと思います。頂点数が多いとパフォーマンスにも影響が出てしまうため、余分な頂点数がありそうなラインを抽出するのに役立つかと思います。

実行環境

事前準備

  1. 「Line」フィーチャクラスを作成します。「頂点数」格納用のカラムも用意します。
  2. 「Line」フィーチャクラスにラインを作図します。

f:id:sanvarie:20160302133304p:plain

f:id:sanvarie:20160302142157p:plain

サンプルコード

ラインの頂点数を取得して属性テーブルに格納するサンプルです。

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

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

#対象のフィーチャクラス
infc = "Line"

dictTop = {}
cnt = 0

#ラインフィーチャの頂点数を検索
for row in arcpy.da.SearchCursor(infc, ("OID@","SHAPE@")):
    for part in row[1]:
        for pnt in part:
            if pnt:
                cnt += 1
            else:
                print("error")
                continue
        #ディクショナリにOIDと頂点数を格納
        dictTop[row[0]] = cnt
        cnt = 0
del row

#各フィーチャの頂点数を更新
for i in dictTop:
    cursor = arcpy.da.UpdateCursor(infc, ("OID@","Top"))
    for row in cursor:
        if i == row[0]:
            row[1] = dictTop[i]
            cursor.updateRow(row)
del row

こんな感じになれば成功です。今回は属性に頂点数を格納しましたが、CSV出力したりしてもいいかもしれませんね。
f:id:sanvarie:20160302142614p:plain

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

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

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

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

ダウンロード

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

ダウンロードページ>Large scale data, 1:10m>Cultural からダウンロードします。
f:id:sanvarie:20160228162730p:plain

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

実行環境

事前準備

  • 「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)

こんな感じになれば成功です。
f:id:sanvarie:20160229111608p:plain

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

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

対象サイト

外務省のサイトがよさげだったので今回はこのサイトから情報を取得しようと思います。世界の国々 | 外務省
ただし、世界と日本のデータを見る(世界の国の数,国連加盟国数,日本の大使館数など) | 外務省で確認すると、「世界の国の数=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])

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

サンプルコード④

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だと思います。
f:id:sanvarie:20160302101018p:plain

課題

  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)。f:id:sanvarie:20160228113422p:plain
  4. 「HighShcool」フィーチャクラスを作成します。

f:id:sanvarie:20160228105537p:plain

f:id:sanvarie:20160228105618p:plain

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

「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)

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

「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)

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

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

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

サンプルコード

各県ごとに付与されたポイントを以下レンジでレンダリングを行います。数値が大きいほど、その県が強いということがわかりますね。
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()

このような感じになりました。
f:id:sanvarie:20160228132129p:plain

結果の分析

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

f:id:sanvarie:20160228134230p:plain

課題

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

☆番外編☆

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

テーブルの作成

「HighSchool_Record」テーブルを作成します。
f:id:sanvarie:20160228141940p:plain

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

サンプルコード

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

■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

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

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

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

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