日々適当

hibitekitou

Shapefileを読んでみる

cg |2016-02-20

Esri社が提唱したベクトル形式のフォーマットであるところのShapefileはこちら [ESRIジャパン] によると正しい表記はシェープファイル・Shapefileであり、拡張子はSHPだそうだが、こちら、国土交通省で配布している国土数値情報ダウンロードサービスから落とせるファイルフォーマットの一つだったりします。
で、行政区域の平成27年1月現在のデータを読み込んでみようと考えるわけです。

Shapefileを読み込もうとするとき、Pythonには幾つかモジュールがあるようです。んで、ググって真っ先に見つかったpyshp[GitHub]を使うことにしました。

import shapefile

でモジュールをインポートします。
そいでもってReaderオブジェクトを用意します。

sf = shapefile.Reader("ShapeFilePath")

ちなみに、今回想定しているのは、前述の国土数値情報ダウンロードサービスから落とせる行政区域データです。例えば広島県ならN03-15_34_150101.shp というファイルでした。なのでデスクトップにこのファイルが置いてあったら、

sf = shapefile.Reader("/Users/username/Desktop/N03-15_34_150101.shp")

でございます。実際は付随する幾つかのファイルがあるので、.shpだけデスクトップに直置きってことはないと思いますが。ここの書き方は幾つかあるようなので、ドキュメントをご覧になってください。

Readerオブジェクトを介して、ジオメトリを読み込んでいきます。
ジオメトリのリストはshapes()メソッドで取得できます。取得されるシェープレコードは幾つかの属性を持っているようで、今回の目的にはpointsとpartsを利用します。
pointsは位置情報のリスト。partsは後述します。
ちなみに、shapeTypeって属性もあって、これは汎用的にいろんなShapefileを読もうとするなら考慮すべき値となりましょう。ジオメトリがどんなタイプのものかを示す値が格納されています。今回読んでいる行政区域データの場合、基本的にPOLYGON型ということで5という値が返ってきます。これがPOINTだったりPOLYLINEだったらその他いろいろな型が規定されているようです。

というわけで、各々の行政区域を構成している境界線（閉曲線。つまりはポリゴン）を示す点群の情報を取得してみます。例えば10個目のやつについては

points = shapes[10].points

となります。シンプルですね。
しかして、ここで問題になるのが、ある区域に穴が空いている場合です。ドーナツ状と言いますか、そんな形状の場合です。pointsはそんなドーナツの穴を構成する点の情報もいっしょくたに羅列してしまっているのですな。なので例えば34213という名前の行政区画を読み込むとこんな形になってしまいます。

ちなみにこちらは廿日市市でございまして、隣の大竹市の飛び地が廿日市市から見ると穴となって見えていることになります。

というわけで、ドーナツの外側と内側を別個の線にするための指標となるのが parts っていう属性です。shapesの10個目のやつについては

parts = shapes[ 10 ].parts

となりますね。

廿日市市の穴が空いている区域のparts属性を見ると

[0, 9740, 9900, 9911, 9923, 9988, 10022, 10920, 10945, 11233, 11469, 11493]

となっています。ちなみに穴が空いていない区域は [0] と返ってきます。これは点群の何番目の要素で線が分離するかを示した値となります。つまり

0～9739
9740～9899
9900～9910
・
・
・
11493～11581

という範囲の値がそれぞれ区域を作成していることを示している、と。ちなみに最後の11581という値は、pointsの要素の数が11582であることから導いています。

ところで、先ほど、市町村コード34213とかいう数字を出しましたが、これはshapes()で取得されるジオメトリのレコードには入っていない情報です。それは

recs = sf.records()

で取得しています。このファイルの場合、各区域について4つの値が入っています。先の廿日市市の問題の区域は435番目の区域なので例えばそこについての値を見てみますと、

recs[435][0]：広島県
recs[435][1]：廿日市市
recs[435][2]：（空白）
recs[435][3]：（空白）
recs[435][4]：34213

ってなっています。～県～群～町とか～県～市～区とか、そんな風に最大4つの単位を入れらるようになっているようですね。
で、最後の数字（5番目の値）は市町村コード [総務省] のようでした。

というわけで、Softimageのnurbsカーブとして書くためのデータを書き出してみました。一行につき一つの閉曲線を構成する点のリストと市町村コードを記しています。

市町村コード,posx,posy,posz,1.0,posx,posy,posz,1.0....

ってなってます。意外と処理に時間がかかるので、以下は市町村コード34213だけを対象としていますw

# -*- coding: utf-8 -*-
import shapefile
import re

sf = shapefile.Reader("pathto/N03-15_34_150101.shp")
shapes = sf.shapes()
recs = sf.records()

forSIStr = ""

for i in range(len(shapes)):
	points = shapes[i].points
	parts = shapes[i].parts
	if recs[ i ][4] == '34213':
		for k in range( len( parts ) ):
			startNum = parts[k]
			if k < len( parts ) - 1:
				endNum = parts[k + 1]
			else:
				endNum = len( points ) 
			print("list from %d to %d"%(startNum,endNum))
			for j in range( startNum , endNum ):
				if j == startNum:
					forSIStr = forSIStr + recs[ i ][4]
				forSIStr = forSIStr + "," + str( points[ j ][0] ) + ",0.0," + str(-points[ j ][1]) + ",1.0"
				if j == endNum - 1:
					forSIStr = forSIStr + "\r\n"

f = open( "/Users/username/Desktop/Hiroshima34213.txt", "w")
f.write( forSIStr.encode('utf-8') ) 
f.close()

んで、書き出したテキストファイルをSoftimageで読み込んでみたのがこれ。