geopandasを使ってみよう！

さて、今回はgeopandasですね！以前から気になっていたライブラリなので、これを機に使ってみようと思います。

geopandasとは

geopandasはpandasの拡張で、地理データを含むデータをpandasのように表形式で扱うことができるpythonのライブラリです。こんな便利なものを作っていただけるとはなんともありがたい話ですね。

ちなみにこちらのエントリーでより詳しくgeopandasについて紹介していますので、興味のある方はぜひ読んでみてください。
www.gis-py.com

インストール

公式サイトによるとhttp://geopandas.org/install.html、GeoPandasをインストールする前に以下ライブラリをインストールする必要があるそうです。その後、pip install geopandas ですね！

• numpy
• pandas (version 0.15.2 or later)
• shapely
• fiona
• six
• pyproj
• descartes

あとは必須ではないですが、オプションとして以下も入れておくといいと思います。

• matplotlib
• geopy

もしpipでエラーが出るようなことがありましたらこのサイトhttps://www.lfd.uci.edu/~gohlke/pythonlibs/からWheelをダウンロードしてインストールしてみてください。

実行環境

Windows7 64bit
ArcGIS10.4.1
Python2.7.3

データ

geopandasをインストールすると自動的にサンプルデータもついてくるみたいです。今回はそれを使用してみます。データの内容は以下です。

1. 世界地図
2. 世界の都市

また、今回はjupyterを使用してデータを展開してみようと思います。

■データの取得

import geopandas

world = geopandas.datasets.get_path('naturalearth_lowres')      #世界地図
cities = geopandas.datasets.get_path('naturalearth_cities')         #都市

これだけでＯＫです。簡単ですね。

■データの読み込み

df_world = geopandas.read_file(world)
df_cities = geopandas.read_file(cities)

こちらも簡単。取得したデータをデータフレーム化するだけですね。

片方のデータフレームの中身を確認します。

df_world.head()

ちゃんとデータが入っていることが確認できました。

■データのプロット
公式ドキュメントによるとdataflameに対して.plot()を使用することで、matplotlibが走るようです。

base = df_world.plot(color='white', edgecolor='black')
df_cities.plot(ax=base, marker='o', color='red', markersize=5);

おぉー素晴らしい。サンプルデータが描画されました。

サンプルコード

上記をまとめたものを記載します。

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

world = geopandas.datasets.get_path('naturalearth_lowres')      #世界地図
cities = geopandas.datasets.get_path('naturalearth_cities')         #都市

df_world = geopandas.read_file(world)
df_cities = geopandas.read_file(cities)

base = df_world.plot(color='white', edgecolor='black')
df_cities.plot(ax=base, marker='o', color='red', markersize=5);

このgeopandasですが、かなり面白いライブラリですね。今回は基礎的なことしか触れなかったのですが、geopandasでShapeファイルを作成するエントリーもありますので、よかったら読んでみてください。
www.gis-py.com

年々GISデータの重要性が高まる中、こういった便利なライブラリがもっと登場してくれるといいですね。

以上、「geopandasを使ってみよう！」でした。

【PythonとGIS】GDALを使ってShapeファイルを作成してみる

皆さん、お久りぶりです。ずっとブログ更新をさぼってしまっていましたが、またちょくちょく書いていこうと思います。さて、今回はGDALですね！GDALは以下記事のように当ブログで紹介はしているのですが、今回はShapeファイルの作成にチャレンジしてみたいと思います。

www.gis-py.com

www.gis-py.com

sanvarie.hatenablog.com

sanvarie.hatenablog.com

また、今までのブログではArcPyでオープンストリートマップのデータ（.osm)をShapeに変換してみた - GIS奮闘記のようにArcpyを使用してのShapeファイルの作成方法は紹介したことがあるのですが、これだとArcGISに依存した作成方法なので、ArcGISを持っていない方でもShapeファイルの作成が手軽にできるようにGDALのようなフリーのライブラリでの作成方法を紹介します。

GDALとは

GDAL - Wikipedia
色々小難しいことが書いてありますが、要はフリーのGISライブラリです。

インストール

以下のようなサイトを参考にインストールをお願いします。
http://blog.godo-tys.jp/2013/06/19/2465/
https://github.com/domlysz/BlenderGIS/wiki/How-to-install-GDAL

実行環境

Windows7 64bit
ArcGIS10.4.1
ArcMap10.4.1
Python2.7.3

サンプルコード

ポイントShapeを作成するサンプルです。

# -*- coding: utf-8 -*-
import osgeo.ogr as ogr
import osgeo.osr as osr

# 今回作成するポイントの緯度経度
lat_lon = [[35.679933,139.714465],[35.1231433,139.2533465]]

# shapefileドライバ
driver = ogr.GetDriverByName("ESRI Shapefile")

# Shape出力先を設定
data_source = driver.CreateDataSource(ur"D:\gis-py\gdal\sample.shp")

# 投影法を設定（projectionファイルがうまく作成されないので、コメントアウト）
#srs = osr.SpatialReference()
#srs.ImportFromEPSG(4612)

# レイヤ作成
#layer = data_source.CreateLayer("sample", srs, ogr.wkbPoint) # projetionファイル作るときに使う
layer = data_source.CreateLayer("sample",  geom_type= ogr.wkbPoint)

# フィールド追加
field_name = ogr.FieldDefn("column", ogr.OFTString)
field_name.SetWidth(24)
layer.CreateField(field_name)

layer.CreateField(ogr.FieldDefn("Latitude", ogr.OFTReal))
layer.CreateField(ogr.FieldDefn("Longitude", ogr.OFTReal))

# Shape作成
for l in lat_lon:
  # フィーチャ作成
  feature = ogr.Feature(layer.GetLayerDefn())

  # 属性を設定
  feature.SetField("column", "test")
  feature.SetField("Latitude", l[0])
  feature.SetField("Longitude", l[1])

  # 作成するポイントに緯度経度をセット
  wkt = "POINT(%f %f)" %  (float(l[1]) , float(l[0]))

  # ポイント作成
  point = ogr.CreateGeometryFromWkt(wkt)
  feature.SetGeometry(point)
  layer.CreateFeature(feature)
  feature = None

data_source = None

おー、ちゃんとShapeファイルができあがっていますね！

しかし、projectionファイルがありません。ファイルの作成はできたのですが、なぜか中身が空！解決策がわからず、今回は泣く泣く作成せずに終わりました。

ArcMapでShapeを読み込んでみました。ポイントが作成されていることがわかります。

属性もしっかり入っていますね。

GDALのAPIについてはあまり日本語化されておらず、情報量もあまり多くないのですが、使いこなすことができたらとても便利なライブラリだと思います。これを機にぜひ使ってみてください！

以上、「GDALを使ってShapeファイルを作成してみる」でした！

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（対象のワークスペースに存在するフィーチャクラスを取得するメソッド）とか色々あるので、ぜひ使ってみてください。

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