第04回バッファ・最近傍分析¶

科目：ビッグデータ応用Ⅱ〔特徴量エンジニアリング〕
担当：小野原彩香
対象：3年生（後期）

この授業の進め方¶

この授業は個人作業・自己ペース型です。以下を前提として設計しています。

ハンズオン手順は授業時間内に全部終わらせる必要はありません。続きは持ち帰って自分のペースで進めてください。
演習課題の締め切りは2週間後です。
グループワークや口頭発表はありません。全作業は個人で行います。
質問はメールで受け付けています（返信は数日以内）。
各手順の後に「うまくいかないとき」の対処法を記載しています。
躓いたままにせず、途中経過でも提出してください。どこまでできたかを書いてくれれば評価します。

学習目標¶

バッファ分析の概念を説明し、ogr2ogr で実装できる
最近傍分析の概念を説明し、2地点間の距離を計算できる
等時圏（徒歩・自転車・車）の概念を理解し、バッファで近似できる
GeoPandasとmatplotlibを使ってバッファ・最近傍の結果を可視化できる
施設データからバッファを作成し、ポイントとの重複を検出できる
発展：「度」で作ったバッファが円にならない理由を、自分で導ける

1. バッファ分析とは¶

1.1 バッファとは¶

バッファ（Buffer）とは、ある地物（点・線・面）から指定した距離内の領域を示すポリゴンです。

バッファの例（テキスト図解）:

      　　　　　━━━━━━━━━━━━━
      　　　　╱ ← バッファ領域 ╲
      　　   ╱      ┌──────┐    ╲
      　　  │       │ 施設 │     │
      　　   ╲      └──────┘    ╱
      　　　　╲               ╱
      　　　　　━━━━━━━━━━━━━
            ←  半径 500m  →

1.2 バッファの用途¶

用途	説明
到達圏分析	施設から徒歩X分、半径Xmに到達できる範囲の推定
影響圏分析	騒音・汚染源からX m 以内の住民数を推計
近接性の特徴量化	地点から最寄り施設まで何mか（機械学習の特徴量）
重複検出	バッファ内に含まれる別のデータを抽出する

1.3 バッファ計算における座標系の重要性¶

バッファを正確な距離（メートル）で作るには、投影座標系（単位: m）を使う必要があります。

WGS84（EPSG:4326、単位: 度）でのバッファ:
  ST_Buffer(geometry, 0.005)  ← 「0.005度」の円。距離として意味不明

平面直角座標系（EPSG:6680、単位: m）でのバッファ:
  ST_Buffer(geometry, 500)    ← 500メートルの円。正確

必ず平面直角座標系に変換してからバッファを作成してください。

導出①：「度」で作ったバッファは、円ではなく南北に長い楕円になる

上の ST_Buffer(geometry, 0.005) が「距離として意味不明」なのは、単位が度だからというだけではありません。そもそも円になっていないことが問題です。

第02回で見たとおり、緯度1度と経度1度の長さは違います。札幌（\(\varphi = 43.062096°\)）では経度1度が緯度1度の \(\cos\varphi\) 倍ほどに縮みます。半径を「0.005度」と指定すると、南北と東西で実際の長さが変わります。

向き	0.005度の実際の長さ
南北	555.5 m
東西	407.3 m

東西は南北の 73.3% しかありません。できあがる図形は、南北に長い楕円です。

逆向きに考えるともっとはっきりします。「500mのバッファ」を度で作ろうとすると、

\[\text{南北}: \frac{500}{111{,}093.9} = 0.004501°, \qquad \text{東西}: \frac{500}{81{,}458.8} = 0.006138°\]

と2つの違う値が必要になり、ST_Buffer に渡す1つの半径では表せません。投影座標系（単位 m）に変換してからバッファを作れば、この問題は起きません。

使った数学：第02回の「経度1度は \(\cos\varphi\) で縮む」。割り算による単位の読み替え。

2. 最近傍分析とは¶

2.1 最近傍分析の概念¶

最近傍分析（Nearest Neighbor Analysis）は、ある地点からみて最も近い地物を見つける分析です。

問い: 「各ユーザーの問い合わせ地点から、最も近いバス停はどこか？」

ポイント（問い合わせ地点）:  ●  ●  ●  ●  ●
バス停:                        ▲      ▲      ▲

結果: 各●に、最も近い▲とその距離が付与される

2.2 距離の種類¶

距離の種類	説明	本講義での対応
ユークリッド距離	直線距離（「カラスが飛ぶ距離」）	`ST_Distance()`, `.distance()` で計算
マンハッタン距離	格子状の道路網での距離	今回は扱わない
ネットワーク距離	実際の道路/経路に沿った距離	今回は扱わない

本講義ではユークリッド距離（直線距離）を使います。実際の移動距離より短くなりますが、近接性の特徴量としては十分です。

導出②：平面の直線距離で近似してよいのは、範囲が狭いときだけ

ユークリッド距離は平面上の直線距離です。地球は曲がっているので、これは近似になります。どのくらいずれるのかを確かめます。

札幌（\(\varphi = 43.062096°\)）を起点に東へ動いたとき、緯線に沿って測った長さと、球面上の最短経路（大円距離）を比べます。

東へ	緯線に沿った長さ	大円距離	差
0.005°	407.3 m	406.2 m	1.09 m（0.27%）
0.05°	4,072.9 m	4,062.0 m	10.91 m（0.27%）
0.5°	40,729.4 m	40,620.3 m	109.11 m（0.27%）

差の割合はこの範囲ではほぼ一定で、0.3%に届きません。500mのバッファなら1m程度のずれです。

比べるべきは、導出①で見た 27% のゆがみのほうです。座標系を間違えたときの誤差は、曲面を平面で近似したときの誤差より桁違いに大きい。だから「まず投影座標系に直す」が先で、「直線距離で近似する」のはそのあとの話になります。

どの投影法を選ぶか、測地線（楕円体上の最短経路）を厳密に計算するとどうなるか、といった話はこの授業で扱う範囲外です。

使った数学：距離の比較と割合。大円距離の式そのものは使うだけで導きません。

2.3 等時圏（イソクロン）¶

等時圏（Isochrone）は、ある地点から一定時間以内に到達できる範囲です。簡易的にはバッファで近似できます。

簡易的な等時圏の計算（直線距離で近似）:

  徒歩10分:
    平均歩行速度 ≈ 80 m/分
    到達圏の半径 ≈ 80 × 10 = 800 m

  自転車10分:
    平均速度 ≈ 250 m/分
    到達圏の半径 ≈ 250 × 10 = 2,500 m

  車10分:
    平均速度 ≈ 500 m/分（市街地）
    到達圏の半径 ≈ 500 × 10 = 5,000 m

注意: 実際の等時圏は道路ネットワークを考慮するため、円形にはなりません。バッファはあくまで近似です。

3. ogr2ogrによるバッファ生成¶

3.1 基本的なバッファ生成¶

ステップ1: WGS84 → 平面直角座標系XII系（メートル単位）に変換する

ogr2ogr \
  -f "GeoJSON" \
  -t_srs EPSG:6680 \
  data/processed/facilities_6680.geojson \
  data/raw/facilities.geojson

ステップ2: 500mバッファを生成する

ogr2ogr \
  -f "GeoJSON" \
  data/processed/facilities_buffer_500m.geojson \
  data/processed/facilities_6680.geojson \
  -dialect SQLite \
  -sql "
    SELECT
      id,
      name,
      category,
      ST_Buffer(geometry, 500) AS geometry
    FROM facilities_6680
  "

# 確認（Extentがバッファ分だけ広がっているはず）
ogrinfo -al -so data/processed/facilities_buffer_500m.geojson

3.2 複数距離のバッファを一度に生成する¶

# 徒歩5分（400m）・10分（800m）・15分（1200m）の3種類のバッファ

ogr2ogr \
  -f "GeoJSON" \
  data/processed/facilities_buffer_multi.geojson \
  data/processed/facilities_6680.geojson \
  -dialect SQLite \
  -sql "
    SELECT id, name, '5min' AS walk_time, ST_Buffer(geometry, 400) AS geometry
    FROM facilities_6680
    UNION ALL
    SELECT id, name, '10min' AS walk_time, ST_Buffer(geometry, 800) AS geometry
    FROM facilities_6680
    UNION ALL
    SELECT id, name, '15min' AS walk_time, ST_Buffer(geometry, 1200) AS geometry
    FROM facilities_6680
  "

3.3 バッファとポイントの重複検出¶

# バッファ内に含まれるポイントを抽出する（ST_DWithin で距離条件）

ogr2ogr \
  -f "CSV" \
  results/points_in_buffer.csv \
  data/processed/llm_log_6680.geojson \
  -dialect SQLite \
  -sql "
    SELECT
      p.id,
      p.category,
      f.name AS facility_name,
      ROUND(ST_Distance(p.geometry, f.geometry), 1) AS distance_m
    FROM llm_log_6680 AS p
    JOIN 'data/processed/facilities_6680.geojson'.facilities_6680 AS f
      ON ST_DWithin(p.geometry, f.geometry, 500)
    ORDER BY f.name, distance_m
  "

# 結果の確認
head -20 results/points_in_buffer.csv

4. 距離・最近傍の計算¶

4.1 GeoPandasでの距離・最近傍計算¶

import geopandas as gpd
import pandas as pd

# データの読み込み（平面直角座標系で読む）
gdf_points = gpd.read_file("data/processed/llm_log_6680.geojson")
gdf_facilities = gpd.read_file("data/processed/facilities_6680.geojson")

print(f"ポイント数: {len(gdf_points)}")
print(f"施設数: {len(gdf_facilities)}")

# 各ポイントから最寄り施設までの距離を計算する関数
def get_nearest_distance(point, facilities_gdf):
    """
    1つのポイントから全施設への距離を計算し、最小値とその施設を返す
    """
    distances = facilities_gdf.geometry.distance(point)
    nearest_idx = distances.idxmin()
    return distances[nearest_idx], nearest_idx

# 全ポイントに対して最寄り距離を計算
print("最近傍距離を計算中...")
results = []
for idx, row in gdf_points.iterrows():
    dist, fac_idx = get_nearest_distance(row.geometry, gdf_facilities)
    results.append({
        "point_id": row["id"],
        "nearest_facility": gdf_facilities.loc[fac_idx, "name"],
        "distance_m": round(dist, 1)
    })

df_nearest = pd.DataFrame(results)
print(df_nearest.head(10))
print(f"\n平均距離: {df_nearest['distance_m'].mean():.1f} m")
print(f"最大距離: {df_nearest['distance_m'].max():.1f} m")
print(f"最小距離: {df_nearest['distance_m'].min():.1f} m")

うまくいかないとき（処理が遅い）: データ件数が多いと時間がかかります。まず gdf_points = gdf_points.head(100) で100件だけ試してから、全件で実行してください。

5. 可視化¶

5.1 バッファとポイントの基本的な可視化¶

以下のコードを src/visualize.py として保存し、実行してください。

# src/visualize.py

import geopandas as gpd
import matplotlib.pyplot as plt
import os

os.makedirs("results", exist_ok=True)

# データの読み込み
gdf_points = gpd.read_file("data/processed/llm_log_6680.geojson")
gdf_facilities = gpd.read_file("data/processed/facilities_6680.geojson")
gdf_buffer = gpd.read_file("data/processed/facilities_buffer_500m.geojson")
gdf_districts = gpd.read_file("data/raw/districts.geojson").to_crs(epsg=6680)

# 図の作成
fig, ax = plt.subplots(figsize=(12, 10))

# 1. 区ポリゴン（背景）
gdf_districts.plot(
    ax=ax,
    color="lightgray",
    edgecolor="white",
    linewidth=0.8,
    alpha=0.5
)

# 2. バッファポリゴン
gdf_buffer.plot(
    ax=ax,
    color="skyblue",
    edgecolor="steelblue",
    linewidth=1.0,
    alpha=0.4,
    label="施設バッファ（500m）"
)

# 3. 問い合わせポイント
gdf_points.plot(
    ax=ax,
    color="orange",
    markersize=3,
    alpha=0.6,
    label="問い合わせ地点"
)

# 4. 施設ポイント
gdf_facilities.plot(
    ax=ax,
    color="red",
    markersize=60,
    marker="*",
    zorder=5,
    label="施設"
)

# 凡例・タイトル
ax.legend(loc="upper right", fontsize=10)
ax.set_title("施設バッファ（500m）と問い合わせ地点の分布", fontsize=14, pad=15)
ax.set_xlabel("X座標（m）", fontsize=10)
ax.set_ylabel("Y座標（m）", fontsize=10)

plt.tight_layout()
plt.savefig("results/buffer_map.png", dpi=150, bbox_inches="tight")
print("地図を保存しました: results/buffer_map.png")
# plt.show()  # 画面に表示したい場合はコメントアウトを外す

python3 src/visualize.py

うまくいかないとき（日本語が文字化けする）: matplotlibでフォントが見つからない場合、以下を visualize.py の先頭に追加してください。
import matplotlib
matplotlib.rcParams['font.family'] = 'Hiragino Sans'  # macOSの場合
# Windowsの場合: matplotlib.rcParams['font.family'] = 'MS Gothic'

5.2 距離分布のヒストグラム¶

# src/plot_distance.py

import matplotlib.pyplot as plt
import pandas as pd
import os

os.makedirs("results", exist_ok=True)

# 最近傍距離のデータ（前節で作成済み）
df_nearest = pd.read_csv("results/nearest_facility.csv")

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# ヒストグラム
axes[0].hist(
    df_nearest["distance_m"],
    bins=30,
    color="steelblue",
    edgecolor="white",
    alpha=0.8
)
axes[0].axvline(
    df_nearest["distance_m"].mean(),
    color="red",
    linestyle="--",
    linewidth=2,
    label=f"平均: {df_nearest['distance_m'].mean():.0f}m"
)
axes[0].axvline(
    500,
    color="orange",
    linestyle="--",
    linewidth=2,
    label="バッファ半径: 500m"
)
axes[0].set_title("最寄り施設までの距離分布", fontsize=13)
axes[0].set_xlabel("距離 (m)", fontsize=11)
axes[0].set_ylabel("頻度", fontsize=11)
axes[0].legend(fontsize=10)

# 累積分布
sorted_dist = df_nearest["distance_m"].sort_values().values
cumulative = [i / len(sorted_dist) for i in range(len(sorted_dist))]

axes[1].plot(sorted_dist, cumulative, color="steelblue", linewidth=2)
axes[1].axhline(0.5, color="gray", linestyle=":", label="50%")
axes[1].axhline(0.8, color="gray", linestyle="-.", label="80%")
axes[1].axvline(500, color="orange", linestyle="--", label="500m")
axes[1].set_title("最寄り施設までの距離の累積分布", fontsize=13)
axes[1].set_xlabel("距離 (m)", fontsize=11)
axes[1].set_ylabel("累積割合", fontsize=11)
axes[1].legend(fontsize=10)
axes[1].set_ylim(0, 1)

plt.tight_layout()
plt.savefig("results/distance_distribution.png", dpi=150, bbox_inches="tight")
print("グラフを保存しました: results/distance_distribution.png")

6. ハンズオン手順¶

所要時間の目安: 全手順で60〜90分程度です。授業時間内に終わらなければ持ち帰ってください。

手順1: 座標系を確認してから変換する¶

cd ~/bigdata2/pipeline_practice

# 施設データのCRS確認
ogrinfo -al -so data/raw/facilities.geojson

# 施設データを平面直角座標系に変換
ogr2ogr \
  -f "GeoJSON" \
  -t_srs EPSG:6680 \
  data/processed/facilities_6680.geojson \
  data/raw/facilities.geojson

# LLMログも変換
ogr2ogr \
  -f "GeoJSON" \
  -t_srs EPSG:6680 \
  data/processed/llm_log_6680.geojson \
  data/raw/llm_log.geojson

手順2: 500mバッファを生成する¶

ogr2ogr \
  -f "GeoJSON" \
  data/processed/facilities_buffer_500m.geojson \
  data/processed/facilities_6680.geojson \
  -dialect SQLite \
  -sql "
    SELECT
      id,
      name,
      category,
      ST_Buffer(geometry, 500) AS geometry
    FROM facilities_6680
  "

# バッファが生成されたことを確認
ogrinfo -al -so data/processed/facilities_buffer_500m.geojson

手順3: バッファ内のポイントを検出する¶

ogr2ogr \
  -f "CSV" \
  results/points_in_buffer.csv \
  data/processed/llm_log_6680.geojson \
  -dialect SQLite \
  -sql "
    SELECT
      p.id,
      p.category,
      f.name AS facility_name,
      ROUND(ST_Distance(p.geometry, f.geometry), 1) AS distance_m
    FROM llm_log_6680 AS p
    JOIN 'data/processed/facilities_6680.geojson'.facilities_6680 AS f
      ON ST_DWithin(p.geometry, f.geometry, 500)
    ORDER BY f.name, distance_m
  "

# 件数確認
wc -l results/points_in_buffer.csv
head -10 results/points_in_buffer.csv

手順4: 最近傍距離を計算してCSV出力する¶

以下のスクリプトを src/nearest_analysis.py として保存し、実行してください。

# src/nearest_analysis.py

import geopandas as gpd
import pandas as pd
import os

os.makedirs("results", exist_ok=True)

# データ読み込み（平面直角座標系）
gdf_points = gpd.read_file("data/processed/llm_log_6680.geojson")
gdf_facilities = gpd.read_file("data/processed/facilities_6680.geojson")

print(f"ポイント数: {len(gdf_points)}")
print(f"施設数: {len(gdf_facilities)}")

# 最近傍計算
results = []
for idx, point_row in gdf_points.iterrows():
    distances = gdf_facilities.geometry.distance(point_row.geometry)
    nearest_idx = distances.idxmin()
    results.append({
        "point_id": point_row["id"],
        "category": point_row.get("category", ""),
        "nearest_facility_name": gdf_facilities.loc[nearest_idx, "name"],
        "distance_m": round(distances[nearest_idx], 1),
        "within_500m": distances[nearest_idx] <= 500
    })

df_result = pd.DataFrame(results)
df_result.to_csv("results/nearest_facility.csv", index=False, encoding="utf-8-sig")

print("\n--- 統計サマリー ---")
print(f"平均距離:          {df_result['distance_m'].mean():.1f} m")
print(f"中央値:            {df_result['distance_m'].median():.1f} m")
print(f"最大距離:          {df_result['distance_m'].max():.1f} m")
print(f"500m以内の件数:     {df_result['within_500m'].sum()} 件")
print(f"500m以内の割合:     {df_result['within_500m'].mean()*100:.1f} %")

python3 src/nearest_analysis.py

手順5: 結果を可視化する（余裕があれば）¶

python3 src/visualize.py
python3 src/plot_distance.py
ls results/

7. 演習課題¶

締め切り：2週間後（次の次の授業の前日23:59）
提出方法：作成したファイルをzipにまとめてメールで送付
再提出：締め切り後も1週間以内であれば再提出可能です。完成していなくても、どこまでできたかを書いて提出してください。

課題1（必須）：バッファの作成と重複検出¶

施設データを平面直角座標系（EPSG:6680）に変換する
各施設に半径500mのバッファを作成する（ogr2ogr または GeoPandas）
バッファ内に含まれるLLMログのポイントを抽出してCSVに保存する
施設ごとに「バッファ内のポイント数」を集計してCSVに保存する

課題2（必須）：最近傍分析¶

各LLMログのポイントから、最寄り施設までの距離（メートル）を計算する
以下の統計量を求めてREADMEに記録する：
平均距離・中央値・最大距離・最小距離
500m以内に施設があるポイントの割合

課題3（任意・加点対象）：可視化¶

バッファ・施設・LLMログを重ね合わせた地図を作成し、PNGファイルとして保存する。

提出物¶

src/buffer_analysis.py（バッファ・重複検出のスクリプト）
src/nearest_analysis.py（最近傍分析のスクリプト）
results/points_in_buffer.csv
results/nearest_facility.csv
README.md（実行手順、気づいた点、躓いた箇所）
（任意）results/buffer_map.png

まとめ¶

項目	内容
バッファ分析	地物から指定距離内の領域を作成する
等時圏（簡易）	徒歩/自転車/車の速度 × 時間 = 半径として円形バッファで近似
ST_Buffer(geom, r)	半径 r のバッファを生成するSQL関数
ST_DWithin(a, b, d)	a と b の距離が d 以内かどうかを判定するSQL関数
ST_Distance(a, b)	a と b の距離（m）を返すSQL関数
geometry.distance()	GeoPandas/Shapely で距離を計算するメソッド
座標系の注意	バッファ・距離の計算は必ず投影座標系（m単位）で行う
可視化	GeoDataFrame.plot() と matplotlib を組み合わせて地図を描く

これで第01〜04回のベースとなるGIS操作（再現性・CRS変換・空間結合・バッファ最近傍）をひととおり学びました。次回以降はこれらの技術を組み合わせ、機械学習用の特徴量マトリックスを構築していきます。

ビッグデータ応用Ⅱ〔特徴量エンジニアリング〕 / 担当：小野原彩香

第04回 バッファ・最近傍分析¶