GIS奮闘記

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

広告

GISエンジンおよびサービスのまとめ

今日は私が日常的に使用しているGISエンジンおよびサービスについてまとめてみようと思います。

GISって何?という方は以下のエントリーを読んでみてください。
www.gis-py.com

1.ArcGIS

言わずと知れたGIS界の巨人ESRIが展開するGISエンジンです。様々な製品を展開しており、正直、私もすべてを把握できてはいません。

ArcGIS製品で私がよく使用するものは以下です。
①ArcGIS Desktop
GISデータを参照、編集することができるソフトです。ArcGIS Desktop | ESRIジャパン

f:id:sanvarie:20181016102438p:plain

おそらくこれがあれば、ほとんどのGIS関連の業務は事足りるのではないでしょうか。分析や色分け機能などがとても充実しており、綺麗な地図を作成することができます。


②ArcGIS Engine
ArcGIS Desktopで大抵のことはできてしまうのですが、例えば、特定の業務を行う際、ArcGIS Desktopだと少し物足りない、もしくは、作業手順が煩雑になるという場合、必要なGISアプリを開発することもあるかと思います。その際に使用するのがArcGIS Engineです。ざっくりいうとライブラリですね。これを使用すれば痒い所に手が届くアプリを開発することができます。
ArcGIS Engine | ESRIジャパン


③ArcPy
製品というわけではないと思うのですが、ArcGISをインストールすると一緒にインストールされるのがこのArcPyです。作業の自動化はしたいが、ArcGIS Engineを使うほどではないという場合に非常におすすめです。名前から推測できるかと思いますが、Pythonを使ってArcGISデータを扱うことができます。ArcGIS Desktop上で使用することができるので、とても便利です。

f:id:sanvarie:20181016103731p:plain

ArcPyについては色々なエントリーを書いていますので、ぜひ読んでみてください。
www.gis-py.com


※使ってみたいサービス
④ArcGIS Online
近年、ESRIが力を入れているクラウド GISです。ユーザーがいつでもどこでも手軽に利用できるというメリットがあります。ただ、使用料が従量課金制というのがいまいちかなと個人的には思っています(ただし、参照のみの場合は無償)。
ArcGIS Online | ESRIジャパン

2.SIS

ArcGISほどメジャーなGISエンジンではありませんが、こちらも私はよく使用します。Cadcorpという会社が販売しているGISエンジンです。https://www.cadcorp.com/

こちらの会社もESRIほどではないですが、色々な製品を展開しています。

SIS製品で私がよく使用するものは以下です。
①SIS Map Modeller
GISデータを参照、編集することができるソフトです。
f:id:sanvarie:20181016104339p:plain

この製品の特長としては色々なファイルを開くことができるということでしょうか。この点に関してはArcGIS Desktopに勝っているのではないかと思います。

あと、ArcGIS Desktopは大容量のShapeファイルを開こうとしたりすると、頻繁に落ちたりするのですが、SIS Map Modellerではまずそういったことはありません。また、ArcGISのデータ形式であるFGDBはよく壊れたりするのですが、SISのデータ形式であるbdsは壊れることはあまりありません。こういったところは本当にSISは非常によくできており、もう少しメジャーになってもいいのではないかと思うエンジンですね。

②SIS SDK
ArcGIS Engineと同じように、SISをベースとしたアプリケーションを開発する際に使用します。

3.QGIS

ArcGISやSISと異なり、こちらは無償のGISエンジンです。無償といってもその機能はとても充実しています。参照、編集に関しては有償のエンジンとの差はほとんどないかもしれません。ただ、細かい分析などをしたい場合はArcGISに軍配が上がるかと思います。個人でGISエンジンを利用する場合はQGISを選ぶ方が大多数ではないでしょうか。

①QGIS Desktop
GISデータを参照、編集することができるソフトです。ArcGIS DesktopやSIS Map Modellerと同じですね。

f:id:sanvarie:20181016110753p:plain

4.Google Map

こちらは言うまでもないでしょうか。世界ナンバーワンのGISサービスですね。
Google マップ

f:id:sanvarie:20181016113000p:plain


Google Maps APIを使用すればGoogle Mapを使用してGISアプリを作ることができます。Google Maps APIに興味がある方は以下を読んでみてください。

www.gis-py.com

5.Google Earth

Google Mapほどではないにしろ、こちらも人気のGISサービスですね。Google Earth

f:id:sanvarie:20181016115313p:plain

うーん、かっこいい。

6.地理院地図

国土地理院が展開するGISサービスですね。少しマニアックになってしまいますが、GIS関係の仕事をしている方はよく使われるのではないでしょうか。
地理院地図

Google Mapとは違った見た目ですが、渋くて自分は大好きですね。
f:id:sanvarie:20181016115609p:plain

7.OpenStreetMap

OpenStreetMap(OSM)は、道路地図などの地理情報データを誰でも利用できるよう、 フリーの地理情報データを作成することを目的としたプロジェクトです。地図版のWikipediaといったところでしょうか。

OpenStreetMap

f:id:sanvarie:20181016122753p:plain

OpenStreetMapに興味がある方はぜひこちらを読んでみてください。
www.gis-py.com


個人的な評価はざっとこんな感じでしょうか。やはり機能面は有償のGISエンジンが優れていると思いますが、価格や手軽さを考慮すると無償のGISの存在というのはとても大きいと思います。ただ、どのエンジンやサービスを使うべきというのは無くて、各々の目的に適したものを使用されるのがいいかと思います。

GIS 無償 機能 手軽さ
ArcGIS × ×
SIS × ×
QGIS ×
Google Map
Google Earth
地理院地図
OSM

geopandasを使ってShapeファイルを作成しよう! ~Airbnbのデータを可視化してみよう~

さて、今回はgeopandas第二弾です。前回はgeopandasの簡単な使い方を紹介しましたが、今回はもう少し掘り下げてShapeファイルの作成をしてみようと思います。

geopandas第一弾に関しては以下をご参照ください。
www.gis-py.com

Airbnbとは?

Airbnb(エアビーアンドビー)は、宿泊施設・民宿を貸し出す人向けのウェブサイトである。世界192カ国の33,000の都市で80万以上の宿を提供している。そうです。私は興味はあるのですが、一度も利用したことがありません・・・いつか使ってみたいですね!
Airbnb - Wikipedia

データ

以下のサイトからAirbnbのオープンデータをダウンロードすることができます。しかし、日本のデータがない!日本ではあまり人気がないからでしょうか?
Get the Data - Inside Airbnb. Adding data to the debate.

なので、今回はブログ主に馴染み深い西オーストラリアのデータを使用したいと思います。

赤枠が今回使用したデータです。上の赤枠が宿泊施設データを格納したCSV、下の赤枠が地形データを格納したgeojsonですね。
f:id:sanvarie:20181012105442p:plain

実行環境

Windows7 64bit
ArcGIS Desktop10.4.1
Python2.7.3

サンプルコード

ダウンロードしたCSVとgeojsonからShapeファイルを作成するサンプルです。

# -*- coding: utf-8 -*-
import pandas as pd
import geopandas
from shapely.geometry import Point

#csvの読み込み
df = pd.read_csv(r"D:\gis-py\geopandas\airbnb\listings.csv")

#ポイント作成
geometry = [Point(xy) for xy in zip(df.longitude, df.latitude)]
geo_df = geopandas.GeoDataFrame(df, geometry=geometry)

#Shape出力
geo_df.to_file(driver='ESRI Shapefile', filename=r"D:\gis-py\geopandas\airbnb\accommodation.shp")

#geojsonの読み込み
df2 = geopandas.read_file(r"D:\gis-py\geopandas\airbnb\neighbourhoods.geojson")

#Shape出力
df2.to_file(driver='ESRI Shapefile', filename=r"D:\gis-py\geopandas\airbnb\western_Australia.shp")

ものすごいあっさり作れてしまいました。本当に便利なライブラリですね。ArcMapを使用して作成したShapeを確認してみます。

f:id:sanvarie:20181012110339p:plain

おぉーまさしく西オーストラリアですね!そして、宿泊施設は都市部に集中していることがわかりますね。たしかオーストラリアの真ん中らへんは砂漠地帯だったと思うので、宿泊施設がほとんどないことも納得がいきます。

属性はこんな感じです。各宿泊施設のWEBサイトに載っているような情報(値段とか場所とか)が入っているみたいですね。

f:id:sanvarie:20181012111238p:plain

ちなみに西オーストラリア最大の都市であるパースの宿泊施設が何個あるか空間検索して調べてみました。

f:id:sanvarie:20181012111644p:plain

なんと710個!やはり都市部は違いますね。

今回の記事ではgeopandasを使用すれば、本当に簡単にShapeを作成することができるということがわかりました。ただ、もっと色々な機能を有しているライブラリですので、応用すればもっと色々なことができるかと思います。これを機にぜひ使ってみてはいかがでしょうか。

以上、今回は「geopandasを使ってShapeファイルを作成しよう! ~Airbnbのデータを可視化してみよう~」でした。

geopandasを使ってみよう!

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

geopandasとは

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

インストール

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

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

もし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')         #都市

これだけでOKです。簡単ですね。

■データの読み込み

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

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

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

df_world.head()

f:id:sanvarie:20181009105523p:plain


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


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

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

f:id:sanvarie:20181009111006p:plain


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

サンプルコード

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

# -*- 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ですが、かなり面白いライブラリですね。今回は基礎的なことしか触れなかったのですが、今後、もっと深いところまで使用した内容を紹介できればと思います。
また、年々GISデータの重要性が高まる中、こういった便利なライブラリがもっと登場してくれるといいですね。

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

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

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


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ファイルができあがっていますね!
f:id:sanvarie:20181005130142p:plain


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

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


属性もしっかり入っていますね。
f:id:sanvarie:20181005130459p:plain


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フィーチャークラスを作成します。カラムは以下のように作成します。

f:id:sanvarie:20170108153149p:plain

f:id:sanvarie:20170108153155p:plain

f:id:sanvarie:20170108182218p:plain

このような感じになるかと思います。
f:id:sanvarie:20170108160203p:plain

株式市場の種類について

参考までに

日本の株式市場は、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()

このようにポイントがプロットされたら成功です!
f:id:sanvarie:20170110083034p:plain

属性もばっちり入っていますね。
f:id:sanvarie:20170110083105p:plain

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

f:id:sanvarie:20170110083351p:plain

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

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

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

環境

Windows7 64bit
ArcGIS10.2.0
Python2.7.2

データ

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

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

f:id:sanvarie:20160719082703p:plain

f:id:sanvarie:20160719083300p:plain

サンプルコード

対象の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

f:id:sanvarie:20160719090517p:plain

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

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

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

サブタイプとは?!

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

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

f:id:sanvarie:20160609083037p:plain

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

環境

Windows7 64bit
ArcGIS10.2.0
Python2.7.2

データ

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

f:id:sanvarie:20160609083919p:plain

f:id:sanvarie:20160609083924p:plain


f:id:sanvarie:20160609083927p:plain

サブタイプの設定

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

f:id:sanvarie:20160609084058p:plain

サンプルコード

とても簡単ですね。まずは「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"道路")

結果を見てみるとばっちりですね!
f:id:sanvarie:20160609090343p:plain

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