GIS奮闘記

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

スポンサーリンク

【Python と GIS】GDAL を使って ファイルジオデータベースを操作してみよう!

さて、本日も GDAL について書いてみようと思います。前回はGDAL で Shapeファイルの読み込みと属性検索について書きましたが、今回はGDAL でのファイルジオデータベースの操作について書いてみようと思います。興味のある方は前回のエントリーもぜひ読んでみてください。

www.gis-py.com

GDALとは

GDAL はフリーの GIS ライブラリです。GDAL — GDAL documentation
これがあればShapeの読み込みからジオメトリの作成まで基本的なことは何でもできるかと思います。

ファイルジオデータベースとは

以下エントリーでも紹介していますが、ファイルジオデータベースとは ArcGIS で使用されるデータフォーマットです。

www.gis-py.com

インストール

公式サイトの手順に従ってインストールをお願いします。https://gdal.org/download.html

実行環境

Windows10(64bit)
Python3.6.5
Jupyter Notebook

サンプルデータ

ArcMap で以下を作成しました。

  1. GDB・・・Sample.gdb
  2. フィーチャクラス・・・GDAL_TEST
  3. カラム・・・テスト1、テスト2

f:id:sanvarie:20190713163137p:plain

f:id:sanvarie:20190713163224p:plain

サンプルコード

GDAL_TEST フィーチャクラスのテスト1、テスト2カラムの値を取得するサンプルです。

from osgeo import ogr

# GDBオープン
driver = ogr.GetDriverByName("OpenFileGDB")
gdb = driver.Open(r"D:\python\data\Sample.gdb", 0)

list1, list2 = [], []

# GDBに格納されているフィーチャクラス数分ループ
for featsClass_idx in range(gdb.GetLayerCount()):

    # フィーチャクラスを取得
    featureClass = gdb.GetLayerByIndex(featsClass_idx)

    layerDefinition = featureClass.GetLayerDefn()

    # GDAL_TESTフィーチャクラスの場合のみ処理を実行
    if featureClass.GetName() == "GDAL_TEST":

        while featureClass:

            feat = featureClass.GetNextFeature()

            if feat is None:
                break

            else:
                for column_num in range(layerDefinition.GetFieldCount()):

                    # カラム名取得
                    column_name = layerDefinition.GetFieldDefn(column_num).GetName()

                    #テスト1、テスト2カラムのときのみ属性を取得する
                    if column_name == "テスト1":
                        list1.append(feat.GetField(column_name))
                    elif column_name == "テスト2":
                        list2.append(feat.GetField(column_name))

print(list1)
print(list2)

del gdb


結果を見ると、想定通りテスト1、テスト2カラムの値を取得できたことがわかりました。
f:id:sanvarie:20190713163525p:plain

処理のポイント

以下にポイントとなる箇所を記述します。

GDBオープン
# GDBオープン
driver = ogr.GetDriverByName("OpenFileGDB")
gdb = driver.Open(r"D:\python\data\Sample.gdb", 0)
フィーチャクラス取得

GetLayerByIndex(index) を使用することでフィーチャクラスを取得することができます。

# GDBに格納されているフィーチャクラス数分ループ
for featsClass_idx in range(gdb.GetLayerCount()):

    # フィーチャクラスを取得
    featureClass = gdb.GetLayerByIndex(featsClass_idx)
フィーチャクラスの走査

GetNextFeature() を使用することで、フィーチャクラス内のフィーチャにアクセスできます。

while featureClass:
    feat = featureClass.GetNextFeature()
属性の取得

基本的には前回のエントリー(https://www.gis-py.com/entry/gdal-shape)と同じですが、今回はリスト内包表記は使用していません。ぱっと見はこちらの方がわかりやすいかもしれないですね。

for column_num in range(layerDefinition.GetFieldCount()):

    # カラム名取得
    column_name = layerDefinition.GetFieldDefn(column_num).GetName()

    #テスト1、テスト2カラムのときのみ属性を取得する
    if column_name == "テスト1":
        list1.append(feat.GetField(column_name))
    elif column_name == "テスト2":
        list2.append(feat.GetField(column_name))


本日は以上になります。また GDAL について書いてみようと思いますが、今度はもう少し突っ込んだ内容にしてみようかと思います。

【Python と GIS】 GDAL を使って Shapeファイルを読み込んで属性検索してみよう!

さて、本日は GDAL について書いてみようと思います。GDAL については以下のように以前に何個か記事を書いているのですが、紹介していない機能がまだまだたくさん残っています。今回は属性検索について紹介します。

www.gis-py.com

www.gis-py.com

www.gis-py.com

www.gis-py.com

GDALとは

GDAL はフリーの GIS ライブラリです。GDAL — GDAL documentation
これがあればShapeの読み込みからジオメトリの作成まで基本的なことは何でもできるかと思います。

インストール

公式サイトの手順に従ってインストールをお願いします。https://gdal.org/download.html

サンプルデータ

以下エントリーで紹介した全国市区町村界データを使用します。

www.gis-py.com

サンプルデータの属性情報

フィールド名

  • JCODE・・・市区町村コード
  • KEN・・・都道府県名
  • SICHO ・・・支庁名・振興局名
  • GUN・・・郡名(町村部のみ)
  • SEIREI・・・政令指定都市の市名
  • SIKUCHOSON・・・市区町村名
  • CITY_ENG・・・市区町村名(英語)
  • P_NUM・・・人口
  • H_NUM・・・世帯数

実行環境

Windows10(64bit)
Python3.6.5
Jupyter Notebook

GDAL で属性検索をするには

以下のような手順を踏む必要があります。

  1. Shape ファイルの読み込み
  2. レイヤの取得
  3. フィルター設定
  4. 属性の取得

サンプルコード

KEN = '神奈川県'のフィーチャを属性検索するサンプルです。

from osgeo import ogr

daShapefile = r"D:\data\japan_ver81.shp"
dataSource = ogr.Open(daShapefile)
layer = dataSource.GetLayer(0)
layerDefinition = layer.GetLayerDefn()
layer.SetAttributeFilter("KEN = '神奈川県'")

for feature in layer:
    attribute_list = [feature.GetField(layerDefinition.GetFieldDefn(column).GetName()) for column in range(layerDefinition.GetFieldCount())]
    print(attribute_list)


結果を見ると、想定通り神奈川県のデータを取得できたことがわかりました。
f:id:sanvarie:20190713143258p:plain

処理のポイント

サンプルコードの処理を重要なポイントごとに分けると以下のようになります。

Shapeファイルの読み込み
daShapefile = r"D:\data\japan_ver81.shp"
dataSource = ogr.Open(daShapefile)
レイヤーの取得
layer = dataSource.GetLayer(0)
フィルター設定
layer.SetAttributeFilter("KEN = '神奈川県'")
属性の取得
attribute_list = [feature.GetField(layerDefinition.GetFieldDefn(column).GetName()) for column in range(layerDefinition.GetFieldCount())]


ここが少し複雑なので、この部分の処理を分解して解説します。

  • レイヤーのフィールド数の取得
for column in range(layerDefinition.GetFieldCount())
  • フィールド名の取得

上記の「レイヤーのフィールド数の取得」で取得した要素を引数にします。

layerDefinition.GetFieldDefn(column).GetName()
  • 属性の取得
feature.GetField("フィールド名")

でそのフィールドの属性を取得することができます。

なので、上記の「フィールド名の取得」で取得した要素、つまり、フィールド名を引数にしています。

feature.GetField(layerDefinition.GetFieldDefn(column).GetName())

少し長いのですが、上記の処理をリスト内包表記にしたものをサンプルコードに載せました。

本日は以上になりますが、まだまだ紹介しきれていない使い方がたくさんあるので、 また GDAL について書いてみようと思います。

Pythonで位置参照情報ダウンロードサービスの Web API を使ってみよう

本日はみんな大好きな位置参照情報ダウンロードサービスについて書いてみようと思います。

位置参照情報ダウンロードサービスとは

以下エントリーでも紹介していますが、日本全国の大字町丁目または街区単位の座標を CSV形式でダウンロードできるWEBサイトです。

www.gis-py.com

位置参照情報ダウンロードサービス WEB APIとは?

このAPIを利用することで、位置参照情報を組み込んだ Webサイトやアプリケーションの開発を行うことが可能になります。また、データをダウンロードするには以下のような画面で一つ一つデータを手動でダウンロードする必要がありますが、APIを使うことでデータを一括でダウンロードしたりすることもできます。

f:id:sanvarie:20190706230310p:plain

API 仕様

API 仕様はここに公開されています。http://nlftp.mlit.go.jp/isj/specification_api_isj.pdf

以下のような形でリクエストを投げるとXML形式で結果が返ってきます。
f:id:sanvarie:20190706231349p:plain

例えば、このようなリクエストを投げると
http://nlftp.mlit.go.jp/isj/api/1.0b/index.php/app/getISJURL.xml?appId=isjapibeta1&areaCode=13000&fiscalyear=%27%E5%B9%B3%E6%88%9025%E5%B9%B4%27&posLevel=0

このような結果が返ってきます。
f:id:sanvarie:20190706231554p:plain

パラメータ

リクエストするURLに入れるパラメータはこのようになっています。必須は以下の二つのみですね。

  1. appId
  2. areaCode

f:id:sanvarie:20190706232101p:plain

areaCode コードですが、http://www.tt.rim.or.jp/~ishato/tiri/code/kanto.htmのようなサイトで一覧を手に入れることができます。

実行環境

Windows(64bit)
Python3.6.5

インストールが必要なライブラリ

lxml を使っていますので、pip install lxml でインストールしてください。

サンプルコード

今回は以下のデータをダウンロード、解凍する処理を書いてみました。

  1. 横須賀市のデータを使用
  2. 街区レベルのデータを使用
  3. すべての年度のデータをダウンロード
  4. 各年度毎にフォルダを作ってそこにデータ(zip)を格納して解凍
# -*- coding: utf-8 -*-
import zipfile
import os.path
import requests
from lxml import etree

class WebAPI(object):

    def __init__(self):
        self.appId = "isjapibeta1" #固定
        self.areaCode = "14201"    #固定
        self.posLevel = "0"        #街区レベル
        self.directory = "D:\data" #フォルダ作成先
        self.new_directory = ""    #ダウンロードしたデータの解凍先

    def get_url(self):
        """
        位置参照情報のWEBAPIを使用して各情報をXMLで取得します。
        """
        url = "http://nlftp.mlit.go.jp/isj/api/1.0b/index.php/app/getISJURL.xml" \
              "?appId=%s&areaCode=%s&posLevel=%s" %(self.appId, self.areaCode, self.posLevel)

        resp = requests.get(url, timeout=10)
        tree = etree.fromstring(resp.content)
        tree = tree.xpath('/ISJ_URL_INF/ISJ_URL/item/child::node()')

        for t in tree:
            if t.tag == "fiscalyear":
                self.make_directory(t.text)

            if t.tag == "zipFileUrl":
                self.download_zip(t.text)

    def make_directory(self,fiscalyear):
        """
        新しいフォルダを作ります
        """
        self.new_directory = os.path.join(self.directory, fiscalyear)
        if os.path.exists(self.new_directory) == 0: #なかったら作る
            os.makedirs(self.new_directory)

    def download_zip(self,text):
        """
        指定したURLからZIPをダウンロードします。
        """

        filename = text.split('/')[-1]
        r = requests.get(text, stream=True)
        with open(filename, 'wb') as f:
            for chunk in r.iter_content(chunk_size=1024):
                if chunk:
                    f.write(chunk)
                    f.flush()

            self.uncompress_zip(filename)

    def uncompress_zip(self,filename):
        """
        ZIP ファイルを指定したディレクトリに展開します。
        """
        zfile = zipfile.ZipFile(filename)
        zfile.extractall(self.new_directory)

if __name__ == '__main__':
    w = WebAPI()
    w.get_url()

このようにファイルがダウンロードされました。
f:id:sanvarie:20190707111150p:plain

各フォルダの中身を確認するとデータが解凍されていることが確認できました
f:id:sanvarie:20190707111246p:plain

大字・町丁目レベルでデータをダウンロードするには

上記プログラムの self.posLevel = "0" をself.posLevel = "1" に変更すればダウンロード可能です。

都道府県ごとのデータをダウンロードするには

self.areaCode = "14201" の部分を各都道府県ごとのコードに変更すればダウンロードできます。APIの使用を見る限りhttp://www.tt.rim.or.jp/~ishato/tiri/code/kanto.htmの中の市区町村コードの上二桁+000 が各都道府県コードになっているかと思います。

感想

けっこうあっさりとデータのダウンロードができました。位置参照情報をよく使うという方はぜひこのAPIを使ってみてください。

Windows環境のPython3でMeCabを使ってみよう

さて、本日はPython3にMeCabを入れてみようと思います。MeCabについては以下のエントリーで紹介しているのですが、Python2系だったのと、あまり詳しくは紹介できませんでした。なので、今回はMeCabをメインにエントリーを書いてみました。

www.gis-py.com

MeCab とは

オープンソースの形態素解析エンジンです。開発者の方が和布蕪(めかぶ)が好物なのでこの名前になったみたいです。公式サイトもありますので、興味のある方はこちらも参考にしてみてください。MeCab: Yet Another Part-of-Speech and Morphological Analyzer

形態素解析 とは

検索エンジンなどにも使用されている自然言語処理の手法の一つで、ある文章・フレーズを「意味を持つ最小限の単位(=単語)」に分解し、それらの品詞等を判別する作業のことです。例えば、「私は横浜に住んでいます。」という文章を形態素解析した結果が以下です。

f:id:sanvarie:20190706135256p:plain

この技術を利用すれば様々なデータを分析することができます。例えば、Twitter などでしょうか。

  1. 特定のハッシュタグを含むツイートを抽出
  2. 「1」のツイートを形態素解析して、その中から地域名を含むツイートを抽出
  3. そこからジオコーディングでツイートされた地域を可視化

こんな感じで使えばトレンド分析みたいなことができるのではないでしょうか?それ以外にも色々な分野で使えるかと思います。

実行環境

Windows10(64bit)
Python3.6.5

手順

色々やることはあるのですが、大きく分けて二つです。

  1. Windows 環境で MeCab のインストール
  2. MeCabをPython上で使えるようにする

※MeCab をインストールしただけではPython上からは実行できませんので、ご注意ください。

Windows 環境で MeCab のインストール

選択肢は以下の二つです。

  1. 公式サイトからMeCabの実行プログラム本体をダウンロード(32bit版)http://taku910.github.io/mecab/#download
  2. 有志がビルドした64bit版をダウンロードhttps://github.com/ikegami-yukino/mecab/releases/tag/v0.996

私の環境のPythonが64bitのため、今回は64bit版を使用します。mecab-0.996-64.exeをダウンロードして、MeCabをインストールします。

環境変数の設定

Pathに「(MeCabのインストール先)/bin(例えば、C:\Program Files\MeCab\bin)」を追加します。

MeCab を使ってみる

コマンドプロンプトでmecabと入力すればMeCab が起動します(環境変数を設定していない方は「(MeCabのインストール先)/bin」のディレクトリに移動してmecabを入力してください)。そして、「私は東京に住んでいます。」という文章を入力してみました。しかし、このように文字化けしてしまいました。
f:id:sanvarie:20190706101145p:plain

ただ、これはコマンドプロンプトがutf-8に対応していないからで特に問題はありません(MeCab 自体は問題なく動いています)。もしコマンドプロンプトをutf-8対応したいという方はコマンドプロンプト起動時、自動的に文字コードをUTF-8にして日本語もちゃんと表示できるようにする方法 - Qiitaのようなサイトを参考にしてみてください。

MeCabをPython上で使えるようにする

MeCabのインストールができたら、それをPython上で動かせるようにする必要があります。以下二つの手順が必要です。

  1. Python上のMeCabバインディングの導入
  2. libmecab.dll をコピペ

まずは、Python上のMeCabバインディングの導入ですね。今回はWindows用のmecab-python-windowsというバインディングを導入します。以下二つのモジュールをpip でインストールします。

pip install ipykernel
pip install mecab-python-windows

※コメントでご指摘いただいたのですが、Microsoft Visual C++ 14.0 が入っていない場合、ここでエラーがでるようです。もしエラーが出た場合は Microsoft Visual C++ 14.0 が入っているかどうか確認してみてください。

インストール完了後、にもう一つ作業(libmecab.dll をコピペ)が必要です。

(MeCabのインストール先)/bin

にある"libmecab.dll"というファイルを

(Pyhthonのインストール先)/Lib/site-packages

にコピペしてください。これでPython上でMecab が使えます。

import MeCab

もしインポートできないという場合は以下を確認してください。

  • ImportError: DLL load failed: %1 は有効な Win32 アプリケーションではありません。→64bit Pythonに対して 32bitのMeCabがインストールされていないか
  • ImportError: DLL load failed: 指定されたモジュールが見つかりません。 →MeCabのbin\ に環境変数PATHが通っているか

Pythonで形態素解析をしてみる

ようやく環境が整いました。それでは形態素解析をしてみましょう。

import MeCab

tagger = MeCab.Tagger("-Ochasen")
result = tagger.parse("私は東京に住んでいます。")
print(result)

素晴らしい。きちんと文章が品詞ごとに分別されていますね。
f:id:sanvarie:20190706124711p:plain

出力フォーマットの変更

上記のコードで-Ochasenを使用していますが、これは出力フォーマットのオプションの一つです。公式サイトによるとオプションは以下の三つのようです。

  1. -Oyomi (ヨミ付与)
  2. -Ochasen (ChaSen互換)
  3. -Odump (全情報を出力)
  • Oyomiを使用するとこのような結果になります。

f:id:sanvarie:20190706125109p:plain

  • Odumpを使用するとこのような結果になります。

f:id:sanvarie:20190706125159p:plain

未知語推定

MeCab は, 辞書に単語が未登録の場合でも適当にその品詞を推定してくれるみたいです。試しに「ひょっこりはん」で形態素解析してみようと思います。果たしてどのような結果になるのでしょうか。

import MeCab

tagger = MeCab.Tagger("-Ochasen")
result = tagger.parse("ひょっこりはん")
print(result)

このような結果になりました。たしかにこれを人名と分類するのは難しいですが、何かしらの結果は返してくれるみたいですね。
f:id:sanvarie:20190706125641p:plain

Pythonでもっと形態素解析をしてみる

これだけではつまらないので、もう少し突っ込んだところまでやってみようと思います。

今回使用する文章

今度の打ち合わせは永田町にある先方のオフィスで行うことになりました。


なぜ永田町なのかは秘密です。文章に地名が含まれている場合、その地名を抽出するプログラムを書いてみます。ちなみに単純に形態素解析の結果を出力するとこのようになります。

f:id:sanvarie:20190706141128p:plain

以下が作成したプログラムです。少しだけ複雑になりましたが、ポイントさえつかめばすぐに使いこなすことができるかと思います。

import MeCab

tagger = MeCab.Tagger("-Ochasen")
node = tagger.parseToNode("今度の打ち合わせは永田町にある先方のオフィスで行うことになりました。")

while node:
    if node.surface != "":  # ヘッダとフッタを除外
        if (node.feature.split(",")[1] == "固有名詞") and (node.feature.split(",")[2] == "地域"):
            print(node.feature.split(",")[6])
    node = node.next

少しだけポイントを解説します。

parseToNode メソッド

このメソッドを使用することにより、「文頭」という特別な形態素が MeCab::Node クラスのインスタンスとして 取得できます。要するにこのメソッドを使えば形態素解析の結果で色々できるということですね。

その他のポイント

  1. parseToNode メ ソッドで取得した node をwhile でループ。node.next で次の行に行きます。
  2. node.surface は形態素解析の結果のヘッダとフッタが取得できます。上記画像の一番上と下の「BOS BOS/EOS」の行ですね。
  3. node.feature で品詞の取得ができます。

結果

「永田町」のみを抽出することができました。
f:id:sanvarie:20190706141316p:plain

本当に素晴らしい技術ですね。おそらくAI分野でこういった技術が数多く使われているんじゃないかと思います。自分ももう少しこういったことを勉強して GISとAI を使って何か面白いことをしてみたいですね。

本エントリーの後に以下を書いてみました。Mecab の具体的な使用例になるのではないかと思います。興味のある方はぜひ読んでみてください。
www.gis-py.com

本日は以上です。

ArcGIS Pro でサッカーのヒートマップを作ってみよう

本日はArcGIS Pro でサッカーのヒートマップを作ってみようと思います。久しぶりのサッカーネタですね。以前は高校サッカーでしたが今度はプロの試合を分析してみようと思います。

www.gis-py.com

ただ、今回はプログラムなしで、ArcGIS Pro オンリーで行きたいと思います。

データに関して

かなり前からサッカー関連のデータを色々探していたのですが、なかなか見つからず途方に暮れていたところ、ついにいい感じのデータを見つけてしまいました!Soccer video and player position datasetというサイトがそのデータを公開していて、学術関係の論文発表の際に作成されたデータのようです。論文に関してはこちらを参照してください。http://home.ifi.uio.no/paalh/publications/files/mmsys2014-dataset.pdf

QGISで同じようなことをしていらっしゃる方がいましたので、そちらを参考にさせていただきました。

qiita.com

出典

"Soccer video and player position dataset": S. A. Pettersen, D. Johansen, H. Johansen, V. Berg-Johansen, V. R. Gaddam, A. Mortensen, R. Langseth, C. Griwodz, H. K. Stensland, and P. Halvorsen, in Proceedings of the International Conference on Multimedia Systems (MMSys), Singapore, March 2014, pp. 18-23

ちなみに、このデータを使用する場合は出典を明記する必要がありますので、お気を付けください。

データの概要

  1. サッカーの試合でプレイヤートラッキングシステムを使ってボールとプレイヤーの位置関係を抽出したもの
  2. 場所はアールヴヘイム・スタディオン(ノルウェー一部リーグに所属するトロムソILのホームスタジアム)
  3. データは3試合分公開(前半と後半にデータが分かれている)
  4. CSVで提供されている

使用するデータ

データはこのように公開されています。上から三つ目までがCSVになっています。20Hzが0.05秒間隔で1Hzが1秒間隔で選手とボールの位置を抽出しています。今回はZXY data raw (20Hz)を使用します。
f:id:sanvarie:20190705223233p:plain

※ZXY data raw (20Hz)はセンサーから受け取ったそのままのデータ、ZXY data interpolated (20Hz)はその後処理加工されたデータみたいですね。

試合について

使用するデータですが、トロムソIL vs FCアンジ・マハチカラのものです。

※FCアンジ・マハチカラはロシア・プレミアリーグ所属。かつてロベルト・カルロスも所属していたチームですね。

データを開いてみた

データを開いてみると前半部分に355553件、後半部分に360677件のデータがあることがわかりました。
f:id:sanvarie:20190705225119p:plain

ヘッダーがないので、以下のようにヘッダーをつけました。

f:id:sanvarie:20190705225929p:plain

ArcGIS Pro でCSVの取り込み

XYテーブル→ポイントのジオプロセシングツールを使います。
f:id:sanvarie:20190705231110p:plain

そうすると、このように取り込みができました。それっぽい感じですね。
f:id:sanvarie:20190705231414p:plain

特定の選手を抽出

tag_idというカラムが選手を表しているらしいのですが、どのtag_idがどの選手と紐づいているかは匿名性のため公表してないそうです。なので、とりあえず今回はtag_idが10の選手にフォーカスをあててみたいと思います。

まずはレイヤプロパティ画面でtag_id=10の選手のみが抽出されるように設定します。
f:id:sanvarie:20190705234446p:plain

それっぽい感じになりましたね。
f:id:sanvarie:20190705234529p:plain

ヒートマップの作成

次にシンボル設定でヒートマップを選択します。
f:id:sanvarie:20190705234626p:plain

このようになりました。
f:id:sanvarie:20190705234824p:plain

背景の作成

参考サイトと同様にFile:Football pitch metric.svg - Wikimedia Commonsの画像を使用します。

f:id:sanvarie:20190705235918p:plain

ArcGIS Proに画像を取り込みましたが、位置があっていませんね。
f:id:sanvarie:20190706000231p:plain

そこでジオリファレンスですね。ジオリファレンスについてはArcGIS Pro の画像処理機能紹介 その 1 ~ジオリファレンス~ | ArcGISブログを参照してください。

f:id:sanvarie:20190706000317p:plain


かなりそれっぽくなりましたね!
f:id:sanvarie:20190706000821p:plain

3D上でのヒートマップ

ArcGIS Proの特徴の一つでもある3Dの機能を使ってヒートマップを確認してみようと思います。

ローカル シーンに変換します。
f:id:sanvarie:20190706001128p:plain

そして、「3D シンボルを実世界の単位で表示」をオフにします。
f:id:sanvarie:20190706002235p:plain

ばっちりですね!
f:id:sanvarie:20190706003715p:plain



後半も表示してみました。
f:id:sanvarie:20190706003631p:plain

感想

思ったよりも簡単にヒートマップを作成することができました。さすがはArcGIS Proといったところでしょうか。できればこういったデータがもっとたくさん転がっているいいのですが。好きな選手の動きとかを自分で分析したりするのはとても楽しそうですしね。また何かいいデータがあったらサッカーネタで何か書きたいと思います。本日は以上です。

pyodbcを使ってPythonでDBを操作しよう ~SQL Server編~

本日はPythonでODBC経由でDBを操作する方法を紹介したいと思います。RDBMSといってもOracleやMySQLなど色々あるかと思いますが、今回はSQL Serverを使ってみようと思います。

使用するライブラリ

pyodbcというライブラリをインストールする必要があります。

pip install pyodbc でインストールできると思いますが、インストールがうまくいかない場合は以下サイトでwheelをダウンロードしてそれを使ってインストールしてください。
pypi.org

実行環境

Windows10
SQL Server2016
Python3.6.5

GISとRDBMS

余談にはなるのですが、GISを使われている方の中でどのくらいの方がGISデータの管理にRDBMSを使われているのかすごく気になります。

  1. データはあくまでShapeやGDBのような形式でのみ管理している
  2. ShapeやGDBなどのGISデータとそれに紐づく属性情報などをRDBMS上で管理している
  3. ShapeやGDBなどのGISデータとそれに紐づく属性情報などをRDBMSではなく、MongoDBなどのNoSQLで管理している
  4. ArcGIS Enterprise Geodatabaseで図形も属性もRDBMS上管理している

大雑把ではあるのですが、このようにカテゴライズできるかと思います。皆様はどのような形でGISデータを管理されているのでしょうか?ちなみに私は「2」です。

テストデータ

以下のようなデータを使用します。

  • DB・・・Sample
  • テーブル・・・TEST

f:id:sanvarie:20190629160312p:plain

レコードは一つだけ作りました。
f:id:sanvarie:20190629161353p:plain

pyodbcの使用方法

pyodbcを使用して以下のような操作をするための方法を書いてみました。

  1. ログイン
  2. 検索
  3. 更新
  4. 挿入
  5. 削除

ログイン

接続文字列を作成してpyodbc.connectを実行すればログインができます。

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

def login():
    # インスタンス
    instance = "インスタンス"

    # ユーザー
    user = "ユーザー"

    #パスワード
    pasword = "パスワード"

    #DB
    db = "Sample"

    connection = "DRIVER={SQL Server};SERVER=" + instance + ";uid=" + user + \
                 ";pwd=" + pasword + ";DATABASE=" + db


    return pyodbc.connect(connection)

if __name__ == '__main__':
    con = login()

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

検索

検索はログイン処理で取得したコネクションとSQLをselect_executeメソッドで処理しています。そして、取得したレコードをループし(1行しかありませんが)、各カラムの中身をprintしています。

def select_execute(con, slq):
        cursor = con.cursor()
        cursor.execute(slq)
        rows = cursor.fetchall()
        cursor.close()

        return rows

if __name__ == '__main__':
    con = login()
    sql =  '''select *
                from TEST'''
    res = select_execute(con, sql)
    for r in res:
        print(r[0])
        print(r[1])
        print(r[2])
        print(r[3])

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

f:id:sanvarie:20190629162754p:plain

更新

NAMEカラムの中身を「ライン」から「ラインテスト」となるようにデータを更新します。要領は検索と同じですが、更新処理なので最後にcommitする必要があります。

def updatet_execute(con, slq):
    cursor = con.cursor()
    cursor.execute(slq)
    con.commit()

if __name__ == '__main__':
    con = login()
    sql = """UPDATE TEST
                SET NAME = NAME + 'テスト'"""

    # データ更新
    updatet_execute(con, sql)

想定通りの結果になりました。
f:id:sanvarie:20190629165149p:plain

挿入

1レコード挿入する処理を書いています。

def insert_execute(con, slq):
    con.execute(sql,
                2,
                "ポイント",
                "マルチ",
                "特になし"
                )

    con.commit()

if __name__ == '__main__':
    con = login()
    sql = """INSERT
             INTO TEST(FEATURE_NUM,
                       NAME,
                       TYPE,
                       NOTE
                       )
                  VALUES (?,
                          ?,
                          ?,
                          ?
                          )"""

    # データ登録
    insert_execute(con, sql)

ちゃんと挿入した分のレコードがDB側に反映されていました。
f:id:sanvarie:20190629164159p:plain

削除

上記処理で挿入したレコードを削除しています。

def delete_execute(con, slq):
    cursor = con.cursor()
    cursor.execute(slq)
    con.commit()

if __name__ == '__main__':
    con = login()
    sql = """DELETE
               FROM TEST
              WHERE FEATURE_NUM = 2"""

    # データ削除
    delete_execute(con, sql)

想定通りレコードが削除されました。
f:id:sanvarie:20190629165718p:plain

pyodbcを使えばとても簡単にPythonからSQL Serverを操作できることがわかりました。今回はあくまでサンプルなので、クラスは作りませんでしたが、実際にDB操作をする際はDBアクセス用のクラスなどを作ってそこに各処理を記載するのがいいかと思います。今後はGISデータとDBのデータを一緒に扱って何かしらの処理をするようなことも紹介できればと思います。

本日は以上です。

Pythonで作ったモジュールのパッケージ化をしてみよう

GISはあまり関係なくPython一般のお話になってしまうのですが、本日はPythonでのパッケージ化について書いてみようと思います。

このブログを読んでくださっている方の中でPython経験者は多数いらっしゃるかと思います。ただ、以下のような理由でパッケージ化をしたことがないという方も一定数いらっしゃるのではと思っています。

  1. パッケージ化するほどの大きなプログラムを書く必要がない
  2. パッケージ化しなくてもやりたいことができるので覚える必要がない
  3. そもそもパッケージ化できることを知らなかった
  4. etc

たしかに、開発するプログラムが小規模の場合、パッケージ化の効果は少ないかもしれませんが、大規模化して同じメソッドやクラスを再利用する回数が増えると、パッケージ化の効果は大きくなります。もし、将来的に大規模開発を行う、もしくは、行いたいという場合は覚えていておいて損はないかと思います。

パッケージの作成

以下のようなディレクトリを作ってみました。
f:id:sanvarie:20190629135525p:plain

main.pyが今回メインとなる処理のファイルですね。

libフォルダの中身はこのようになっています。
f:id:sanvarie:20190629134926p:plain

featureclass.pyがmain.pyでインポートされるモジュールです。

パッケージの初期化ファイル (__init.py__)

もう一つ、__init__.pyというのがあるのですが、これは空のファイルです。Python 3.3 より前のバージョンでは、パッケージディレクトリに __init__.py というファイルを置かなければ、そのディレクトリをパッケージとして認識させることができませんでした(ImportError: No module named ... というエラーが発生する)。Python 3.3 より前のバージョンをお使いの方は__init__.pyを必ず作成してください。

全体の構成はこのようになります。
f:id:sanvarie:20190629140958p:plain

モジュールのインポート

featureclass.pyのインポートをする前にfeatureclass.pyでどのような処理をしているのかを解説します。一応GISがメインのブログなのでArcPyを使用したプログラムにしてみました。

featureclass.py
# -*- coding: utf-8 -*-
import arcpy

class FeatureClass():

    def __init__(self):
        pass

    def get_shape_type(self, fc):
        """対象フィーチャクラスのshapeTypeを取得します。
           アノテーションのFeatureTypeはPolygon→Annotationに変換します

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

        returns:
            geo_type -- ジオメトリタイプ
        """

        fc_type = arcpy.Describe(fc).FeatureType

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

            geo_type = "Annotation"

        else:

            geo_type = arcpy.Describe(fc).shapeType

        return geo_type

get_shape_typeというメソッドのみ実装したシンプルなクラスですね。このメソッドは引数として受け取ったフィーチャクラス名からそのフィーチャクラスのジオメトリタイプをreturnします。これをmain.pyで呼び出すとこのようになります。

main.py
# -*- coding: utf-8 -*-
import arcpy
from lib.featureclass import FeatureClass

class Main():

    def __init__(self):
        pass

    def main(self):
        arcpy.env.workspace = r"D:\python\data\Sample.gdb"
        featureclass = FeatureClass()

        # ワークスペースのフィーチャクラスのジオメトリタイプを取得
        for f in arcpy.ListFeatureClasses():
            print(featureclass.get_shape_type(f))

if __name__ == '__main__':
    m = Main()
    m.main()

ポイントとしては以下3点です。

  1. from lib.featureclass import FeatureClassでモジュールをインポート
  2. featureclass = FeatureClass()でインスタンス化
  3. 作成したインスタンスでget_shape_typeを実行

Sample.gdbにあるフィーチャクラス分ループを行い、get_shape_typeで取得したジオメトリタイプをprintするという処理です。使用したデータは以下のフィーチャクラスです。

f:id:sanvarie:20190629141206p:plain

実行結果はこのようになりました。
f:id:sanvarie:20190629141636p:plain

これだけです。とても簡単ですね。これをすることで共通的な処理などをパッケージ化して作業の効率化などを図ることができるかと思います。また、少しハードルが上がるのですが、PYPIに作成したパッケージを登録することもできるかと思います。夢が広がりますね。

pypi.org

本日は以上です。