GIS に欠かせない Shapefile を深堀ってみる

# 1. はじめに
地理情報システム (GIS) で、使用するデータとして「Shapefile 形式」でデータ提供されていることがあります。この記事では、この Shapefile が、どのような構造となっているのかについて、説明していきます。

## 1.1. トポロジ構造とは
Shapefile について説明していく前に、まずは GIS データ等でよく利用される「トポロジ構造」について説明しておこうと思います。

例えば、トポロジ構造を持つファイルとして、次のようなものがあります。

- **.gdb**: ESRIのファイルジオデータベースは、トポロジを含む複雑な空間データを格納することができます。
- **.topojson**: トポロジを保存するための拡張 GeoJSON 形式で、隣接する地物間の境界を共有することでファイルサイズを削減します。
- **.gml**: 地理マークアップ言語 (Geography Markup Language) は、地理的特徴とそのトポロジを表現するために使用されます。

これらのフォーマットは、地理的データのトポロジ関係（例えば、接続性、隣接性、包含関係など）を保持し、クエリや分析においてこの情報を活用するために使用されます。

トポロジ構造は、GIS や CAD などで使われるデータ構造の一つで、地図上の点、線、ポリゴンなどの地理的要素間の空間的関係を表現します。トポロジ構造は、地理的要素間の接続性、隣接性、面認識などを正確に管理することができるため、多くの地理情報システムで重要な役割を果たしています。

- **接続性**:
    接続性は、地理的要素間の接続関係を表します。これは主に線状の要素に関連しており、例えば道路ネットワーク内での道路同士の接続や、河川システムにおける河川の流れのつながりを示すのに使用されます。
    接続性により、GIS はネットワーク分析を行うことができ、最短経路の計算、到達可能性の分析、ネットワーク内の障害点の影響評価などが可能になります。
    <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/57ae99ab-e72d-def5-87e0-ba0ce6208921.png" width="480px" />

- **面認識**:
    面認識は、地理的要素が形成する空間的領域を認識し、定義する能力を指します。これは、ポリゴン（多角形）を使用して、特定の地域や空間が占める面積を正確に表現することに関連しています。
    面認識を通じて、GIS は土地の区画、行政区域、自然保護区などの面積を計算し、その境界を管理することができます。また、面認識は、地理的要素がどのように空間を占めているかを理解し、空間分析や地理的クエリの実行において重要な役割を果たします。
    <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/bb39d703-a756-2615-463e-6fd307c62cd3.png" width="450px" />

- **隣接性**:
    隣接性は、地理的要素間の隣接関係、つまり互いに隣り合っている関係を指します。
    これは主に面状の要素（土地利用区域、湖、国境など）に適用されます。
    隣接性により、GIS はどの地理的要素が互いに接しているかを認識し、例えば、特定の土地利用区分がどの区分に囲まれているか、または特定の国がどの国と国境を共有しているかを分析することができます。
    <img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/1f7c37e2-347f-55f0-1449-0e4f8bbd54ea.png" width="480px" />


### 1.1.1. メリット
トポロジ構造には以下のようなメリットがあります。

- **データ整合性の保持**: トポロジルールを使用して地理的要素間の関係を厳格に管理することで、データの整合性を保持しやすくなります。例えば、境界線が正確に接続されていることを保証することができます。
- **空間分析の精度向上**: 地理的要素間の正確な関係を利用して、より複雑で高度な空間分析を行うことが可能になります。例えば、ネットワーク分析やトポロジ的なクエリがこの例にあたります。
- **データ圧縮**: 同じ境界線や頂点を共有する地理的要素を一度だけ格納することで、データの重複を避け、データストレージの効率を向上させることができます。

## 1.1.2. デメリット
一方で、トポロジ構造は以下のようなデメリットもあります。

- **処理の複雑さ**: トポロジ構造を構築し、維持するには複雑な処理が必要であり、システムの実装が困難になる場合があります。
- **初期設定の手間**: トポロジ構造を使用するためには、地理的要素間の正確なトポロジ関係を初期に設定する必要があり、これが手間と時間を要する作業となることがあります。
- **更新時のコスト**: 地理的要素の追加や削除があるたびに、トポロジ関係を再計算し、更新する必要があります。これにより、データのメンテナンスコストが上がる可能性があります。

# 2. Shapefile とは

Shapefile は、ESRI によって開発された地理情報のためのポピュラーなファイル形式で、点、線、ポリゴンなどの地理的要素を保存するために広く使用されています。

Shapefile 自体はトポロジ構造ではなく、重ねて表示をしたり、連続性を持たないといった特性を持ちます。そのため、処理がシンプル、かつサイズが軽量となっており、トポロジ構造を持つデータと比較して、以下のようなメリットがあります。

- 描画速度が速い
- ファイルサイズが小さい
- 読み取りおよび書き込みが容易

## 2.1. Shapefile のファイル種別
Shapefile は、複数のファイルで構成されるデータセットです。最低限、.shp ファイル、.shx ファイル、.dbf ファイルの３つのファイルが必要となります。

:::note info
Shapefile は、`map_file.shp`、`map_file.shx`、`map_file.dbf` というように、ファイル名を揃える必要があります。
:::

### 2.1.1. .shp ファイル
各ジオメトリの実際の点、線、ポリゴンデータを格納します。

ファイルは、ファイルヘッダとレコードの集合体で構成されます。ヘッダには、ファイルのバージョン、シェープタイプ、およびデータの空間的範囲（バウンディングボックス）が含まれます。レコードは、個々のジオメトリデータを含みます。

| shapeType | shapeType 名 | 説明 |
|:--|:--|:--|
| 1 | POINT<br/>(ポイント) | 単一の地理的位置を表す。地図上の１つの位置を示し、例えば住所や特定の地点を指すのに使用される。
| 3 | POLYLINE<br/>(ポリライン) | ２つ以上の点が連結されて形成される線を表す。道路、川などの線形の地理的特徴を示すのに使われる。
| 5 | POLYGON<br/>(ポリゴン) | ３つ以上の点が連結され、最初の点と最後の点が接続されることで閉じた形状を形成する。湖、国境などの面的な地理的特徴を表すのに使用される。
| 8 | MULTIPOINT<br/>(マルチポイント) | 複数の点が１つの地理的要素として扱われる場合を表す。複数の小さな島や建物の群れなど、互いに関連する点の集まりを示すのに使われる。
| 11 | POINTZ<br/>(ポイント Z) | POINT 形式に Z 座標 (高度や深度) が付加されたもの。 |
| 13 | POLYLINEZ<br/>(ポリライン Z) | POLYLINE 形式に Z 座標 (高度や深度) が付加されたもの。 |
| 15 | POLYGONZ<br/>(ポリゴン Z) | POLYGON 形式に Z 座標 (高度や深度) が付加されたもの。 |
| 18 | MULTIPOINTZ<br/>(マルチポイント Z) | MULTIPOINT 形式に Z 座標 (高度や深度) が付加されたもの。 |
| 21 | POINTM<br/>(ポイント M) | POINT 形式に M 値 (距離や時間などの測定値) が付加されたもの。 |
| 23 | POLYLINEM<br/>(ポリライン M) | POLYLINE 形式に M 値 (距離や時間などの測定値) が付加されたもの。 |
| 25 | POLYGONM<br/>(ポリゴン M) | POLYGON 形式に M 値 (距離や時間などの測定値) が付加されたもの。 |
| 28 | MULTIPOINTM<br/>(マルチポイント M) | MULTIPOINT 形式に M 値 (距離や時間などの測定値) が付加されたもの。 |
| 31 | MULTIPATCH<br/>(マルチパッチ) | 複数の３次元形状が組み合わされて形成される複雑な形状を表す。通常、地形や他の３次元の地理的特徴を表現するのに使用され、ポリゴンやポリラインなどの異なる形状が組み合わさって構成されている。 |

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/55e011dc-a326-9c5b-5c29-965c1dedc373.png)

### 2.1.2. .shx ファイル
.shp ファイル内のレコードへのインデックスを提供します。このファイルを使って、特定のレコードを迅速に検索できます。

### 2.1.3. .dbf ファイル
属性データを格納します。各レコードは、.shpファイルのジオメトリに対応する属性情報を持ちます。フィールドは、名前、タイプ（文字、数値、日付など）、長さ、小数点以下の桁数などの属性を持ちます。

### 2.1.4. その他の補助ファイル

Shapefile は、他にも以下のような補助ファイルを含めることもあります。

- **.prj ファイル**: プロジェクション情報を格納するファイル。地理座標系や投影座標系を定義する。
- **.sbn ファイル**、**.sbx ファイル**: 空間インデックスファイル。GIS ソフトウェアがジオメトリデータをより効率的に処理できるようにする。
- **.cpg ファイル**: 文字エンコーディングを指定するファイル。これにより、属性データの正しい読み込みと表示が保証される。

### 2.1.5. .shp ファイル、.shx ファイル、.dbf ファイルの関係性
Shapefile は地理的な要素を格納するための人気のある非トポロジカルな形式となっており、通常は複数のファイルにより構成されます。
下図は、地理情報システム（GIS）のコンテキストで使用される、Shapefileの3つの主要なファイルタイプ（.shp、.shx、.dbf）の関係を以下に示したものです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/fb87ae7b-6199-1d3e-ce06-9855874ca621.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/dd8d0350-8333-c64f-8de7-25adca0d81c5.png)

- **.shp ファイル**:
これは「メインファイル」とも呼ばれ、実際のジオメトリデータが含まれています。このファイルは地理的な特徴を定義する点、線、ポリゴンなどの形状データを保持します。画像の左の表では、**shapeType** と **points** という2つの列があります。**shapeType** は、形状の種類（例えばポイント、ライン、ポリゴンなど）を示し、「1」はポイントを意味します。**points** 列には、実際の緯度経度座標がリストされており、これが地理的な位置を定義します。

- **.shx ファイル**:
「インデックスファイル」と呼ばれ、.shp ファイルのレコードへのインデックスを提供します。これにより、プログラムは特定のジオメトリを迅速に検索することができます。画像の中央の表には、**Offset** と **Content Length** という列があります。**Offset** は .shp ファイル内のジオメトリの開始位置をバイト単位で指し、**Content Length** はそのジオメトリのレコードの長さを示しています。

- **.dbf ファイル**:
これは「属性データファイル」として知られ、dBase 形式でテーブルデータを格納します。各レコードは .shp ファイル内の形状に関連付けられた属性情報を持っています。画像の右の表には、**osm_id**、**code**、**fclass**、**population**、**name** といった列があり、これらは人口、名称、分類コードなど、それぞれの地理的な特徴に関連したデータを表します。

画像には、これら3つのファイルタイプ間の具体的な関連性が線で示されています。特定の点の座標は、特定の **Offset** と **Content Length** に関連付けられています。そして、これは特定の属性データにリンクされています。

この連携により、GIS アプリケーションは地理的な位置、形状、およびそれに関連する属性情報を効果的に処理し、例えば地図上での可視化や解析に利用することができます。

# 3. Python を使用した Shapefile の操作
## 3.1. Shapefile の読み込み
### 3.1.1. `pyshp` ライブラリの使用
`pyshp` ライブラリは、Python で Shapefile を扱うためのライブラリです。このライブラリを使うことにより、Shapefile 内でサポートされているすべての形状タイプを作成・読み込み・編集することができます。

以下のコマンドで、`pyshp` ライブラリをインストールします。

```bash
$ pip install pyshp
```

Python コードの冒頭で `pyshp` ライブラリをインポートします。

```python
import shapefile
```

`pyshp` ライブラリを使用して Shapefile ファイルを読み込みます。この際、`encoding` オプションをして、エンコードを指定することもできます。

今回は国土数値情報ダウンロードサイトの「[医療機関データ](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-P04-v2_1.html)」を使用しています。

```python
# Shapefileのパス
shp_path = 'P04-14_13_GML/P04-14_13-g_MedicalInstitution.shp'

# Shapefileの読み込み
sf = shapefile.Reader(shp_path, encoding='cp932')
```


### 3.1.2. ジオメトリデータの読み取り
以下のように記述することで、すべてのジオメトリデータを確認することができます。

ここでは、shape オブジェクト、shapeType、shape オブジェクトの位置情報を表示しています。
```python
print("Geometry:")
for shape in sf.shapes():
    print(shape, '/', shape.shapeType, '/', shape.points)
```

shapeType は `1` となっているので、POINT データであることがわかります（「[2.1.1. .shp ファイル](#211-shp-ファイル)」参照）。
```text
Geometry:
Shape #0: POINT / 1 / [[139.7654047, 35.69851345]]
Shape #1: POINT / 1 / [[139.7666533, 35.69695747]]
Shape #2: POINT / 1 / [[139.7630328, 35.69842448]]
...
```

### 3.1.3. 属性データの読み取り

属性データだけでは、何の値かわからないことも多々あるため、以下のソースコードでは、まず属性データのフィールド名を表示した後、各属性データの表示をしています。
```python
# 属性データのフィールド名の表示
print("\nFields:")
fields = sf.fields[1:]  # 最初のフィールドは削除フラグであるためスキップ
field_names = [field[0] for field in fields]
print(field_names)

# 属性データの表示
print("\nRecords:")
for record in sf.records():
    print(record)
```

```text
Fields:
['P04_001', 'P04_002', 'P04_003', 'P04_004', 'P04_005', 'P04_006', 'P04_007']

Records:
Record #0: ['1', '（医）済安堂井上眼科病院', '東京都千代田区神田駿河台４－３', '眼科\u3000麻酔科', '', '', '4']
Record #1: ['1', '神尾記念病院', '東京都千代田区神田淡路町２－２５', '精神科\u3000形成外科\u3000耳鼻いんこう科\u3000皮膚科\u3000美容皮膚科\u3000麻酔科\u3000美容外科', '', '', '4']
Record #2: ['1', '（公財）佐々木研究所附属杏雲堂病院', '東京都千代田区神田駿河台１－８', '内科\u3000呼吸器内科\u3000消化器内科\u3000肝臓内科\u3000循環器内科\u3000腫瘍内科\u3000リウマチ科\u3000呼吸器外科\u3000消化器外科\u3000乳腺外科\u3000整形外科\u3000形成外科\u3000婦人科\u3000泌尿器科\u3000放射線科\u3000美容皮膚科\u3000精神科\u3000麻酔科', '', '', '6']
...
```

フィールド名 (`Fields`) に表示されている項目の順で、属性データ (`Records`) がリスト形式で表示されています。

データの取得元である国土数値情報ダウンロードサイトの「[医療機関データ](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-P04-v2_1.html)」を確認すると以下のように属性名の定義が記述されており、これと照らし合わせることでどのような値か把握することができそうです。

| shp<br/>属性名　　 | 属性名　　　　 | 説明 |
|:--|:--|:--|
| P04_001 | 医療機関分類 | 病院、診療所、歯科診療所の区分 |
| P04_002 | 施設名称 | 医療施設の名称 |
| P04_003 | 所在地 | 医療機関の所在地の住所 |
| P04_004 | 診療科目１ | 当該施設が有する診療科目 |
| P04_005 | 診療科目２ | 当該施設が有する診療科目（「診療科目１」が254バイトを超える場合使用） |
| P04_006 | 診療科目３ | 当該施設が有する診療科目（「診療科目２」が254バイトを超える場合使用） |
| P04_007 | 開設者分類 | 当該施設（病院のみ）を開設した組織または個人の分類 |

## 3.2. Shapefile の作成
このセクションでは、`pyshp` ライブラリを使って、Shapefile を作成するソースコードを、よく使われる以下の shapeType を例に説明してきます。

- POINT (ポイント)
- POLYLINE (ポリライン)
- POLYGON (ポリゴン)
- MULTIPOINT (マルチポイント)


### 3.2.1. Shapefile の作成
Shapefile を作成するには、Shapefile Writer を定義し、属性やポイント情報などを書き込みます。

例えば、`example_null` という名前の空の Shapefile を作成する場合は、以下のように記述します。
今回は、`name` という文字列型の属性フィールドを定義し、属性データとして 'Example for Empty Shapefile' を挿入しています。

```python
import shapefile

with shapefile.Writer('example_null', shapeType=shapefile.NULL) as w:
    w.field('name', 'C')  # 'C' は文字列型を意味する
    # NULL Shapeを追加
    w.null()
    w.record('Example for Empty Shapefile')
```

実行すると .shp, .shx, .dbf ファイルが生成されます。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/31200e54-650f-5636-9004-d661f39358a0.png" width="320px" />

生成されたファイル確認のため、オープンソースの地理情報システムである [QGIS](https://www.qgis.org/ja/site/) を使って、example_null.shp を読み込んでみると定義通りデータが入っていることが確認できます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/582e3378-2aea-3bc8-8907-9c55d265cf5f.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/37abde76-862e-8987-ce85-e554b3357bc2.png)


### 3.2.1. POINT ジオメトリの作成
新宿区役所の場所を POINT ジオメトリとして描画する場合は、以下のように記述します。

```python
with shapefile.Writer('example_point', shapeType=shapefile.POINT) as w:
    w.field('name', 'C')
    w.point(139.70342179531394, 35.69389388730781)
    w.record('新宿区役所')
```

QGIS で開き、属性データも表示した例を以下に示します。
（わかりやすくするため、バックに OpenStreetMap レイヤを表示しています。）
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/6612500d-1b02-d187-d030-9a16583de983.png" width="500px"/>

属性データとして、数値を使うと、数値の範囲によって、色を変えるといったこともできます。

### 3.2.2. POLYLINE ジオメトリの作成
POLYLINE ジオメトリで、新宿アルタから新宿区役所までの移動経路を表示してみます。
ルート情報や河川、道路などを描画する際によくPOLYLINE ジオメトリが使用されます。

```python
with shapefile.Writer('example_polyline', shapeType=shapefile.POLYLINE) as w:
    w.field('name', 'C')
    # POLYLINE ジオメトリのリスト
    w.line([[[139.701102, 35.6925463], [139.70107, 35.69249], \
             [139.7010706, 35.6924882], [139.70096, 35.69253], \
             [139.70091, 35.69254], [139.7009067, 35.692545], \
             [139.70097, 35.69266], [139.70103, 35.69278], \
             [139.70109, 35.69287], [139.70115, 35.69298], \
             [139.70118, 35.69305], [139.70119, 35.69306], \
             [139.70123, 35.69314], [139.7013, 35.69326], \
             [139.70132, 35.69331], [139.70138, 35.69341], \
             [139.70139, 35.69342], [139.70144, 35.69341], \
             [139.70145, 35.69344], [139.70153, 35.69368], \
             [139.70154, 35.69371], [139.70157, 35.69371], \
             [139.7015685, 35.6937072], [139.70165, 35.69392], \
             [139.7017, 35.69406], [139.70171, 35.6941], \
             [139.70177, 35.69425], [139.70181, 35.69437], \
             [139.7018133, 35.6943694], [139.70213, 35.69429], \
             [139.7022, 35.69427], [139.70238, 35.69423], \
             [139.70247, 35.6942], [139.70252, 35.69419], \
             [139.7026, 35.69417], [139.70267, 35.69416], \
             [139.70279, 35.69412], [139.70282, 35.69412], \
             [139.70289, 35.6941], [139.70295, 35.69408], \
             [139.70301, 35.69407], [139.70313, 35.69404], \
             [139.7031278, 35.6940394]]])
    w.record('新宿アルタ～新宿区役所')
```

頂点をつなぐように POLYLINE ジオメトリが描画されていることが確認できます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/823cb4a9-5b30-b383-c2c0-b2a01d7a6e5e.png)

### 3.2.3. POLYGON ジオメトリの作成

POLYGON ジオメトリで、新宿アルタの区画を表示してみます。
```python
with shapefile.Writer('example_polygon', shapeType=shapefile.POLYGON) as w:
    w.field('name', 'C')
    # POLYGON ジオメトリのリスト
    w.poly([[[139.701047631901, 35.6925422790551], [139.701267572649, 35.6925019792727], \
             [139.701405706472, 35.6927219953152], [139.701482819958, 35.6926991222383], \
             [139.701543840553, 35.6927862572184], [139.701466056452, 35.692816210059], \
             [139.701511654378, 35.6929017113191], [139.70129573596, 35.6929507250553], \
             [139.701164307566, 35.6927317980543], [139.701129438806, 35.6927454128803]]])
    w.record('新宿アルタ')
```

これで、多角形の POLYLINE ジオメトリが表示されます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/26cd93c1-44e2-450c-a950-f5dfcd995ca4.png)

### 3.2.4. MULTIPOINT ジオメトリの作成
今回は、渋谷駅周辺のコンビニ（セブンイレブン、ファミリーマート、ローソン）の位置情報を MULTIPOINT ジオメトリで描画してみます。

```python
with shapefile.Writer('example_multipoint2', shapeType=shapefile.MULTIPOINT) as w:
    w.field('name', 'C')
    w.multipoint([[139.6996402462862, 35.65832798417496], \
                  [139.70642485473059, 35.66018986344715], \
                  [139.70271672168957, 35.657690264955704], \
                  [139.70630733154596, 35.65573297666057]])
    w.record('セブンイレブン')
    w.multipoint([[139.69883722149757, 35.65992599458355], \
                  [139.7020490995507, 35.66079835696538], \
                  [139.70509850422766, 35.657939464255264], \
                  [139.70746239320937, 35.65992593017978]])
    w.record('ファミリーマート')
    w.multipoint([[139.70115591664006, 35.65951552988615], \
                  [139.70371817774168, 35.659170995591275], \
                  [139.70490914834318, 35.66079832939078]])
    w.record('ローソン')
```

MULTIPOINT ジオメトリを使うとカテゴリ値による分類などで、種別ごとに色を変えたりすることもできるようになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/adf93362-05e1-dcba-0608-363dc5f8e74d.png)

# 4. 最後に
Shapefile の実態を知っておくと、GIS 等で Shapefile を活用する際に役立ちますし、`pyshp` と組み合わせることで、Python を使って、DB 上のデータや CSV ファイルなどから独自の Shapefile を生成したりといったことが容易に行えるようになります。

:::note info
今回は説明を割愛しておりますが、～M ジオメトリや ～Z ジオメトリは、それぞれ、地物情報の Z 値、M 値として情報が保存されるようです。
（Z 値は、想定通り認識させるに至りませんでした。）
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/2dab7f74-27bd-9e2c-2daf-1c3959bf2325.png" width="300px"/>
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2659159/f432f9d3-b223-277a-f751-17d3eba15fed.png" width="300px"/>
:::

# 参考資料
- https://saylordotorg.github.io/text_essentials-of-geographic-information-systems/s08-02-vector-data-models.html
- https://www.esri.com/content/dam/esrisites/sitecore-archive/Files/Pdfs/library/whitepapers/pdfs/shapefile.pdf
- https://biz-pasco.jp/help/wordlist/detail/000208
- https://desktop.arcgis.com/ja/arcmap/latest/manage-data/coverages/coverage-topology.htm
- https://pypi.org/project/pyshp/