GIS奮闘記

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

スポンサーリンク

【PythonとGIS】ArcPyでライン接合を行ってみた

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

以上、本日は「【PythonとGIS】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的視点で高校サッカーを分析してみる」でした。

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

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

カラム追加

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

サンプルコード

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

カラムが追加されましたね。
f:id:sanvarie:20160226081351p:plain

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

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

インストール

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)

人口、世帯数を持った県ごとのポリゴンの作成に成功しました。
f:id:sanvarie:20160226082040p:plain

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

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

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

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

データ

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

f:id:sanvarie:20160209161731p:plain

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

■CHARACTER
f:id:sanvarie:20160209162348p:plain

■LANDMARK
f:id:sanvarie:20160209162621p:plain

■PIPE
f:id:sanvarie:20160209162650p:plain

■ROAD
f:id:sanvarie:20160209162730p:plain

インストール

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

結果を確認します。
f:id:sanvarie:20160209170233p:plain

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

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

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

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

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

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

変換ツールの使用

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

ArcMapにこのようなツールを追加することができました。
f:id:sanvarie:20160205200056p:plain

これで変換ですね。
f:id:sanvarie:20160205204355p:plain

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

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

実行した結果、そうそうこれ。これですよ!ESRIさんありがとうございますm(_ _)m
f:id:sanvarie:20160205200551p:plain

データの内容

今回一番欲しい建物のデータ定義はこんな感じになっています。
f:id:sanvarie:20160205200904p:plain

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

サンプルコード

各テーブルやフィーチャクラスを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Ⅱの各テーブル、フィーチャクラスと同じにしました。

■大字
f:id:sanvarie:20160205202107p:plain

■字丁目
f:id:sanvarie:20160205202151p:plain

■街区
f:id:sanvarie:20160205202224p:plain

■建物
f:id:sanvarie:20160205202255p:plain

どのテーブルもコードばかりでわけがわかりませんね。定義書を見た感じ以下を考慮した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に出力したりして利用することができますね!
f:id:sanvarie:20160205202716p:plain

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

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