2
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

シェープファイルの頂点のXYZM座標を取り出す

Last updated at Posted at 2021-05-21

対象読者 シェープファイルを扱う方, pythonユーザ

実行環境 Windows10, Python 3.9.2, pandas 1.2.3, pyshp 2.1.3

はじめに

シングルフィーチャのシェープファイルから、各オブジェクトの頂点の座標を取り出します。座標はXY座標だけでなく、Z座標とM座標も対象とします。

なお、M座標というのはシェープファイルに格納できる4次元目の座標で、ArcGIS Pro のヘルプページには以下のように説明があります。

ライン フィーチャの頂点には、M 値 (メジャー値) を追加することもできます。一部の GIS アプリケーションでは、道路、河川、パイプラインなどの線形フィーチャに沿った距離を補完するために、線形の計測値を使用します。フィーチャの各頂点に M 値を割り当てることができます。

また、ジオメトリの読み取りには pyshp を使います。

シェープファイルのジオメトリタイプは shapefile.Reader クラスの shapeTypeName または shapeType にアクセスして調べます。
今回は以下のラインデータを対象にします。

shapeTypeName shapeType
POLYLINE 3
POLYLINEZ 13
POLYLINEM 23

Z座標とM座標の両方を持っている場合は POLYLINE Z タイプになります。M座標の設定がない頂点のM座標は None 扱いです。
POLYLINE M にはZ座標がありません。

オブジェクト名も取得したいので、ジオメトリとレコードの両方にアクセスできる shapefile.ShapeRecords クラスを利用します。

処理としては、シンプルに shapefile.ShapeRecords クラスを for 文で回して各オブジェクトの座標値とオブジェクト名を取得し、オブジェクト別に座標値をcsv出力します。

コード

import pandas as pd
import shapefile

def output_xyzm(path_shp: str, field: str) -> None:
	# ファイルオープン
	with shapefile.Reader(path_shp) as r:
		# 全オブジェクトの読み込み
		shapeRecs = r.shapeRecords()
		# 各オブジェクトの読み取り
		for shapeRec in shapeRecs:
			# オブジェクト名
			name = str(shapeRec.record[field])
			# 辞書 {軸名: 座標値のリスト}
			dic_point = {}
			# (x, y) のリスト
			list_point_xy = shapeRec.shape.points
			# x のリスト
			dic_point['x'] = [xy[0] for xy in list_point_xy]
			# y のリスト
			dic_point['y'] = [xy[1] for xy in list_point_xy]
			try:
				# z のリスト
				dic_point['z'] = shapeRec.shape.z
			except AttributeError:
				pass
			try:
				# m のリスト
				dic_point['m'] = shapeRec.shape.m
			except AttributeError:
				pass
			# 出力用のpd.DataFrame
			df = pd.DataFrame(dic_point)
			# 出力
			df.to_csv(name + '.csv', index=False)

実行

適当なサンプルデータが見つからなかったので、こちらのコードでラインデータを自作しました。
5つの頂点がXYZM座標をもつAB2本のラインです。

サンプルデータ作成用
import shapefile

list_name = ['A', 'B']
with shapefile.Writer('sample.shp') as w:
	w.field('name', 'C')
	coord = 0
	for name in list_name:
		list_line = []
		for i in range(5):
			list_point = [coord + j for j in range(4) ]
			list_line.append(list_point)
			coord += 4
		w.linez([list_line])
		w.record(name)

この sample.shp の各ラインの座標値を取り出します。

実行
path_shp = 'sample.shp'
field = 'name'

output_xyzm(path_shp, field)

結果

GISソフトでsample.shpの頂点の座標値を表示
A.png
A.PNG

B.png
B.PNG

出力csvファイル(スペースで列幅を調整しています)

A.csv
   x,    y,    z,    m
 0.0,  1.0,  2.0,  3.0
 4.0,  5.0,  6.0,  7.0
 8.0,  9.0, 10.0, 11.0
12.0, 13.0, 14.0, 15.0
16.0, 17.0, 18.0, 19.0
B.csv
   x,    y,    z,    m
20.0, 21.0, 22.0, 23.0
24.0, 25.0, 26.0, 27.0
28.0, 29.0, 30.0, 31.0
32.0, 33.0, 34.0, 35.0
36.0, 37.0, 38.0, 39.0

同じ座標値であることが確認できました。

最後に

ラインタイプのシェープファイルに対してXYZM座標を取得しました。

単に頂点座標を取得するだけのコードですが、頂点間のZ座標やM座標を内挿で計算したい場合にも応用が利くと思います。

注意点

マルチフィーチャになってる場合や、ポイントタイプ、ポリゴンタイプは除外しています。
シェープファイルの読み込みで UnicodeDecodeError が出る場合は、shapefile.Reader() の引数に encoding を指定してください。
同じオブジェクト名が複数ある場合はファイル名が重複してしまうので、os.path.isfile で分岐等が必要です。

2
2
1

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
2
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?