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

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

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

## はじめに

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

なお、M座標というのはシェープファイルに格納できる4次元目の座標で、[ArcGIS Pro のヘルプページ](https://pro.arcgis.com/ja/pro-app/latest/help/data/geodatabases/overview/feature-class-basics.htm)には以下のように説明があります。

> ライン フィーチャの頂点には、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出力します。

## コード

```py
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本のラインです。

```py:サンプルデータ作成用
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 の各ラインの座標値を取り出します。

```py:実行
path_shp = 'sample.shp'
field = 'name'

output_xyzm(path_shp, field)
```

## 結果

`GISソフトでsample.shpの頂点の座標値を表示`
A.png
![A.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1409923/bdbf8b95-a110-d676-05f1-7a0b2c7391f6.png)

B.png
![B.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1409923/f0addfe3-37cd-5651-c044-8a2e1279b3ec.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 で分岐等が必要です。