PyQGISでレイヤを裁断するプログラムを作った話（FOSS4GのLTのリベンジ)

この記事は、「FOSS4G Advent Calendar 2018」の9日目の記事です。
＃FOSS4G Tokyoのライトニングトークでお話させていただきましたが、改めて

今年はシュレッダーが流行ったらしい（と思ってる）

元ネタはイギリスの芸術家バンクシーが描いた絵が、オークション会場で落札された瞬間、会場客の目の前でシュレッダーで切り刻まれたというニュース
https://gigazine.net/news/20181009-banksy-artwork-shredded/

結局はこのシュレッダーは作者であるバンクシー本人が仕組んだものなわけだが、これが各方面でパロディされている。
https://www.mag2.com/p/news/374280

・・・なるほど。いまシュレッダーが流行っているのか。と
これをGIS的にもパロディできないものか。
ということで、QGIS上でレイヤをシュレッダー（裁断）するプログラムを作ったというお話です。

といってもやってることは単純で、下の図のように、シュレッダーにかけたいレイヤを覆うようにいくつかの細長のポリゴンデータを作ります。

その細長ポリゴンフィーチャを順番に選択してはクリップ、次を選択してクリップ・・・・という具合です。

実際にPyqgisで書くと以下のようになりました。
(python3で記述しています）

実行するときは、最初の3行の処理設定だけ書き換えてコピーしてQGISのPythonコンソールにコピーしていただければシュレッダーにかけられます
（少し注意点がありますので実行の前に一番下まで読んでください）

#--------------以下を設定してください-------------------
#シュレッドするレイヤ名は？
layername=u''

#何分割?
shrednum =10

#縦(0)or横(1)
horizon = 1

#------------------------------------------------------


#レイヤをシュレッダーにかけるプログラム
import os
import glob
import processing



#入力レイヤ
layer = QgsProject.instance().mapLayersByName(layername)[0]
#layerin = iface.activeLayer()

(Dir,filename) = os.path.split(layer.dataProvider().dataSourceUri())
os.chdir(Dir)



#getExtent
box = layer.extent()
#print('xMinimum={0:.9f}'.format(box.xMinimum()))
#print('xMaximum={0:.9f}'.format(box.xMaximum()))
#print('yMinimum={0:.9f}'.format(box.yMinimum()))
#print('yMaximum={0:.9f}'.format(box.yMaximum()))
#print('width={0:.9f}'.format(box.width()))
#print('height={0:.9f}'.format(box.height()))
xmin = box.xMinimum()
xmax = box.xMaximum()
ymin = box.yMinimum()
ymax = box.yMaximum()

exta = str(xmin) + ',' + str(xmax) + ',' + str(ymin) + ',' +str(ymax)


#VectorGrid
if horizon == 0:
	shredGrid = processing.runAndLoadResults("qgis:creategrid",{'TYPE':2,'EXTENT':exta,'HSPACING':box.width()/shrednum,'VSPACING':box.height(),'HOVERLAY':0,'VOVERLAY':0,'CRS':layer.crs().authid(),'OUTPUT':'memory:'})
else:
	shredGrid = processing.runAndLoadResults("qgis:creategrid",{'TYPE':2,'EXTENT':exta,'HSPACING':box.width(),'VSPACING':box.height()/shrednum,'HOVERLAY':0,'VOVERLAY':0,'CRS':layer.crs().authid(),'OUTPUT':'memory:'})


shredlayer =  QgsProject.instance().mapLayersByName("グリッド")[0]




if layer.type() == 0:
	
	#Vector---------------------	
	
	#shred
	for i in range(1,shrednum+1,1):
		shredlayer.selectByExpression('"id" =' + str(i),QgsVectorLayer.SetSelection)
		
		
		processing.runAndLoadResults("native:clip", {"INPUT": layer, "OVERLAY":QgsProcessingFeatureSourceDefinition(shredGrid['OUTPUT'],True) , "OUTPUT": Dir + '/shred_'+layer.name() + "_" +  str(i) + '.shp'})
		
	
else:
	
	#Raster-------------------------------
	
	#shred
	for i in range(1,shrednum+1,1):
		shredlayer.selectByExpression('"id" =' + str(i),QgsVectorLayer.SetSelection)
		
		
		rasClip =processing.runAndLoadResults("gdal:cliprasterbymasklayer", {"INPUT": layer, "MASK":QgsProcessingFeatureSourceDefinition(shredGrid['OUTPUT'],True) ,"NODATA":None,"ALPHA_BAND":False, "CROP_TO_CUTLINE":True,"KEEP_RESOLUTION":True, "OPTIONS":'',"DATA_TYPE":0,"OUTPUT":Dir + '/shred_'+ layer.name() + "_" + str(i) + '.tif'})
	
#----------------------------------

#remove layer
QgsProject.instance().removeMapLayer(layer)
QgsProject.instance().removeMapLayer(shredlayer)



#shp一式選択
files = glob.glob(filename.split('shp')[0]+'*')


#Delete
for i in files:
	os.remove(i)

今日のAdvent Calendarではその内部処理をどう記述しているのかを説明しがてら、下の東京都周辺の鉄道線レイヤ（出典：国土地理院　基盤地図情報）にシュレッダーをかけたいと思います。
（自分も含め）これからpyQGIS扱おうって方にとって、「こんな時どう記述するの？」って時に役立てれば光栄です。

※えらそうなこと言いましたがpython初心者なので効率的でない処理をしてるかもしれません。気がついた点があればお知らせいただけると幸いです。

#処理設定
まずはじめの数行は入力レイヤと分割の設定です。
・1つ目はシュレッダにかけたいレイヤ名
・2つ目は何分割に裁断したいかを数字で
・3つ目は縦方向に裁断するか横方向に裁断するかを1or0で指定。

#--------------以下を設定してください-------------------

#シュレッドするレイヤは？
layername='layername'

#何分割?
shrednum =10

#縦(0)or横(1)
horizon = 1

続いて必要なモジュールのインポート

import os
(入力レイヤの保存先にアクセスするため）

import glob
(入力レイヤの一式のファイルを削除するため)

import processing
(QGISのプロセッシングを実行するため)

#入力レイヤの設定

#上で入力した名前のレイヤをレイヤインスタンスとして設定
layer = QgsProject.instance().mapLayersByName(layername)[0]

ちなみにレイヤ名を指定せずに、レイヤウインドウ上で選択しているレイヤを指定する場合は
layer = iface.activeLayer()
でいけるそうですね。

#入力レイヤの東西南北の範囲(extent)を取得

box = layer.extent()
xmin = box.xMinimum()
xmax = box.xMaximum()
ymin = box.yMinimum()
ymax = box.yMaximum()

#クリップ用のベクタグリッドを作成

ベクタ＞調査ツール＞グリッドの作成を使います。

通常メッシュデータなどを作成する時に使用するツールですが、これの縦幅を入力レイヤの南北の長さにするとシュレッダーをかけたみたいな短冊状のポリゴンができるというわけです。
コードで書くと以下のとおり

processing.runAndLoadResults("qgis:creategrid",{'TYPE':2,'EXTENT':exta,'HSPACING':box.width()/shrednum,'VSPACING':box.height(),'HOVERLAY':0,'VOVERLAY':0,'CRS':layer.crs().authid(),'OUTPUT':'memory:'})

ベクタグリッドの詳しい入力プロパティはalgorithmHelpしていていただくと幸いですが、こだわったポイントとしては以下2点。

・EXTENT: 作成するグリッドの東西南北の領域。先のextent処理で引っ張ってきた入力レイヤの領域が入るようにしています。
・HSPACING（水平方向間隔）とVSPACING（垂直方向間隔）
　-　最初の処理設定で入力した間隔の分だけ分割されるように水平or垂直方向の間隔が決定されるように、「入力レイヤの東西（もしくは南北）範囲の距離/分割数」にしています。

図でいうと以下のような関係

ちなみにこの短冊ポリゴン、FOSS4GのTLでお話させてもらった時に聞いてくださった方の一人から
「横方向にもできるの？」と聞かれたのでできるようにしました。
最初の設定でhorizonを0にすれば縦方向に、1にすれば横方向のポリゴンが作成されます。

（横の短冊ポリゴンの例）

#短冊レイヤのフィーチャを1つずつ選択

for i in range(1,shrednum+1,1):
	shredlayer.selectByExpression('"id" =' + str(i),QgsVectorLayer.SetSelection)

属性検索でforループをつかってフィーチャを順番に選択します。

あとは断裁したい入力レイヤにたいして選択したフィーチャで順番にクリップしていくだけになりますが・・・

#裁断したいレイヤがベクターならベクタクリップ、ラスターならラスタクリップ（当たり前　笑）

if layer.type() == 0:

これでレイヤのデータタイプが判別できます。0ならベクタ、1ならラスタになります

if文で入力レイヤがベクタならベクタクリップ

processing.runAndLoadResults("native:clip", {"INPUT": layer, "OVERLAY":QgsProcessingFeatureSourceDefinition(shredGrid['OUTPUT'],True) , "OUTPUT": Dir + '/shred_'+layer.name() + "_" +  str(i) + '.shp'})
		

ラスタクリップならラスタクリップ

processing.runAndLoadResults("gdal:cliprasterbymasklayer", {"INPUT": layer, "MASK":QgsProcessingFeatureSourceDefinition(shredGrid['OUTPUT'],True) ,"NODATA":None,"ALPHA_BAND":False, "CROP_TO_CUTLINE":True,"KEEP_RESOLUTION":True, "OPTIONS":'',"DATA_TYPE":0,"OUTPUT":Dir + '/shred_'+ layer.name() + "_" + str(i) + '.tif'})

が実行されるようになっております。
処理結果は以下のようになります

レイヤ一覧にクリップされたレイヤがたくさんあるのと、鉄道線が短冊状に色わけされていることがわかるかと思います。

もちろん一部を非表示にするとその部分だけ表示されなくなることからも裁断されていることがわかります。

ラスタもご覧のとおり

#これだけでは終わりません
作りたいのはシュレッダーにかけるプログラム。
シュレッダーなので不可逆です。
つまりもとのデータは帰ってきません

ということで、基データを削除するように記述しています。

QgsProject.instance().removeMapLayer(layer)
QgsProject.instance().removeMapLayer(shredlayer)

まず、↑でシュレッダーかける前の入力レイヤとシュレッダー用の短冊レイヤをレイヤ一覧から削除します。

files = glob.glob(filename.split('shp')[0]+'*')

#Delete
for i in files:
	os.remove(i)

次に入力レイヤの関連ファイル一式をコンピュータ上からも削除します
実行していただけるとわかるかと思いますが、基ファイルがディレクトリから完全に削除されます（ゴミ箱にも残りません）
実行の際は自己責任でお願いします。

#こんなときにお使いください
まぁ完全に自己満足なプログラムなわけで需要なんてないかもしれませんが、せめてこんなときにお使いください。

・レイヤをカラフルに見せたいとき
・「このレイヤは自動的に削除される」的な指令を出すとき
・データをオークションにだすとき

などなど
使ってみた方はご感想をお聞かせいただけると幸いです

以上、お粗末さまでした