《Python/GIS》国土数値情報のシェープファイルから施設の分布や個数を調べる

郵便局や鉄道駅のような公共施設の分布状況を可視化したり、個数を数えたりする方法を紹介します。

↓noteでも同様の記事を執筆しております。

https://note.com/my_channel7200/n/nb0e5f14f5fed

# 概要

国土交通省が提供している国土数値情報のデータを用います。

**国土交通省より、商用での利用は禁止されています。**

https://nlftp.mlit.go.jp/ksj/

下にスクロールし、調べたい施設をクリックしてください。
ここで、「ポイント」というのは施設を**点**で表示し、「ポリゴン」とは施設を**図形**で表示することを指します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3753744/02e2cfec-2f69-4475-9598-e51e5cff2295.png)

今回は例として、「郵便局」を扱います。クリックすると、ダウンロードするデータの対象を選択する画面が出てきます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3753744/a6ced2ab-b121-4551-b6a4-77f6e33dadb2.png)

本記事では、「全国」を例として扱います。

# Pythonを用いて郵便局の分布や個数を可視化
### データの中身
ダウンロードしたデータは、以下のような複数の拡張子のデータを持ちます。
・DBFファイル
・PRJファイル
・SHPファイル
・SHXファイル
これらは「**シェープファイル**」と呼ばれ、一つのファイルを削除するとうまく作動しないようになっています。

この中でも特に重要な.shpファイルを見てみましょう。私の場合は以下のようにシェープファイルを整理し、「公共施設個数取得コード.ipynb」というPythonファイルを作成しました。
（ご自身で整理される場合は、パスの部分を各自変更してください。）

＞作業用フォルダ
　＞公共施設個数取得コード.ipynb
　＞data
　　＞国土数値情報
　　　＞全国郵便局P30-13
　　　　＞P30-13
　　　　　＞P30-13.shp　など一式

公共施設個数取得コード.ipynbに、以下のようなコードを書きます。

``` python
import geopandas as gpd
import matplotlib.pyplot as plt

# 1. シェープファイルの読み込み
shapefile_path = "data/国土数値情報/全国郵便局P30-13/P30-13/P30-13.shp"  # 郵便局のシェープファイル
gdf = gpd.read_file(shapefile_path)

# 2. データの基本情報を確認
print("データの基本情報:")
print(gdf.head())  # 最初の5行を表示
print("\nカラム一覧:")
print(gdf.columns)  # カラム（列名）を確認

# 3. 施設の個数を確認
total_facilities = gdf.shape[0]
print(f"\n郵便局の総数: {total_facilities}件")


# 4. 施設分布の地図プロット
gdf.plot(figsize=(10, 10), color="blue", markersize=5, alpha=0.5)
plt.title("distribution of post office in Japan")
plt.xlabel("longtitude")
plt.ylabel("latitude")

plt.show()
```

2で、.shpファイルの中身を確認します。
3で、全国の郵便局数を確認します。（まだ市町村別にはできません）
4で、施設分布を地図にプロットします。

コードを実行すると、以下のように.shpファイルの中身を確認できます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3753744/481e052d-8b18-443b-9967-2b2b20f6e3be.png)

P30_001からP30_007までの列名と緯度経度を意味する列名などがあります。意味については、先ほどのダウンロード画面に書いてあります。

また、#4の実行結果として、日本の郵便局の分布が描画されます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3753744/65b1a4f1-d646-422c-b359-7a69fce339ba.png)

もう少し細かい分布を見たい場合は、QGISというソフトウェアと、行政区域データのシェープファイルを用いるのがおすすめです。

### 市町村別の郵便局数を取得
先ほどの.shpファイルの「**P30_001**」という列（1列目）は、5桁の市町村コードを意味しています。これをもとに、市町村名と照らし合わせます。

まずは、「P30_001」列でgroupbyすることで、各市町村コードのデータが何行存在するかを計算できます。
*groupbyは、( )内で指定した列に対して同じ値が登場する回数を返します。

``` python
#市町村別の施設の個数を確認
city_count = gdf.groupby("P30_001").size().reset_index(name="郵便局数")
city_count.shape
print(city_count)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3753744/c325caa1-f207-45c5-9e29-7ec0246dc55f.png)

これではまだどの市町村かわからないので、ダウンロード画面の「**行政区域コード**」をダウンロードします。これは、市町村コードと市町村名が一覧になったExcelファイルです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3753744/571dd5c5-cf29-4199-8edc-d3d983a552a7.png)

ダウンロードした「行政区域コード」の中身は以下のようになっています。今後pandasで作業を行う際、列名を取得したいので、上２行を削除します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3753744/7bd3e965-6009-4162-aef6-3e8756c365d4.png)

削除したらファイルを保存し、先ほどの**作業用フォルダの中**に入れます。

### 行政区域コードと郵便局シェープファイルを結合
今の行政区域コードと、先ほど扱った郵便局のシェープファイルは、ともに５桁の市町村コードの列を持っています。pandasのmergeを用いて、この列を基準に結合します。

１点注意点として、市町村コードは文字列 (str) 指定をしないと、先頭に0をもつコードは4桁の数字として扱われてしまいます。

不要な列は取り除き、必要な「都道府県名」「市区町村名」「郵便局数」の３列のみを抽出し、最後にExcelファイルに出力して完成です！！

``` python
import pandas as pd

# ファイルをDataFrameとして読み込み
#5桁の市町村コードを文字列として扱うため，dtype=strで読み込み（これをしないと先頭が0のコードは4桁と認識されてしまう）
city_code = pd.read_excel("AdminiBoundary_CD.xlsx",dtype={"改正後のコード": str})

city_count["P30_001"] = city_count["P30_001"].astype(str)


# 5桁の市町村コードを基準に結合する場合
merged_df = pd.merge(city_code, city_count, left_on="改正後のコード",right_on ="P30_001", how="left")

#不要な列は除き、「都道府県名」「市区町村名」と「郵便局数」の3列のみを抽出
selected_merged_df = merged_df.iloc[:, [1,2,12]]

#print(merged_df)
print(selected_merged_df)

# 結果をExcelファイルに出力
output_excel = "都道府県市区町村名と郵便局数.xlsx"
selected_merged_df.to_excel(output_excel, index=False, engine="openpyxl")
```

↓完成した「都道府県市区町村名と郵便局数.xlsx」の中身
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3753744/96a6bfae-28eb-4ca6-8919-3c90296d757b.png)
ちなみに政令指定都市は、各区ごとに郵便局数が与えられているため、市の行は空欄になっています。

### Pythonコードまとめ
``` python
from mip import * # 数理最適化モデリング言語Python-MIP
import shapefile # shapeファイルを扱うためのモジュール

import matplotlib.pyplot as plt
%matplotlib inline
import pandas as pd
import json
import itertools
import folium
import time
import os

import geopandas as gpd
import matplotlib.pyplot as plt

# 1. シェープファイルの読み込み
shapefile_path = "data/国土数値情報/全国郵便局P30-13/P30-13/P30-13.shp"  # 郵便局のシェープファイル
gdf = gpd.read_file(shapefile_path)

# 2. データの基本情報を確認
print("データの基本情報:")
print(gdf.head())  # 最初の5行を表示
print("\nカラム一覧:")
print(gdf.columns)  # カラム（列名）を確認

# 3. 施設の個数を確認
total_facilities = gdf.shape[0]
print(f"\n郵便局の総数: {total_facilities}件")


# 4. 施設分布の地図プロット
gdf.plot(figsize=(10, 10), color="blue", markersize=5, alpha=0.5)
plt.title("distribution of post office in Japan")
plt.xlabel("longtitude")
plt.ylabel("latitude")
plt.show()
```
``` python
#市町村別の施設の個数を確認
city_count = gdf.groupby("P30_001").size().reset_index(name="郵便局数")
city_count.shape
print(city_count)

# ファイルをDataFrameとして読み込み
#5桁の市町村コードを文字列として扱うため，dtype=strで読み込み（これをしないと先頭が0のコードは4桁と認識されてしまう）
city_code = pd.read_excel("data/AdminiBoundary_CD.xlsx",dtype={"改正後のコード": str})

city_count["P30_001"] = city_count["P30_001"].astype(str)


# 5桁の市町村コードを基準に結合する場合
merged_df = pd.merge(city_code, city_count, left_on="改正後のコード",right_on ="P30_001", how="left")

#不要な列は除き、「都道府県名」「市区町村名」と「郵便局数」の3列のみを抽出
selected_merged_df = merged_df.iloc[:, [1,2,12]]

#print(merged_df)
print(selected_merged_df)

# 結果をExcelファイルに出力
output_excel = "都道府県市区町村名と郵便局数.xlsx"
selected_merged_df.to_excel(output_excel, index=False, engine="openpyxl")
```


今回の記事は以上になります。この記事が少しでもご参考になったと思って頂けたら、いいねやフォローをしていただけるととても力になります。