GIS奮闘記

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

【PythonとGIS】ESRIデータコレクションを分析する ~ArcPyを使用して詳細地図の秘密に迫る~

さて、本日はESRIデータコレクションを分析してみようと思います。「ESRIデータコレクション」とはESRIさんが販売している背景地図のことです。以前、ESRIの方とお話しした際に「データコレクションは涙ぐましい努力をして少しでも速度が出るようにしています」とお聞きしました。実際にそのデータをいじってみると、信じられないくらい早いというのが正直な感想でした。そして、ArcGISは何だか描画が遅いと感じていた私はデータコレクションの中身を知れば速度改善の方法がわかるのでは?!と思い、今回の記事を書くに至りました。

本記事で紹介するサンプルコードはGistにも公開しました。ESRIデータコレクションの詳細地図.lyrの情報をCSV出力するスクリプトです · GitHub

データコレクションの種類

以下のように色々あります。詳細はArcGIS データコレクション | 製品 | ESRIジャパンを参照してください。

・スタンダードパック
・詳細地図
・道路網
・住居レベル住所
・地形
・統計

今回はその中でも「詳細地図」について分析をしてみようと思います。

詳細地図とは

「詳細地図」とはESRIさん曰く、株式会社ゼンリンの地図データと国土交通省基盤地図情報などを加工して開発された背景地図データベースのことらしいです。ちょっとざっくりすぎるので、仕様書を参考にしてひも解いてみようと思います。
http://www.esrij.com/cgi-bin/wp/wp-content/uploads/documents/Shosai_Data_Spec.pdf

データソース

以下を元に「詳細地図」が作られているみたいですね。

①株式会社ゼンリン Zmap-AREA II(2014-3 版)
→Zmap-TOWNIIでないのがちょっと残念(私は戸番データをよく使うので)

②株式会社ゼンリン 地図データ
→一体何を指すのか?!うーん、わからん。。

国土交通省 国土地理院 数値地図(国土基本情報)
→恥ずかしながら使ったことがありません。お金が必要な国土地理院のデータですね。「都道府県界」「鉄道路線」などに使っているみたいです。

国土交通省 国土地理院 基盤地図情報(数値標高モデル)
→以前ブログで紹介しましたね。Pythonで基盤地図情報の数値標高モデルを解析してGeoTiffに変換する(横須賀編) ~Numpy、GDAL、Arcpyを駆使して横須賀の標高(地盤高)を視覚化する~ - GIS奮闘記

総務省統計局 平成 22 年国勢調査
→統計情報として使用しているみたいです。「人口集中地区」に使用しているみたいですね。

データ形式

色々用意されていますが、今回は「レイヤー ファイル形式(.lyr) 」を分析します。
①地図データベース・・・ファイル ジオデータベース形式
②地図に関する設定・・・ArcMap ドキュメント形式(.mxd)、マップ ファイル形式(.mapx)、レイヤー ファイル形式(.lyr)

構成

以下のような構成になっています。上述したように今回はレイヤファイル、つまり、詳細地図.lyrを分析します。

f:id:sanvarie:20160213153013p:plain

フィーチャクラス

「詳細地図.lyr」は三つのgdbを参照しています。今回は①②について分析します。

①L_SCALE.gdb・・・大縮尺で地図表示を行うためのフィーチャクラスを格納しています。
②MS_SCALE.gdb・・・小〜中縮尺で地図表示を行うためのフィーチャクラスを格納しています。
③PLACE_FINDER.gdb・・・郵便番号、電話番号などを検索できるロケーターを収録しています。

詳細地図.lyrの構成

以下グループレイヤにまとめられています。

①詳細図
②広域地図
③背景

さらにグループレイヤの中にグループレイヤがあり、けっこう複雑な構成になっています。今までこういったやり方をしたことがなかったので、非常に勉強になりますね。
f:id:sanvarie:20160213153928p:plain

実際に見るとこんな感じですね。
f:id:sanvarie:20160213154129p:plain

仕様書を見ると各フィーチャクラスの情報が載っているのですが、これだけだと分析には物足りないですね(しかもPDFですし)。ですので、以前行ったように「詳細地図.lyr」の情報をCSVに出力しようと思います。

インストール

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

サンプルコード

「詳細地図.lyr」の情報をCSVに出力するサンプルコードです。

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

fields = ""
dataSource = ""
subfields = ""

def output_csv(info,f,flg,subTypeCd):
    global dataSource
    global fields
    global subfields

    if flg == 0:
        #フィーチャをnumpyに変換してpandasに格納
        numFeature = arcpy.da.FeatureClassToNumPyArray(dataSource,("OID@"),null_value=-9999)
        dfFeature = pd.DataFrame(numFeature)
        outputstring = info[0] \
                       + "," + info[1] \
                       + "," + info[2] \
                       + "," + info[3] \
                       + "," + info[4] \
                       + "," + info[5] \
                       + "," + str(info[6]) \
                       + "," + str(info[7]) \
                       + "," + "" \
                       + "," + str(len(dfFeature)) + "\n"

        f.write(outputstring.encode("SHIFT-JIS"))
    else:
        #フィーチャをnumpyに変換してpandasに格納
        numFeature = arcpy.da.FeatureClassToNumPyArray(dataSource,(fields),null_value=-9999)
        dfFeature = pd.DataFrame(numFeature)

        #サブタイプごとのデータを抽出
        if fields == "LAYERCODE":
            dfFeatureValue = pd.DataFrame(dfFeature[dfFeature.LAYERCODE == subTypeCd])
        elif fields == "TYPE":
            dfFeatureValue = pd.DataFrame(dfFeature[dfFeature.TYPE == subTypeCd])

        outputstring = info[0] \
                       + "," + info[1] \
                       + "," + info[2] \
                       + "," + info[3] \
                       + "," + info[4] \
                       + "," + info[5] \
                       + "," + str(info[6]) \
                       + "," + str(info[7]) \
                       + "," + subfields \
                       + "," + str(len(dfFeatureValue)) + "\n"

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

def get_layer_info(layer,f,group1,group2,group3,group4):
    global fields
    global dataSource
    global subfields

    info = []
    dataSource = layer.dataSource

    #フィーチャレイヤの場合
    if layer.isFeatureLayer:

        #フィーチャタイプ取得
        fType = arcpy.Describe(dataSource).FeatureType

        #アノテーションの場合(アノテーションのジオメトリタイプはポリゴンになってしまうのでこの処理を追加)
        if fType == "Annotation":
            gType = "Annotation"
        else:
            #ジオメトリタイプ取得
            gType = arcpy.Describe(dataSource).shapeType
    else:
        gType = arcpy.Describe(dataSource).format

    info.append(group1)
    info.append(group2)
    info.append(group3)
    info.append(group4)
    info.append(layer.datasetName)
    info.append(gType)
    info.append(layer.maxScale)
    info.append(layer.minScale)

    #フィーチャレイヤの場合
    if layer.isFeatureLayer:
        #サブタイプを取得
        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 == "":
                        output_csv(info,f,0,"")

                    #サブタイプが設定されている場合
                    else:
                        for stkey in list(groupDict.keys()):
                            if stkey == 'Name':
                                subfields = groupDict[stkey]
                                subDict[subfields] = name
                                output_csv(info,f,1,subDict[subfields])
    else:
        output_csv(info,f,0,"")


if __name__ == '__main__':
    #レイヤファイルのパスを指定
    filepath = ur"D:\python\datacollection\詳細地図\詳細地図.lyr"

    #CSV出力先
    f = open(r"D:\python\datacollection/datalist.csv", "w")

    # CSVにヘッダを出力
    outputstring = u"グループ1,グループ2,グループ3,グループ4,フィーチャクラス,ジオメトリタイプ,最大縮尺,最小縮尺,サブタイプ,アイテム数\n"
    f.write(outputstring.encode("SHIFT-JIS"))

    #レイヤオブジェクト取得
    layer = arcpy.mapping.Layer(filepath)
    for lay in layer:
        if lay.name == u"詳細図" or lay.name == u"広域地図" or lay.name == u"背景":
            if lay.isGroupLayer:
                for ly in lay:
                    if ly.isGroupLayer:
                        for l in ly:
                            if l.isGroupLayer:
                                for ll in l:
                                    get_layer_info(ll,f,lay.name,ly.name,l.name,ll.name)
                            else:
                                get_layer_info(l,f,lay.name,ly.name,l.name,"")
                    else:
                        if ly.longName.find("\\") == -1:
                            get_layer_info(ly,f,lay.name,ly.name,"","")


    f.close()

結果を確認します。各種情報を出力することができましたね!それぞれのジオメトリタイプやサブタイプごとのアイテム数などを出力してみました。ぱっと見、仕様書通りに出てるっぽいですが、間違っていたらすみません(笑)。

グループ1~グループ4とありますが、その中で最後に表示されているものがレイヤ名です(例えばグループ3まで埋まっていて、グループ4が空の場合はグループ1、2がグループレイヤ名、3がレイヤ名となります)。

f:id:sanvarie:20160213155030p:plain

分析

さて、一体どんな秘密が隠されているのか。

アノテーションがない
「L_ANNOCHO」とかアノテーションっぽいフィーチャクラスがありますが、これは実は全部ポイントで、ラべリングしているみたいですね(私は業務の都合上、泣く泣く大量のアノテーションを使っています)。

f:id:sanvarie:20160213160054p:plain

②スケールの設定が細かい
とにかくスケールの設定が細かいですね。ここをもっと研究すべきかもしれません。ざっくりですが、以下のような感じです。当たり前のことですが、拡大時=フィーチャ数多、縮小時=フィーチャ数少、と言う傾向が見られます。やはり基本に忠実にということですね。L_SCALE.gdb、MS_SCALE.gdbをうまく使い分けている気がします。にしてもフィーチャ数が百万以上あってもさくさく動いてくれるのは感動的ですね。

f:id:sanvarie:20160213164707p:plain

③サブタイプを設定
これも当たり前なのかもしれませんが、サブタイプを上手く活用している印象があります。

f:id:sanvarie:20160213163213p:plain

ざっとこれくらいでしょうか。私の分析力の問題かもしれませんが、特別なことをしているようには見えませんでした。ただ、スケールの設定はものすごく緻密だと感じました。この辺をもっとシビアにすることによってスピードアップが出来そうな気がします。

☆驚いたこと☆

建物に影があるなぁと思ってみてみたら、これは「目標物・一般家枠」と「目標物・一般家枠(影)」というレイヤを重ねているだけ、ということがわかりました。なるほど、と言う感じですね。ただ、フィーチャ数が増えるからちょっと厳しいかもしれないですね。。。

f:id:sanvarie:20160213165350p:plain

かなり駆け足でしたが、ESRIデータコレクションというのがどんなものか、というのが少しだけわかった気がします。今後はさらに深く、また、他のデータも分析してみたいと思います。

以上、「【PythonGISESRIデータコレクションを分析する ~ArcPyを使用して詳細地図の秘密に迫る~」でした!