読者です 読者をやめる 読者になる 読者になる

GIS奮闘記

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

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の日本地図(全国市区町村界データ)を加工してみる」でした!

【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を使用して詳細地図の秘密に迫る~」でした!

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~」でした。

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

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

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

 

■対象者

ArcGISに興味のある方全員

 

■対象技術

・ArcMap

・ArcObjects

・ArcPy

 ・その他

 

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

 

■今後の展望

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

ArcGISハッカソンの開催

 

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

 

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

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

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

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

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

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

このページをいったんHTMLとしてローカルに保存します(動的に作られているページなので)。そのHTMLの表をスクレイピングしてみます。

インストール

pandasだけでなくlxml、html5lib、BeautifulSoup4もインストールする必要があります。

サンプルコード①

HTMLの表をスクレイピングして結果をCSVに出力するサンプルです。

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

#保存したHTML
html = 'D:\python\earthquake\earthquake.html'

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

#表の部分を取得
table = pd.DataFrame(dataframes[5])
table.columns = ['No','Time','Name','Latitude','Longitude','Depth','Magnitude','SeismicIntensity']

#CSV出力
table.to_csv('D:\python\earthquake\earthquake.csv' ,encoding = 'Shift-JIS')

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

データ準備

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

ArcMapで読み込みます。
f:id:sanvarie:20160126122232p:plain

ポイントのフィーチャークラスを作成します。
f:id:sanvarie:20160126122349p:plain

フィーチャークラスの属性を以下の地理座標系にします。
f:id:sanvarie:20160126122421p:plain

サンプルコード②

「サンプルコード①」で出力したCSVを元にマップ上にポイントをプロットするサンプルです。色々やっているので、要点だけ書いておきます。

1.CSVから読み込んだ緯度経度を60進数から10進数に変換
2.「1」の緯度経度からポイントをプロット
3.「earthquake」フィーチャクラスにカラム追加(Arcpy)
4.「3」のカラムにCSVから読み込んだ情報を付与(Arcpy)

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

pointGeometryList = []
def main():

    #CSV読込
    table = pd.read_csv(r"D:\python\earthquake\earthquake.csv",encoding="SHIFT-JIS")

    data_list = []
    for key,rowD in table.iterrows():
        if key > 0 :
            #緯度経度を度分秒に分解
            latDo = rowD.Latitude[0:rowD.Latitude.find(u"°")]
            lonDo = rowD.Longitude[0:rowD.Longitude.find(u"°")]
            latFun = rowD.Latitude[rowD.Latitude.find(u"°")+1:rowD.Latitude.find(u".")]
            lonFun = rowD.Longitude[rowD.Longitude.find(u"°")+1:rowD.Longitude.find(u".")]
            latByou = rowD.Latitude[rowD.Latitude.find(u".")+1:rowD.Latitude.find(u"′")]
            lonByou = rowD.Longitude[rowD.Longitude.find(u".")+1:rowD.Longitude.find(u"′")]

            #緯度経度を60進数から10進数に変換
            lat,lon = latLong(latDo,lonDo,latFun,lonFun,latByou,lonByou)
            data_list.append([rowD.No,rowD.Time,rowD.Name,lat,lon,rowD.Depth,rowD.Magnitude,rowD.SeismicIntensity])
            CreatePont(lat,lon)
    table2 = pd.DataFrame(data_list)
    table2.columns = ["No","Time","Name","Latitude","Longitude","Depth","Magnitude","SeismicIntensity"]
    arcpy.CopyFeatures_management(pointGeometryList, r"C:\ArcPySample\test.gdb\earthquake")

    arcpy.env.workspace = ur"C:\ArcPySample\test.gdb"

    #earthquakeレイヤにカラムを追加
    arcpy.AddField_management("earthquake", "Time", "Text")
    arcpy.AddField_management("earthquake", "Name", "Text")
    arcpy.AddField_management("earthquake", "lat", "Double")
    arcpy.AddField_management("earthquake", "lon", "Double")
    arcpy.AddField_management("earthquake", "Depth", "Text")
    arcpy.AddField_management("earthquake", "Magnitude", "Text")
    arcpy.AddField_management("earthquake", "SIntensity", "Text")

    for key,rowD in table2.iterrows():
        cursorJ = arcpy.UpdateCursor("earthquake")
        for rowJ in cursorJ:
            if rowD.No == rowJ.OBJECTID:
                rowJ.setValue("Time", rowD.Time)
                rowJ.setValue("Name", rowD.Name)
                rowJ.setValue("lat", rowD.Latitude)
                rowJ.setValue("lon", rowD.Longitude)
                rowJ.setValue("Depth", rowD.Depth)
                rowJ.setValue("Magnitude", rowD.Magnitude)
                rowJ.setValue("SIntensity", rowD.SeismicIntensity)
                cursorJ.updateRow(rowJ)

#ポイントを配置
def CreatePont(lat,lon):
    global pointGeometryList

    point = arcpy.Point()
    point.X = lon
    point.Y = lat

    #日本地図の空間参照を利用
    dataset = r"D:\python\earthquake\japan_ver80.shp"
    spatial_reference=arcpy.Describe(dataset).spatialReference
    pointGeometry = arcpy.PointGeometry(point,spatial_reference)
    pointGeometryList.append(pointGeometry)

#緯度経度を60進数から10進数に変換
def latLong(latDo,lonDo,latFun,lonFun,latByou,lonByou):
    lat1 = float(latDo)
    lat2 = float(latFun) / 60.0
    if type(latByou) == float:
        lat3 = float(latByou) / 60.0 / 60.0
    else:
        lat3 = float(0)
    lat1 = lat1 + lat2 + lat3

    long1 = float(lonDo)
    long2 = float(lonFun) / 60.0
    if type(lonByou) == float:
        long3 = float(lonByou) / 60.0 / 60.0
    else:
        long3 = float(0)
    long1 = long1 + long2 + long3
    return lat1,long1

if __name__ == '__main__':
    main()

ポイントがプロットされました。位置もあっているっぽいですね。
f:id:sanvarie:20160126132631p:plain

ポイントに属性も入っていますね。
f:id:sanvarie:20160126141215p:plain

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

サンプルコード③

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

# -*- coding: utf-8 -*-
import arcpy
mxd = arcpy.mapping.MapDocument("current")
lyr = arcpy.mapping.ListLayers(mxd, "earthquake")[0]
if lyr.symbologyType == "UNIQUE_VALUES":
  lyr.symbology.valueField = "SIntensity"
  lyr.symbology.addAllValues()
arcpy.RefreshActiveView()
arcpy.RefreshTOC()
del mxd

ポイントの色分けができました。次にシンボルの大きさを分けようと思ったのですが、どうやらこれもArcpyではできない模様。うー貧弱すぎる・・・
f:id:sanvarie:20160126141932p:plain

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

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

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

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

【ArcMapとArcpy】オープンデータを使って自分だけの地図を作ろう! ~阪神・淡路大震災編~

さて、本日は神戸市のオープンデータを利用してみようと思います。阪神・淡路大震災「1.17の記録」阪神・淡路大震災の記録写真が公開されています。前回行ったようにこの写真にジオタグを追加して、その緯度経度にポイントを配置してみます。

データのダウンロード

上記サイトにおいて、画像データの一括ダウンロード>CSVでダウンロードを行います。
こんな感じのCSVがダウンロードできました。緯度経度の項目があるにも関わらず、ここは空白になっています・・・仕方ないので以前行ったジオコーディングを使用したいと思います。ただ、場所が空白、もしくは、ジオコーディングできない(衣掛町5丁目2「こうべ地震災害対策広報」みたいな記述)写真は取込対象外とします。
f:id:sanvarie:20160121082028p:plain

ポイントとなる技術

主に以下のような技術を使用します。
①画像のダウンロード
②画像にジオタグ(geotag)の追加
③ジオコーディング
ジオタグを追加した画像をポイントに変換

日本地図のダウンロード

こちらも以前行ったようにESRIのサイトからダウンロードを行ってください。
全国市区町村界データ | 製品 | ESRIジャパン

f:id:sanvarie:20160121082956p:plain

サンプルコード

上述した「ポイントとなる技術」を盛り込んだサンプルコードです。パスなどはご自身の環境にあったものに書き換えてください。

# -*- coding: utf-8 -*-
import pandas as pd
import urllib2
import cookielib
import os
import pyexiv2
from PIL import Image
import geocoder
import arcpy

def main():
    #ダウンロードしたcsvの読込
    df = pd.read_csv(r"D:\python\geotag\20160120142702.csv",encoding="SHIFT-JIS")

    #各種設定
    folder = "D:\python\geotag"
    outFeatures = r"C:\ArcPySample\ArcPySample.gdb\geotag"
    badPhotosList = r"C:\ArcPySample\ArcPySample.gdb\empty"
    photoOption = "ONLY_GEOTAGGED"
    attachmentsOption = ""

    cj = cookielib.CookieJar()
    opener = urllib2.build_opener(urllib2.HTTPCookieProcessor(cj))

    #CSVの行数分ループ
    for key,rowD in df.iterrows():

        #リスエスト
        req = urllib2.Request(rowD.image)

        #出力先を指定
        with open(os.path.join(folder,rowD.filename), 'wb') as f:
            f.write(opener.open(req).read())

        #場所が無いデータは無視
        if type(rowD.place) != float:

            #緯度経度取得
            loc = get_coordinate(rowD.place)

            #井戸経度が取得できた場合
            if loc.lat is not None:
                #ジオタグを画像に追加
                set_gps_location(os.path.join("D:\python\geotag",rowD.filename),loc.lat,loc.lng)
            loc = None

    #ジオタグ付き写真をポイントに変換
    arcpy.GeoTaggedPhotosToPoints_management(folder, outFeatures, badPhotosList, photoOption, attachmentsOption)

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

    return ret

def to_deg(value, loc):
        if value < 0:
            loc_value = loc[0]
        elif value > 0:
            loc_value = loc[1]
        else:
            loc_value = ""
        abs_value = abs(value)
        deg =  int(abs_value)
        t1 = (abs_value-deg)*60
        min = int(t1)
        sec = round((t1 - min)* 60, 5)
        return (deg, min, sec, loc_value)

def set_gps_location(file_name, lat, lng):

    try:
        lat_deg = to_deg(lat, ["S", "N"])
        lng_deg = to_deg(lng, ["W", "E"])

        # 緯度、経度を10進法→60進法(度分秒)に変換
        exiv_lat = (pyexiv2.Rational(lat_deg[0]*60+lat_deg[1],60),pyexiv2.Rational(lat_deg[2]*100,6000), pyexiv2.Rational(0, 1))
        exiv_lng = (pyexiv2.Rational(lng_deg[0]*60+lng_deg[1],60),pyexiv2.Rational(lng_deg[2]*100,6000), pyexiv2.Rational(0, 1))
        metadata = pyexiv2.ImageMetadata(file_name)
        metadata.read()

        #メタデータを付与
        metadata["Exif.GPSInfo.GPSLatitude"] = exiv_lat
        metadata["Exif.GPSInfo.GPSLatitudeRef"] = lat_deg[3]
        metadata["Exif.GPSInfo.GPSLongitude"] = exiv_lng
        metadata["Exif.GPSInfo.GPSLongitudeRef"] = lng_deg[3]
        metadata["Exif.Image.GPSTag"] = 654
        metadata["Exif.GPSInfo.GPSMapDatum"] = "WGS-84"
        metadata["Exif.GPSInfo.GPSVersionID"] = '2 0 0 0'

        #メタデータを上書き
        metadata.write()
    except:
        print u"エラー。たまに出る。とりあえず飛ばす"

if __name__ == '__main__':
    main()

結果を確認してみます。それっぽい箇所にポイントが配置されましたが、全然関係ない場所にも配置されてしまっています。これは例えば「若松町」という場所が他の都道府県にも存在するという原因かと思います。完璧に行う場合はもう少し工夫が必要ですがとりあえず気にせず進みます(笑)
f:id:sanvarie:20160121090210p:plain

うまくいったっぽい箇所をGoogleMapと比較して確認してみます。「ポートターミナル周辺」でジオコーディング、GoogleMapで検索した結果です。
f:id:sanvarie:20160121085744p:plain

f:id:sanvarie:20160121085813p:plain

大体あってますかね。ちなみにこんな写真です。ジオタグもちゃんと追加されていますね。
f:id:sanvarie:20160121085942j:plain

f:id:sanvarie:20160121090021p:plain

ArcMapでポイントの個別属性から画像を参照することができます。
f:id:sanvarie:20160121090413p:plain

f:id:sanvarie:20160121090426p:plain

これだけで終わるのもあれなので、ここで一工夫してみようと思います。まずフィールドの追加を行います。

arcpy.AddField_management("geotag", "HTML_Path", "String")

こうなったらOKです。
f:id:sanvarie:20160121090941p:plain

その後、以下のようにフィールド演算を行います。

codeblock = """
def setHTMLPath(Path):
    return "<img src='"  + Path + "'" + "Width='300'>"
"""

expression = "setHTMLPath(!Path!)"
arcpy.CalculateField_management("geotag", "HTML_Path", expression,"PYTHON_9.3",codeblock)

こうなったらOKです。
f:id:sanvarie:20160121091827p:plain

ArcMapでHTMLポップアップ機能を使用します。
f:id:sanvarie:20160121091904p:plain

任意のポイントを選択するとポップアップが起動します。属性にもっと詳細な情報を入れたらより充実したマップになりますね。
f:id:sanvarie:20160121091958p:plain

今回わかったことは
①ダウンロードしたCSVの情報がかなり大雑把。そのため、ジオコーディングに一工夫が必要。
ジオタグ追加→ポイント変換のコンボはかなり面白い!GPSの精度がさらに向上すれば、こういった案件も増えてきそうな予感。
③写真を見るとあらためて震災の怖さというものを感じた。

もっとオープンデータが充実すればさらに精度の高いマップを作成することができると思うので、今後そこを期待したいですね。

以上、「【ArcMapとArcpy】オープンデータを使って自分だけの地図を作ろう! ~阪神・淡路大震災編~ 」でした。