PythonとJSMで東証一部上場企業を可視化してみる

ものすごく久しぶりのブログ更新になります。さぼり気味で元々しょぼい腕前がさらに錆びついた気が・・・今後はまた定期的に更新しようと思いますので、どうぞよろしくお願いします。さて、今回は「PythonとJSMで東証一部上場企業を可視化してみる」です。たくさんの企業が東証一部に上場していることはご存知かと思いますが、どの地域に一部上場企業が集まっているかはあまり意識したことがないのではないでしょうか？今回は上場企業の本社所在地に対してポイントをプロットしてみようと思います。とはいってもやっぱりほとんど東京のような気もしますが。。。

データ

データは一体どこからとってくるの？ってことになると思いますが、「Yahoo!ファイナンス」からとってくることができます。「jsm」というライブラリを使えば簡単に取得できます。作っていただいた方に感謝ですね。使用方法などはPyPIに記述されているので、そちらをご覧になってください。jsm 0.19 : Python Package Index

インストール

jsm, Pandas、geocoderのインストールをお願いします。また、今回はPandasを使用してスクレイピングをしますので、lxml、html5lib、BeautifulSoup4も必要です。インストールをお願いします。

環境

Windows7 64bit
ArcGIS10.2.0
Python2.7.3

事前準備

1.Company.gdbを作成します。
2.前回作成した日本地図を「Company.gdb」にコピーします（フィーチャクラス名は「Japan」）。
3.CompanyInfoフィーチャークラスを作成します。カラムは以下のように作成します。

このような感じになるかと思います。

株式市場の種類について

参考までに

日本の株式市場は、4つの証券取引所で成っています。今回は東証一部に限定しましたが、株式市場に上場している企業すべてを分析してみても面白そうですね。

・東京証券取引所
1部・2部・マザーズ・ジャスダック・上場投信信託（ETF）・不動産投資信託（J-REIT）

・名古屋証券取引所
1部・2部・セントレックス

・札幌証券取引所
上場株式市場・アンビシャス

・福岡証券取引所
上場株式市場・不動産投資信託（J-REIT）・Q-Board

サンプルコード

東証一部上場企業の本社所在地にポイントをプロットするサンプルコードです

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

#定数
INDUSTRY_CODE = "INDUSTRY_CODE"
INDUSTRY_NAME = "INDUSTRY_NAME"
BRAND_CODE = "BRAND_CODE"
BRAND_MARKET = "BRAND_MARKET"
BRAND_NAME = "BRAND_NAME"
LOCATION = "LOCATION"

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

arcpy.env.workspace = r"C:\ArcPySample\Company.gdb"

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

    return ret

def create_point(df_merge):

    i = 0
    with arcpy.da.InsertCursor(infc, ["OID@","SHAPE@",INDUSTRY_CODE,INDUSTRY_NAME,BRAND_CODE,BRAND_MARKET,BRAND_NAME,LOCATION]) as cursor:

        for _,row in df_merge.iterrows():

            if row.LOCATION:
                loc = get_coordinate(row.LOCATION) #ジオコーディング

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

                    cursor.insertRow((oid,pointGeometry,row.INDUSTRY_CODE,row.INDUSTRY_NAME,row.BRAND_CODE,row.BRAND_MARKET,row.BRAND_NAME,row.LOCATION))
                    i = i + 1

def get_finace_data():

    q = jsm.Quotes()
    b = jsm.Brand()
    IDS = b.IDS

    data_list = []

    for industry_code in IDS.keys():

        industry_name = IDS[industry_code]
        brand_data = q.get_brand(industry_code) #企業情報を取得

        for brand in brand_data:
            if brand.market == u"東証1部":
                data_list.append([industry_code,industry_name,brand.ccode,brand.market,brand.name])

    brand_data_table = pd.DataFrame(data_list) #データフレームに格納
    brand_data_table.columns = [INDUSTRY_CODE,INDUSTRY_NAME,BRAND_CODE,BRAND_MARKET,BRAND_NAME]

    get_location(brand_data_table) #本社所在地を取得

def get_location(brand_data_table):

    location_list = []
    market_code = ""

    for _ ,row in brand_data_table.iterrows():

        #東証一部に限定したのでこの分岐は不要だが一応残しておく
        #if row.brand_market.find(u"札") > -1:
        #    market_code = u"S"
        #elif row.brand_market.find(u"名") > -1:
        #    market_code = u"N"
        #elif row.brand_market.find(u"福岡") > -1:
        #    market_code = u"F"
        #else:
        #    market_code = u"T"

        url = "http://stocks.finance.yahoo.co.jp/stocks/profile/?code=%s.T" % row.BRAND_CODE
        print url
        fetched_dataframes = pd.io.html.read_html(url) # pandasで本社所在地をスクレイピング

        try:
            location_data = fetched_dataframes[1][1][2] # 本社所在地の要素を取得

            if location_data.find(u"[") > -1: # [周辺地図]もしくは[主要事業所]という文言が本社所在地に含まれてしまっている会社があるので
                location_data = location_data[:location_data.find("[")]

            if location_data.find(u"〒") > -1: #郵便番号ははずす
                location_data = location_data[10:]
        except Exception as e:
            print "error"
        else:
            location_list.append([row.BRAND_CODE,location_data])

    location_data_table = pd.DataFrame(location_list) #データフレームに格納
    location_data_table.columns = [BRAND_CODE,LOCATION]

    df_merge = pd.merge(brand_data_table, location_data_table, on=[BRAND_CODE],how='left') #データフレームを結合

    create_point(df_merge) #ポイントをプロット

if __name__ == '__main__':
    get_finace_data()

このようにポイントがプロットされたら成功です！

属性もばっちり入っていますね。

やはり予想通り、東京に集中していますね。あとは神奈川、大阪などの都市圏にもかなり存在しているようです。そりゃたくさん人がいるわけだぁ。

「jsm」は株価もとってこれるので、東京の地価と企業の時価総額の関連性などを調べてみても面白そうですね。あとは東証一部だけではなく、全ての上場企業のデータを抽出し、都道府県ごとの人口と時価総額を比べてみても面白そうですね。以上、今回は「PythonとJSMで東証一部上場企業を可視化してみる」でした！

ArcPyレシピ集④ ～フィーチャクラスの座標系を取得～

さて、また間が空いてしまいましたが、今回もArcPyレシピを書いてみたいと思います。テーマは「フィーチャクラスの座標系の取得です」。例えば、GDBに複数のフィーチャクラスが入っている際、座標系が異なるフィーチャクラスが入ってしまっていた、というご経験があるのではないでしょうか。今回のレシピはそういったトラブルを回避するために役立ちます。

環境

Windows7 64bit
ArcGIS10.2.0
Python2.7.2

データ

• GDBの作成
• フィーチャクラスの作成

→座標系は右記にしました（平面直角座標系9系(JGD2000)、平面直角座標系8系(JGD2000)、平面直角座標系7系(JGD2000)）

サンプルコード

対象のGDBに存在するフィーチャクラスの座標系のコードとフィーチャクラス名を多次元リストで取得するサンプルです。座標系のコードはこのサイトをご参照願います。
GISのための測地成果、測地系、楕円体、投影座標系、EPSGコードのまとめ - 自然環境保全のための周辺技術

# -*- coding: utf-8 -*-
import arcpy
from os.path import join

def get_projection_code(path):
    """フィーチャクラスの座標系のリストを取得します

    arguments:
    path -- 任意フォルダを指定してください。
　　　　　　指定したフォルダにprojection.csvを出力します。

    """
    arcpy.env.workspace = path
    projection_list = []

    try:

        for gdb in arcpy.ListWorkspaces():

            arcpy.env.workspace = join(path,gdb)

            if gdb.find('Raster') > -1:

                featureclasses = arcpy.ListDatasets()

            else:

                featureclasses = arcpy.ListFeatureClasses()

            for fc in featureclasses:

                projection = arcpy.Describe(fc).spatialReference.factoryCode

                if projection not in projection_list:
                    projection_list.append([fc,projection]) #ラスタの場合でも、ラスタ座標系ではなくXY座標系にしているので要注意

        return projection_list

    except Exception as e:
        raise

if __name__ == '__main__':
    lists = get_projection_code(ur"D:\test")

    for l in lists:
        print l

結果を見てみると、見た目は文字コードの問題でこのようになっていますが（見た目だけのものなので、特に問題なし）、ばっちり座標系のコードは取得できました。自分はよく座標系の設定をミスってしまうので、こういうチェックとかがあると非常に便利だと思います。簡単ではありますが、本日は以上です！

ArcPyレシピ集③ ～サブタイプの取得方法～

さて、ブログの更新が滞っていましたが再開させようと思います。今回はArcPyでサブタイプの取得方法を紹介したいと思います。

サブタイプとは？！

そもそもサブタイプとは何なのでしょうか？

サブタイプは、フィーチャクラスのデータを分類する場合に使用します。例えば、道路を表現する場合に高速道路、国道、地方道と別々のフィーチャクラスを作成するのではなく、道路という単一のフィーチャクラスの属性値で種類を高速道路、国道、地方道のように分類します。（←ESRIさんのサイトから丸パクリしました）

つまり、一つのレイヤ（道路）の中で複数の地物（高速道路や国道など）を管理することができるということですね。例えば、DMなんかだとレイヤ数が200以上になってしまいます。これだとレイヤ数が多くなって管理しにくいので、サブタイプを使用する必要があります。そうすればレイヤ数を減らして管理しやすくなりますね（Zmapなども同様かと思います）。

環境

Windows7 64bit
ArcGIS10.2.0
Python2.7.2

データ

・FGDBの作成(Sutype.gdbという名称にしました)
・フィーチャクラスの作成(ライン。道路という名称にしました)

サブタイプの設定

以下のように設定しました。

サンプルコード

とても簡単ですね。まずは「ListSubtypes」という関数を使用してサブタイプの情報を取得します。その後、ごにょごにょやってサブタイプ名とサブタイプコードを取得しています(←おい、ちゃんと説明しろよという方はメッセージください)。

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

def get_subtype_field(fc):
    """サブタイプをディクショナリで取得します。

    arguments:
        fc -- フィーチャクラス

    returns:
        subDict -- {サブタイプ名：サブタイプコード}

    """

    subDict = {}
    try:
        subtypefields = arcpy.da.ListSubtypes(fc)

        for name in subtypefields.keys():

            group_dict = subtypefields[name]

            for stkey in list(group_dict.keys()):

                if stkey == 'SubtypeField':

                    fields = group_dict[stkey]

                    for stkey in list(group_dict.keys()):

                        if stkey == 'Name':

                            subfields = group_dict[stkey]
                            subDict[subfields] = name

        print "\n".join("%s: %s" % i for i in subDict.items())

    except Exception:

        return None

if __name__ == '__main__':
    arcpy.env.workspace = u"C:\ArcPySample\Sutype.gdb"
    get_subtype_field(u"道路")

結果を見てみるとばっちりですね！

本サンプルに関してですが、例えば、サブタイプごとのアイテム数を抽出したい場合、まずはサブタイプの取得が必要になります。また、多くのフィーチャクラスが存在する中で、どのフィーチャクラスがサブタイプを使用しているかの判定にも使えるかと思います。サブタイプはとても便利なので、まだ使ったことがないという方はぜひ使ってみてください！本日は以上です！！

ArcPyレシピ集② ～ワークスペースのフィーチャクラス名とジオメトリタイプの取得方法～

さて、GWも終わり憂鬱な気分となっておりますが、本日もArcPyレシピを書いてみたいと思います。

本日のお題

• ワークスペースのフィーチャクラス名とジオメトリタイプの取得方法

おそらく業務で使用するワークスペース（FGDBなど）には多数のフィーチャクラスが存在しているのではないでしょうか。状況によってはそれらの名称とジオメトリタイプを取得したい場合があるかと思います。今回はこれを題材にしたいと思います。

対象のワークスペースとフィーチャクラス

以下のようなサンプルデータを作成しました。

①ワークスペース
・WorkSpaceSample.gdb

②フィーチャクラス
・HOUSE・・・ポリゴン
・HOUSE_ANO・・・アノテーション
・ROAD・・・ライン
・SHOP・・・ポイント

これらをディクショナリで取得しようと思います。{'ROAD': 'Polyline', 'SHOP': 'Point'}のような感じですね。

環境

ArcGIS10.2
Python2.7.3

サンプル

今回も簡単ですね。

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

def get_all_featureclass_type():

    geo_type_dict = {}
    fc_list = arcpy.ListFeatureClasses()

    for fc in fc_list:

        geo_type_dict[fc] = get_featureclass_type(fc)

    return geo_type_dict

def get_featureclass_type(fc):

    fc_type = arcpy.Describe(fc).FeatureType

    #アノテーションのジオメトリタイプはポリゴンになってしまうのでこの処理を追加
    if fc_type == "Annotation":

        geo_type = "Annotation"

    else:

        geo_type = arcpy.Describe(fc).shapeType

    return geo_type

arcpy.env.workspace = ""  #gdbを指定
print get_all_featureclass_type()

以下のような結果になりました。

解説

①まず、「ListFeatureClasses()」で「arcpy.env.workspace」で指定したワークスペースに存在するフィーチャクラスのリストを取得します（上記「②フィーチャクラス」参照）。
②取得したフィーチャクラスごとにループして「get_featureclass_type」を呼び出しています。
③「get_featureclass_type」内でジオメトリタイプを取得しています。

※1.アノテーションのジオメトリタイプはポリゴンとして認識されてしまうので、「Annotation」という文字列を代入しています。
※2.「Describe」とは、地理データのプロパティにアクセスする手段です。例えばフィーチャクラスの「Describe」は以下です。

「Describe」はけっこう使える関数なので、ぜひ積極的に使ってみてください！簡単ですが、本日は以上です。

ArcPyレシピ集① ～フィーチャクラスのカラム名の取得方法～

さて、本日から簡単なArcPyのレシピ集を書いてみたいと思います。もっといい方法があるよ！という方はご連絡いただけると助かりますm(_ _)m

本日のお題

• フィーチャクラスのカラム名の取得方法

フィーチャクラスを扱うにおいて、カラム名をとってきたいことはよくあると思います。今回はこれを題材にしたいと思います。

対象フィーチャクラス

こんな感じのフィーチャクラスを作ってみました。

環境

ArcGIS10.2
Python2.7.3

サンプル

超簡単です。フィーチャクラスのカラム名をリストにつっこんでいます

def get_column_name(in_feature):
    """フィーチャクラスのカラム名をリストで返します"""

    fields = arcpy.ListFields(in_feature)
    return [field.name for field in fields]

arcpy.env.workspace = ""      #gdbを指定
columns = get_column_name("") #フィーチャクラスを指定

print columns

以下のような結果になりました。超簡単ですね！

解説

解説するほどでもありませんが、arcpy.ListFields()で対象フィーチャクラスのFieldオブジェクトを取得することができます。Fieldオブジェクトのプロパティは以下のようになっています。

「name」というプロパティがありますね。今回はそれをリスト内包表記で取得しています。本当にこれだけですね。
ListFieldsだけではなく、ListFeatureClasses（対象のワークスペースに存在するフィーチャクラスを取得するメソッド）とか色々あるので、ぜひ使ってみてください。

このレシピ集、どこまで続けられるか（ブログ主の知識不足）。。簡単ですが、本日は以上です！

ArcPyでオープンストリートマップのデータ（.osm)をShapeに変換してみた

さて、本日はArcPyでオープンストリートマップのデータ（.osm)をShapeに変換してみようと思います。最近はまっているオープンストリートマップですが、出力できるデータ形式が.osm(実質xml)というものなので、なんだかなぁという感じです。とりあえずShapeで出力してArcMapでも読み込めるようにしようと思います。

仕様

JA:OSM XML - OpenStreetMap Wiki　に仕様が載っているのですが、大事なところだけをピックアップします。

データ形式

• ノード ブロック。属性として緯度経度を持つ

→各ノードの属性（名称とか）

<node id="1420826646" visible="true" version="1" changeset="9191805" timestamp="2011-09-02T10:28:00Z" user="watao" uid="393606" lat="35.5526160" lon="139.6462920">
  <tag k="amenity" v="police"/>
  <tag k="KSJ2:AdminArea" v="神奈川県横浜市港北区"/>
  <tag k="KSJ2:ADS" v="日吉本町1-23-21"/>
  <tag k="KSJ2:PubFacAdmin" v="都道府県"/>
  <tag k="name" v="港北警察署日吉駅前交番"/>
  <tag k="note" v="National-Land Numerical Information (Public Facility) 2006, MLIT Japan"/>
  <tag k="note:ja" v="国土数値情報（公共施設データ）平成19年　国土交通省"/>
  <tag k="source" v="KSJ2"/>
  <tag k="source_ref" v="http://nlftp.mlit.go.jp/ksj/jpgis/datalist/KsjTmplt-P02-v2_0.html"/>
 </node>

• ウェイ ブロック

→各ウェイごとに、そのノードへの参照（ノードへの参照から緯度経度を取得できるみたい）
→各ウェイの属性（名称とか）

<way id="378405415" visible="true" version="3" changeset="35164763" timestamp="2015-11-08T08:44:40Z" user="Madoka" uid="421690">
  <nd ref="3818193608"/>
  <nd ref="3818193607"/>
  <nd ref="3818193610"/>
  <nd ref="3824079880"/>
  <nd ref="3818193605"/>
  <nd ref="3818193606"/>
  <nd ref="3818193608"/>
  <tag k="amenity" v="cafe"/>
  <tag k="indoor" v="room"/>
  <tag k="level" v="0"/>
  <tag k="name" v="タリーズコーヒー　慶應日吉店"/>
  <tag k="opening_hours" v="Mo-Su 07:30-20:00"/>
  <tag k="room" v="yes"/>
  <tag k="source" v="survey"/>
 </way>

• リレーション ブロック

→各リレーションごとに、そのメンバーへの参照
→各リレーションの属性（名称とか）

<relation id="5647223" visible="true" version="3" changeset="36006131" timestamp="2015-12-17T10:28:22Z" user="nakajyou" uid="2435742">
  <member type="way" ref="377939570" role="part"/>
  <member type="way" ref="377939569" role="part"/>
  <member type="way" ref="159628078" role="part"/>
  <member type="way" ref="377939571" role="part"/>
  <member type="way" ref="377939573" role="part"/>
  <member type="way" ref="377939574" role="part"/>
  <member type="way" ref="377939572" role="part"/>
  <member type="way" ref="377939568" role="part"/>
  <member type="way" ref="377939567" role="part"/>
  <member type="relation" ref="5748542" role=""/>
  <member type="relation" ref="5647285" role=""/>
  <tag k="name" v="協生館"/>
  <tag k="type" v="building"/>
 </relation>

ノード＝ポイント、ウェイ＝ライン、ポリゴンですね。ウェイに関しては始点終点が一致している場合、ポリゴン、そうでない場合、ラインと判断できるっぽいです。リレーションというのは、各アイテムの関連性のことでしょうか（例えば、駅ビルと各テナントの関連性とか）？とりあえず、ここは不要と判断したので、本記事では割愛いたします。また、取得する属性は名称だけにしました（細かい仕様が載ってないので、カラムをどうするか判断がつきませんでした・・・）。

出力するShape

• ポイント
• ライン
• ポリゴン

どれも属性は名称のみです。ポイントのカテゴリ（例えば、スーパーマーケットとか）とかでShapeを分けた方がいいと思いますが、とりあえず今回はまとめて出力してみます。

出力するマップ

東急東横線の日吉駅前
↑とてもいいところです

サンプルコード

■exportosm.py

# -*- coding: utf-8 -*-
import arcpy
import os
from xml.etree import ElementTree

class OSM(object):

    def __init__(self,osm_path,out_dir):
        arcpy.env.workspace = out_dir
        self.osm = osm_path
        self.out_dir = out_dir
        self.file_name_point = "osm_point.shp"
        self.file_name_line = "osm_line.shp"
        self.file_name_polygon = "osm_polygon.shp"
        self.out_file_point = os.path.join(out_dir,self.file_name_point)
        self.out_file_line = os.path.join(out_dir,self.file_name_line)
        self.out_file_polygon = os.path.join(out_dir,self.file_name_polygon)
        self.spatial_reference = arcpy.SpatialReference(4612)
        self.node_dict = {}
        self.way_dict = {}

    def create_shape(self):
        """Shapeファイルを作成します"""

        exist_out_feature = arcpy.ListFeatureClasses(self.file_name_point)

        if len(exist_out_feature) == 0: #指定したディレクトリにosm_point.shpが存在しない場合は作成
            arcpy.CreateFeatureclass_management(self.out_dir,
                                                self.file_name_point,
                                                "POINT",spatial_reference = self.spatial_reference)

        self.add_columns(self.file_name_point) #カラム追加

        exist_out_feature = arcpy.ListFeatureClasses(self.file_name_line)

        if len(exist_out_feature) == 0: #指定したディレクトリにosm_line.shpが存在しない場合は作成
            arcpy.CreateFeatureclass_management(self.out_dir,
                                                self.file_name_line,
                                                "POLYLINE",spatial_reference = self.spatial_reference)

        self.add_columns(self.file_name_line) #カラム追加

        exist_out_feature = arcpy.ListFeatureClasses(self.file_name_polygon)

        if len(exist_out_feature) == 0: #指定したディレクトリにosm_polygon.shpが存在しない場合は作成
            arcpy.CreateFeatureclass_management(self.out_dir,
                                                self.file_name_polygon,
                                                "POLYGON",spatial_reference = self.spatial_reference)

        self.add_columns(self.file_name_polygon) #カラム追加

    def add_columns(self,file_path):
        """Shapeファイルにカラムを追加します"""

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

        if not "name" in field_list: #nameカラムがなかったら作成
            arcpy.AddField_management(file_path, "name", "text")


    def read_osm(self):
        """osmファイルを読み込んで、必要な情報を取得します"""

        tree = ElementTree.parse(self.osm)

        for node in tree.findall(".//node"): #nodeを検索

            for e in list(node):
                if e.attrib["k"] == "name": #名称が存在した場合
                    self.node_dict[node.attrib["id"]] = {"name":e.attrib["v"],"lat":node.attrib["lat"],"lon":node.attrib["lon"]}
                    break

            else: #名称が無い場合はこっち
                self.node_dict[node.attrib["id"]] = {"name":"","lat":node.attrib["lat"],"lon":node.attrib["lon"]}

        tmp_way_dict = {}
        way_list = []
        name = ""

        for way in tree.findall(".//way"): #wayを検索

            for e in list(way):

                if e.attrib.has_key("ref"): #nodeidを取得

                    way_list.append(e.attrib["ref"])

                if e.attrib.has_key("k"):

                    if e.attrib["k"] == "name": #名称がある場合

                        name = e.attrib["v"]

                else:

                    name = ""

            t = set() #重複チェック(wayがlineかpolygonなのか)
            duplicate_list = [x for x in way_list if x in t or t.add(x)]

            if len(duplicate_list) > 0: #重複があったらポリゴンとみなす
                tmp_way_dict[way.attrib["id"]] = {"node_id" :way_list,"name":name,"type":"polygon"}
            else:
                tmp_way_dict[way.attrib["id"]] = {"node_id" :way_list,"name":name,"type":"line"}

            way_list = []

        node_coord = []
        for way in tmp_way_dict.items(): #way_dictに格納したnode_idから実際の座標を取得する

            for c in way[1]["node_id"]:
                node_coord.append([self.node_dict[str(c)]["lat"],self.node_dict[str(c)]["lon"]])

            self.way_dict[way[0]] = {"coord":node_coord,"name":way[1]["name"],"type":way[1]["type"]}

            node_coord = []



    def create_geometry(self):
        """メインメソッド"""

        self.create_shape() #Shape作成
        self.read_osm()     #OSM情報の取得
        self.create_node()  #nodeからポイントを作成
        self.create_way()   #wayからラインとポリゴンを作成

    def create_node(self):
        """nodeからポイントを作成します"""

        array = arcpy.Array()
        with arcpy.da.InsertCursor(self.out_file_point, ["SHAPE@","name"]) as cur:

            for nd in self.node_dict.items():

                point = arcpy.Point()
                point.X = nd[1]["lon"]
                point.Y = nd[1]["lat"]

                pointGeometry = arcpy.PointGeometry(point)

                cur.insertRow((pointGeometry,nd[1]["name"]))

    def create_way(self):
        """wayからラインとポリゴンを作成します"""

        #ラインを作成
        with arcpy.da.InsertCursor(self.out_file_line, ["SHAPE@","name"]) as cur:

            for wd in self.way_dict.items():

                if wd[1]["type"] == "line":

                    array = arcpy.Array()

                    for part in wd[1]["coord"]:

                        array.add(arcpy.Point(part[1],part[0]))

                    cur.insertRow((arcpy.Polyline(array,self.spatial_reference),wd[1]["name"]))


        features = []

        #ポリゴンを作成
        with arcpy.da.InsertCursor(self.out_file_polygon, ["SHAPE@","name"]) as cur:

            for wd in self.way_dict.items():

                if wd[1]["type"] == "polygon":

                    array = arcpy.Array()

                    for part in wd[1]["coord"]:

                        array.add(arcpy.Point(part[1],part[0]))

                    cur.insertRow((arcpy.Polygon(array,self.spatial_reference),wd[1]["name"]))

if __name__ == '__main__':
    osm = OSM("","") #osmのパスとShapeの出力先フォルダを指定
    osm.create_geometry()

結果を確認します。

完全に日吉駅前ですね！素晴らしい！！

位置もぴったりです！ただ、建物の四隅に必ずノードがあることに気づきました。これはオープンストリートマップの仕様なんですかね？！なんか気持ち悪い。あと、プログラム作った後に気づいたのですが、オープンストリートマップを使用する方はみなさんQGISとかを使っているのでしょうか？！（商用のGISには頼らないという思想？！）だとすると、このスクリプトはあまり役に立たないのかなぁと・・・一応、「read_osm」メソッドのところを流用すれば緯度経度や名称を取得することはできますが。。。

とりあえず今回は簡易的な出力を行ってみたのですが、osmファイルの仕様をきちんと理解すればもっと細かい精度でShapeに出すことができそうですね。

簡単ですが、今回は以上です！

OpenStreetMapを使ってみよう！

さて本日はマニアックなArcPyのプログラムからは離れてOpenStreetMapを紹介しようと思います。実は私も最近知りました。GISに携わる者としてこんな面白いものを知らなかったのは実に恥ずかしい。。。

OpenStreetMapとは

「OpenStreetMap(OSM)は、道路地図などの地理情報データを誰でも利用できるよう、フリーの地理情報データを作成することを目的としたプロジェクトです。
誰でも自由に参加して、誰でも自由に編集でき、誰でも自由に利用する事が出来ます。」

要するにウィキペディアの地図版ってとこですかね？！以下サイトから利用することができます。

https://openstreetmap.jp/

地図の見方

以下ボタンをクリックしてください。

地図が表示されます。なかなか綺麗ですね！

拡大するとこんな感じです。ここが拡大の限界なのが微妙ですが、なかなかいいですね！

地図の編集の仕方

以下ボタンをクリックしてください。

編集したい場所に移動します。

ポイント、ライン、エリアのいずれかを選択します。

ポイントを置きます。

必要な情報を書き込みます。

保存します。

Saveして完了です。

これで自分の作成した情報がOpenStreetMapに反映されます。作図だけではなく、既存の情報を修正することもできます（コンビニが閉店してラーメン屋になった場合、ポイントの情報を変更する）。

地図データの出力

エクスポートを押すとデータをエクスポートできます。ただ、形式が「.osm」という形式のため、ArcGISなどで読み込むことができません（なんかプラグイン的なのが公開されていましたが、ArcGIS10.2.0用のはありませんでした・・・）。こいつをShapeに変換できたらすごく便利なのですが、それは今後調べてみようと思います。

地図を閲覧するだけならGoogleMapの方がいいと思いますが、この取組自体がとても素晴らしいと感じますし、色々メリットも多いと思います（例えば、通常地図には載らないであろう情報を自分や自分の仲間むけに載せることができますね。私が勤める会社の業務的にも使えそうな気がします）。

サイトの使い方などを少し勉強する必要がありますが、ぜひ皆さんも活用してみてはいかがでしょうか！