waigani's diary

QGISを中心にFOSS4Gをいじくる

QGIS プロセッシング Scriptsを使いたい その4

QGIS processing

QGisLayers.py

いくつか用意されているprocessingの便利クラスのうち、QGisLayersを使ってみましょう。
動作を確認するには、pythonコンソールから試すのがお手軽です。
まずはPythonコンソールを立ちあげてもらって、準備として必要なモジュールを読み込んでおきます。QGisLayers()インスタンスの作成もしておきましょう。

from processing.core.QGisLayers import QGisLayers
qlayer = QGisLayers()

それではいくつかのメソッドを試してみましょう。

getSupportedOutputVectorLayerExtensions()

アウトプットとして対応しているベクトルファイル形式の拡張子を返してくれます。

qlayer.getSupportedOutputVectorLayerExtensions()

拡張子のリストが返ってきます。

['shp', u'dxf', u'gxt', u'gml', u'gmt', u'geojson', u'itf', u'itf', u'xml', u'csv', u'dgn', u'bna', u'tab', u'gpx', u'gdb', u'sqlite', u'kml', u'000']

qlayer.getSupportedOutputRasterLayerExtensions()

アウトプットとして対応しているラスタファイル形式の拡張子を返してくれます。

qlayer.getSupportedOutputRasterLayerExtensions()

拡張子のリストが返ってきます。

['tif', 'grd', 'sdat', 'bt', 'img', 'ter', 'pix', 'hdr', 'vrt', 'rgb', 'bmp', 'rsw', 'nc', 'gsb', 'gen', 'gtx', 'mpr', 'mpl', '', 'pnm', 'ntf', 'rst']

getSupportedOutputTableExtensions()

アウトプットとして対応しているテーブルの拡張子を返してくれます。

qlayer.getSupportedOutputTableExtensions()

これは

['csv']

と'csv'と入ったリストが返ってくるだけ。ソースを見てもこれしか書いてないです。

getRasterLayers()

読み込む済みのラスタレイヤの一覧を返してくれます。

qlayer.getRasterLayers()

ラスタレイヤの一覧が返ってきます。

[<qgis.core.QgsRasterLayer object at 0x0000000009DE7B70>, <qgis.core.QgsRasterLayer object at 0x0000000009DE7BF8>, <qgis.core.QgsRasterLayer object at 0x0000000009DE7C80>]

getVectorLayers()

読み込む済みのベクトルレイヤの一覧を返してくれます。

qlayer.getVectorLayers()

ベクトルレイヤの一覧が返ってきます。

[<qgis.core.QgsVectorLayer object at 0x0000000009DE7D08>, <qgis.core.QgsVectorLayer object at 0x0000000009DE7D90>, <qgis.core.QgsVectorLayer object at 0x0000000009DE7E18>]

引数としてgeometryTypeをリストにして指定することが出来ます。一致するレイヤのみ返ってきます。
ポリゴンのレイヤのみ貰いたい場合は、下記の指定になります。

qlayer.getVectorLayers([QGis.Polygon])

その他

その他にも便利なメソッドが用意されていますので、QGISインストールディレクトリ/python/plugins/processing/core/QGisLayers.pyの中を覗いて確認してみましょう。

QGIS プロセッシング Scriptsを使いたい その3

QGIS processing

指定したフィールドのみ残したベクトルファイルの作成

簡単なスクリプトをもう１つ。
行っている処理としては、

  • 入力としてベクトルファイルを受け取る(inputに入力ファイル名が入る)
  • 入力ベクトルレイヤから、残したいフィールド名を1つだけ受け取る(fieldNameに入る)
  • processing.getobject()でファイル名→QgsVectorLayer
  • 残したいフィールドの定義だけ抜き出し、出力ファイルに設定するたまにQgsFieldsを用意
  • 残したいフィールド以外は、名前をリストに入れておく、後で削除用に使用
  • processingのVectorWriterクラスを利用してベクトルファイル作成
  • 入力ベクトルレイヤの図形要素に順次アクセス
  • 要らないフィールドは削除した状態で、出力ベクトルファイルへ追加

というものです。

from qgis.core import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from processing.core.VectorWriter import VectorWriter

##[Test scripts]=group
##input=vector
##fieldName=field input
##output=output vector

vectorLayer = processing.getobject(input)
provider = vectorLayer.dataProvider()
fields = provider.fields()
newFields = QgsFields()
newFields.append(fields.field(fieldName))

deleteFields = []
for field in fields:
  if field.name() != fieldName:
    deleteFields.append(field.name())

writer = VectorWriter(output, provider.encoding(), newFields, provider.geometryType(), vectorLayer.crs())

features = vectorLayer.getFeatures()
for feature in features:
  for deleteField in deleteFields:
    feature.deleteAttribute(provider.fieldNameIndex(deleteField))
  writer.addFeature(feature)

del writer

追記 2014.7.24
QGIS2.4では、processing.getobject()がprocessing.getObject()になっているので注意！！

id,nameの2つのフィールドを持ったベクトルデータがあったとして、

スクリプトを実行して、nameのみを残すことを選択します。

実行結果として、フィールドnameのみ残されたベクトルレイヤが作成されます。

processingで用意されているクラス

ラスタファイルを作成するのに使用したRasterWriter、ベクトルファイルを作成するのに使用したVectorWriterはprocessingで用意されたクラスになります。

　QGISインストールディレクトリ/python/plugins/processing/core

以下に様々なファイルがありますので覗いてみましょう。
RasterWriterが入っているRasterWriter.py、VectorWriterが入っているVectorWriter.pyともに短いソースですので内容を確認してみることをお勧めします。引数として何を渡していたかも、すぐにわかると思います。

QGIS プロセッシング Scriptsを使いたい その2

QGIS processing

GROUP

新しいスクリプトを作成すると、"User scripts"以下に追加されるのは前回見てもらったとおり。
これを変更して、既存のグループもしくは新規のグループに入れることも出来ます。
指定方法は、

##[Test scripts]=group

のようにスクリプト内に記載しておくだけです。下図のように"Test scripts"の下に作成したスクリプトtestが表示されます。

出力の指定

スクリプトの実行結果として何を出力するかは、入力に使用したパラメータと同様な記述で定義出来ます。
定義出来る種別としては、

  • output raster
  • output vector
  • output table
  • output html
  • output file
  • output number
  • output string

があります。
書式としては入力パラメータと同じく、

##パラメータ名=パラメータタイプ

となります。
試しに、

##output = output vector

とスクリプトに記述して実行してみましょう。

出力ファイル名を入力するためのフォームが表示されます。
ファイル名を指定しなかった場合、テンポラリのファイル名が与えられますので、入力は必須ではありません。
raster,vector,table,html,fileについては、ファイル名選択のフォームが表示されますが、number,stringについてはフォームは表示されません。

ベクトルレイヤ範囲のラスタ作成

簡単なスクリプトを作ってみます。
処理としては、

  • 入力としてベクトルファイルを受け取る(inputに入力ファイル名が入る)
  • processing.getobject()でファイル名→QgsVectorLayer
  • vectorLayerの範囲を取得しておく
  • processingのRasterWriteクラスを使用して、ベクトルレイヤの範囲を包括するラスタ作成(値は全面に0しか入れていない)

というものです。

from processing.core.RasterWriter import RasterWriter

##[Test scripts]=group
##input=vector
##output=output raster

vectorLayer = processing.getobject(input)
extent = vectorLayer.extent()
xmin = extent.xMinimum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()
ymax = extent.yMaximum()
cellSize = (xmax-xmin)/100.0

writer = RasterWriter(output, xmin, ymin, xmax, ymax, cellSize, 1, vectorLayer.crs())
for x in range(writer.nx):
  for y in range(writer.ny):
    writer.setValue(0.0, x, y)


writer.close()
del writer

追記 2014.7.24
QGIS2.4では、processing.getobject()がprocessing.getObject()になっているので注意！！

QGIS プロセッシング Scriptsを使いたい その1

QGIS processing

プロセッシング Scripts

QGIS 2.0より、SEXTANTEからプロセッシングと名前が変更になった機能を使ってみましょう。もちろんコードをガリガリ書いて使います。
メニューから"プロセッシング"→"ツールボックス"を選択します。
出てきたツールボックスのツリーの一番を開いていくと、"Create new script"というのが現れます。

"Create new Script"をダブルクリックすると、エディタが現れます。ここにプロセッシングの作法を守ってpythonコードを書いていけば、オリジナルのアルゴリズムを足していくことが出来ます。

もしくは、実際に書いたスクリプトが保存されるのは、

 ユーザのホームディレクトリ\.qgis2\processing\scripts

以下になりますので、他のエディタを使ってスクリプトを記述して、こちらにおいてしまってもOKです。

パラメータを受け取る

先頭"##"で記述した行で、入力フォームを作成してパラメータを受け取ることが出来ます。
書式としては、

##パラメータ名=パラメータタイプ 値

となります。

試しに

##inpupt = vector

とだけエディタに記述して保存してみましょう。test.pyとして保存したとします。
プロセッシングツールボックスを見ると、ツリーの一番下"User scripts"に保存したスクリプトが追加されていることが確認出来ます。

追加したスクリプトをダブルクリックして実行してみましょう。
パラメータタイプとして"vector(ベクトルレイヤ)"を指定していますので、読み込み済みのベクトルレイヤがあればレイヤ名が表示されますし、ファイルを指定することも出来ます。

レイヤを選択して"run"をしても、いまは処理の中身を記述していないので何も起きません。内部の処理としては、inputに選択したレイヤのファイル名（レイヤを選択した場合)、もしくは選択したファイル名が入っただけで終わりとなります。

パラメータタイプとして指定できるものとしては、下記があります。

  • vector ベクトルファイル
  • raster ラスタファイル
  • table テーブル
  • number 数値、val = number 10のようにデフォルト値を指定できる
  • string 文字列、str = string hogehogeのようにデフォルト値を指定できる
  • boolean 論理型、val = boolean Trueのようにデフォルト値を指定できる
  • multiple raster 複数のラスタレイヤ
  • multiple vector 複数のベクトルレイヤ
  • folder フォルダ名
  • file ファイル名
  • field フィールド名、ベクトルファイル名をinputとして受け取っていた場合、field_name = filed inputのように指定します
##input=vector
##field_name=field input


DEBUG

さて入力のパラメータは受け取れることがわかりました。スクリプトを書いていく前に、DEBUG方法を確認しておきたいです。調べたのですが、いい方法が良くわからず。ログに逐一書き出すしかないのかな？と思っています。
プロセッシングのログは、

 ユーザのホームディレクトリ\.qgis2\processing\processing_qgis.txt

に書き出されます。ログに追記したい場合は、

from processing.core.ProcessingLog import ProcessingLog

ProcessingLog().addToLog(ProcessingLog.LOG_INFO, 'test output')

のようにしてあげます。

あれ？スクリプトのエディタって日本語の入力受け付けないのかな？？とか、また余計なところが気になってきた。。

QGISでTypographicMap

QGIS

続きです

これで本当に最後の予定です。
せっかく2回やったのでまとめを行います。
QGISのSVG塗りつぶしパターンを属性から作ってみる
QGISのSVGマーカーを属性から作ってラインに置く

使用するデータ

基盤地図情報ダウンロードサービスから縮尺レベル2500の行政区画・軌道の中心線・道路縁を使用しました。
付属のツールでshapefileへの変換を行い、カラム名称を日本語からアルファベットに変換はしています。
行政区画についてはnameカラムに行政名が入るように、軌道の中心線・道路縁についてはsortカラムに種別が入るようにしています。

結果

スケールでフィルタしているので、ひいているときは行政区画だけ表示。


行政区画だけ表示、もう少し拡大。


拡大していくと道路も表示。完全にラインに沿っての表示は出来ていない。ライン上にマーカーを置いているので、はみ出している文字が出てきます。改善したいですが、この方法だと無理ですね。


五稜郭。

書いたソース

QGIS2.0の、pythonコンソールから実行します。"エディタの表示"をしておいて、ソースを貼り付けて実行しています。
前提条件として、埋め込んでしまって、

  • ラインは"sort"カラム、ポリゴンは"name"カラムの属性を描画

としています。
処理としては、

  • レイヤーを順次アクセス
  • ベクトルレイヤのうち、ラインとポリゴンのみが対象
  • 属性値のリストを作っておく
  • 属性値毎のSVGファイルを作成(temporaryファイルですが消さず)
  • 図形種別にあわせたシンボル作成、色はランダム、ポリゴン塗りつぶしについては角度もランダム
  • "ルールに基づいた"スタイルにフィルターと合わせてラインシンボルをセット
  • 全ての属性値に対するルールの設定が終わったら、スタイルをrendererに入れてlayerに設定
  • 再描画

といった流れです。
あいかわらず変なやり方なので、スタイルの設定の参考程度に。

# -*- coding: utf-8 -*-
from PyQt4 import QtGui, QtCore
import tempfile
import random

#描画サイズ
lineSymbolSize = 3
polySymbolSize = 4

def makePolySymbol(attr, svg, size):
    #空のシンボルを作成
    symbol = QgsFillSymbolV2()
    symbol.deleteSymbolLayer(0)
    
    #属性を入れたSVGを指定して、SVG塗りつぶしを作成
    fillLayer = QgsSVGFillSymbolLayer(svg.name, len(attr)*size)
    subSymbol = fillLayer.subSymbol()
    subSymbol.deleteSymbolLayer(0)
    
    #シンボルに塗りつぶしを追加
    symbol.appendSymbolLayer(fillLayer)
    #適当に回転させる
    symbol.setAngle(random.randint(0,90))
    return symbol

def makeLineSymbol(attr, svg, size):
    #SVG読み込み
    markerSymbolLayer = QgsSvgMarkerSymbolLayerV2 (svg.name, len(attr)*size)
    
    #マーカーライン作成、SVGシンボルセット
    markerLine = QgsMarkerLineSymbolLayerV2()
    svgSymbol = markerLine.subSymbol()
    svgSymbol.deleteSymbolLayer(0)
    svgSymbol.appendSymbolLayer(markerSymbolLayer)
    markerLine.setInterval(len(attr)*size)

    symbol = QgsLineSymbolV2()
    symbol.deleteSymbolLayer(0)
    symbol.appendSymbolLayer(markerLine)
    return symbol

#レイヤーを順次アクセス
layers = iface.legendInterface().layers()
for layer in layers:
    if layer.type() != layer.VectorLayer:
        continue

    if layer.geometryType() != QGis.Line and layer.geometryType() != QGis.Polygon:
        continue
   
    #描画する属性フィールド
    key = "sort" if layer.geometryType() == QGis.Line else "name"

    #空の"ルールに基づいた"スタイル作成
    rule = QgsRuleBasedRendererV2.Rule(None)

    #指定の属性フィールドの値リストを作成
    attrList = []
    features = layer.getFeatures()
    for feature in features:
        if feature[key] not in attrList:
            attrList.append(feature[key])

    for attr in attrList:
        #指定の属性を入れたSVGファイルを作成
        svgPattern = """<?xml version="1.0" standalone="no"?>
            <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" 
             "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">"""
        svgPattern += '<svg width="%dpx" height="1px" ' % len(attr)
        svgPattern += 'xmlns="http://www.w3.org/2000/svg" version="1.1">'
        svgPattern += '<rect x="0" y="0" width="%d" height="1" fill="#FFFFFF" />' % len(attr)
        svgPattern += '<text x="0" y="1" font-family="MS-Mincho" font-size="1" fill="rgb(%i,%i,%i)" >' % (random.randint(0,255), random.randint(0,255), random.randint(0,255))
        svgPattern += '%s' % attr.encode('utf-8')
        svgPattern += '</text>'
        svgPattern += '</svg>'
    
        tempSvg = tempfile.NamedTemporaryFile(mode="w+t", delete=False)
        tempSvg.write(svgPattern)
        tempSvg.close()

        symbol = makeLineSymbol(attr, tempSvg, lineSymbolSize) if layer.geometryType() == QGis.Line else makePolySymbol(attr, tempSvg, polySymbolSize)

        #指定の属性カラムによるフィルタを設定して、ルールを作成、追加
        #ラインとポリゴンで表示スケール変更
        scaleMax = 2500 if layer.geometryType() == QGis.Line else 200000
        ruleAdd = QgsRuleBasedRendererV2.Rule(symbol, label=attr, filterExp="\"%s\" = '%s'" % (key, attr), description=attr, scaleMaxDenom=scaleMax)
        rule.appendChild(ruleAdd)
    
    renderer = QgsRuleBasedRendererV2(rule)

    layer.setRendererV2(renderer)

iface.mapCanvas().refresh()

iface.legendInterface().refreshLayerSymbology(layer)

QGISで属性の値に応じてポリゴン分割？ ちょっとつづき

前回の

QGISで属性の値に応じてポリゴン分割？の続きです。折角なので少しだけ絵を作って遊んでみます。

作ったもの

男女比

住民基本台帳人口・世帯数から住民基本台帳人口・世帯数 平成25年9月を取得。
北海道内人口男女比をプロット。良いシンボルが無かったのでWCで失礼します。
男性:2,577,703人 女性:2,868,321人で女性の比率の方が少しだけ高い。図にするとわかりずらいけど。

振興局別人口

男女比と同じデータを使用。振興局別の人口比率をプロット。
石狩が約半分。

天気日数

日本統計年鑑から第1章　国土・気象　1 -10 気象官署別日照時間，天気日数を取得。札幌の晴れ、曇、雪、雨の日数比率をプロット。
晴れ少ない。。

QGISで属性の値に応じてポリゴン分割？

QGIS

うまくタイトル考えられません。
Typography Mapはおいておいて、相変わらずInformation Graphicsを眺めております。

Information Graphics

Information Graphics

  • 作者: Sandra Rendgen,Paolo Ciuccarelli,Richard Saul Wurman,Simon Rogers,Julius Wiedemann
  • 出版社/メーカー: Taschen America Llc
  • 発売日: 2012/05/27
  • メディア: ハードカバー
  • クリック: 2回
  • この商品を含むブログ (1件) を見る
そしてやってみたくなったのが、下記のような行政区画を統計値を使って割って、割合を示すような地図です。

使用するデータ

ちょうど手元にあった北海道の図形を使います。
これに北海道の住民基本台帳人口・世帯数のページから年齢5歳階級別人口を持ってきて、属性として付けました。

北海道に、

こんな感じで5歳区切りの人口の属性を持たせています。

結果

ポリゴン内に一定間隔のポイントを発生させます。
このポイントを、5歳区切りの人口の割合に応じて、属性分けしてあげました。
0-4歳の人口が、総人口の3%なら、発生したポイントのうち3%に"0-5歳"と属性を入れてあげます。
そして属性で色分け。

もう少し階層を少なく、そして色分けをくっきりしないと、なんだかわからないですね。何かいいデータを見つけたら、あらためて絵を作ってみよう。

書いたソース

QGIS2.0の、pythonコンソールから実行します。"エディタの表示"をしておいて、ソースを貼り付けて実行しています。
前提条件として、

  • レイヤーは選択済み
  • 図形は1つしかないこと
  • 入っている属性がすべてintegerであること
  • すべての属性を使います、使わない属性が入っていないこと

としています。
処理としては、

  • 選択済みのレイヤーの属性定義を取っておく、作成したレイヤの属性にカラム名を入れるため
  • 図形を1つ持ってくる、各属性値を属性値の総和で割って、割合を求めておく
  • ポイントを一定間隔で発生
  • 発生したポイントレイヤにカラムを1つ追加
  • 位置でソートされていることを期待してポイントに順次アクセス、同じ属性を、求めておいた割合＊総ポイント数分の図形に入れていく

といった流れです。
色付けは手動でやっています。

# -*- coding: utf-8 -*-
from PyQt4 import QtGui, QtCore
import processing

#作成六角形サイズ
size = 0.1

#属性を入れるカラム名
sort = u'sort'

#レイヤーが選択されていることが前提
layer = iface.activeLayer()
if layer.featureCount() != 1:
    print "Warnnig:this program allow to run on condition layer has one feature"

#1つだけ図形を持ってくる
layer.selectAll()
inputFeature = layer.selectedFeatures ()

#属性定義をとっておく
fields = inputFeature[0].fields()

#各属性の割合を調べておく
attrList = inputFeature[0].attributes()
all = sum(attrList)
attrListFloat = map(lambda x:float(x)/all, attrList)

#processingを使って、六角形作成、中心点作成、必要範囲のみ切り取り
#processingの必要はなし、使ってみたかっただけ
hexResult = processing.runalg('script:hexgridfromlayerbounds', layer, size, None)
centroidsResult = processing.runalg('qgis:polygoncentroids', hexResult['grid'], None)
intersectionResult = processing.runalg('qgis:intersection', centroidsResult['OUTPUT_LAYER'], layer, None)
pointLayer = processing.load(intersectionResult['OUTPUT'])

pointLayer.startEditing ()

#属性を入れるカラムを追加
pointLayer.addAttribute(QgsField(sort, QtCore.QVariant.String, u'string', 10))

#図形はソートされていることを期待して、id順にアクセス
featureCount = pointLayer.featureCount()
fieldIndex = 0
fieldCnt = 0
for i in range(featureCount):
    #調べておいた割合にしたがって、図形に属性を入れていく
    #大体の数ということで
    pointLayer.changeAttributeValue(i, pointLayer.fieldNameIndex(sort), fields.field(fieldIndex).name())
    fieldCnt += 1
    if fieldCnt >= featureCount * attrListFloat[fieldIndex]:
        fieldIndex += 1
        if fieldIndex >= len(attrListFloat):
            break
        fieldCnt = 0
    
pointLayer.commitChanges ()

ToDo

面白い絵が出来そうなデータを見つけて、絵を作ってみる。

以上