第03回空間結合¶

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

この授業の進め方¶

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

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

学習目標¶

空間結合（spatial join）の概念を説明できる
R-treeインデックスの仕組みを理解し、なぜ高速化できるかを説明できる
ogr2ogr のSQL機能を使ってポイント-in-ポリゴン結合を実行できる
GeoPandasの sjoin() を使って空間結合を実行できる
空間結合の結果を集計してCSVに出力できる

1. 空間結合とは¶

1.1 通常の結合と空間結合の違い¶

通常のデータベース結合（JOIN）では、共通のキー（例：区コード、IDなど）を使って2つのテーブルを結合します。

通常のJOIN: テーブルA.district_code == テーブルB.district_code

空間結合（Spatial Join）では、地理的な位置関係をキーとして使います。

空間JOIN: ポイントAがポリゴンBの「内側にある」かどうか
          ポリゴンAがポリゴンBと「交差する」かどうか
          ポイントAがポイントBに「最も近い」かどうか

1.2 主な空間関係（述語）¶

空間関係	説明	使用例
`ST_Within` / `ST_Contains`	AがBの内側にある	地点が市区町村の中にある
`ST_Intersects`	AとBが交差または重なる	道路がバッファ範囲と交差する
`ST_DWithin`	AとBの距離がD以内	駅から500m以内の施設
`ST_Touches`	AとBが境界で接触する	隣接する行政区域

1.3 ポイント-in-ポリゴン検索¶

今回のメイン課題はポイント-in-ポリゴン（PiP）です。

問い: 「このGPS地点は、どの市区町村に属するか？」

入力:
  ・ポイントデータ  → LLMログの問い合わせ地点（緯度・経度）
  ・ポリゴンデータ  → 行政区域の境界（札幌市の区など）

出力:
  ・各ポイントに、そのポイントを含むポリゴンの属性（区名など）が付与される

図解（テキスト）:

  ┌─────────────────────────┐
  │    中央区               │
  │         ●(問い合わせ1)  │
  │  ●(問い合わせ2)        │
  ├─────────────────────────┤
  │    豊平区               │
  │              ●(問い合わせ3) │
  └─────────────────────────┘

結果: 問い合わせ1 → 中央区
      問い合わせ2 → 中央区
      問い合わせ3 → 豊平区

2. 空間インデックスの仕組み（R-tree）¶

2.1 インデックスがない場合の問題¶

インデックスなしで「このポイントがどのポリゴンに属するか」を調べると、すべてのポリゴンとの判定が必要です。

ポイント数:   N = 5,000
ポリゴン数:   M = 10,000
全判定回数:   N × M = 50,000,000 回（= 5千万回）

これでは時間がかかりすぎます。

2.2 R-tree（矩形木）インデックスの仕組み¶

R-treeは、空間データを効率的に検索するためのデータ構造です。各ポリゴンを最小境界矩形（MBR: Minimum Bounding Rectangle）で囲み、それを階層的にグルーピングします。

R-treeのイメージ（テキスト図解）:

レベル0（ルート）:
  ┌──────────────────────────────────────┐
  │          北海道全域のMBR             │
  └──────────────────────────────────────┘

レベル1（中間ノード）:
  ┌──────────────┐ ┌──────────────┐
  │  道南MBR      │ │  道北MBR     │
  └──────────────┘ └──────────────┘

レベル2（リーフ）:
  ┌────┐ ┌────┐ ┌────┐ ┌────┐
  │ 函館区 │ │ 室蘭区 │ │ 旭川区 │ │ 稚内区 │
  └────┘ └────┘ └────┘ └────┘

検索手順:

問い合わせ地点を含まないMBRは枝ごとスキップする
含む可能性のある枝だけを再帰的にたどる
最終的に候補ポリゴンだけで精密判定

比較:
  インデックスなし: 10,000 ポリゴン全部を判定
  R-treeあり:    数個〜数十個の候補だけを判定
  → 100倍以上高速になることも

2.3 R-treeを自動で使う方法¶

GDALのSQLite方言やGeoPandasは、空間インデックスを自動的に活用します。明示的に「インデックスを作れ」と指示しなくても、内部で使われています。

3. ogr2ogrによる空間結合¶

3.1 SQLiteダイアレクトを使った空間結合¶

ogr2ogr の -dialect SQLite オプションを使うと、SQLiteのSpatialite機能（ST_Within, ST_Intersects など）が使えます。

# ポイント-in-ポリゴン結合の基本形
# 各ポイントに、それを含む区の名前を付与する

ogr2ogr \
  -f "CSV" \
  data/processed/llm_log_with_district.csv \
  data/raw/llm_log.geojson \
  -dialect SQLite \
  -sql "
    SELECT
      p.id,
      p.timestamp,
      p.category,
      d.district_name
    FROM llm_log p
    LEFT JOIN 'data/raw/districts.geojson'.districts d
      ON ST_Within(p.geometry, d.geometry)
  "

オプションの説明:

オプション	内容
`-f "CSV"`	出力フォーマット
`-dialect SQLite`	SQLiteのSpatialite拡張を使用
`-sql "..."`	実行するSQL文
`LEFT JOIN`	マッチしないポイント（区域外）もNULLとして保持

3.2 区ごとの集計を同時に行う¶

# 空間結合と集計を1つのSQLで行う
ogr2ogr \
  -f "CSV" \
  results/count_by_district.csv \
  data/raw/llm_log.geojson \
  -dialect SQLite \
  -sql "
    SELECT
      d.district_name,
      COUNT(p.id) AS point_count
    FROM llm_log p
    LEFT JOIN 'data/raw/districts.geojson'.districts d
      ON ST_Within(p.geometry, d.geometry)
    GROUP BY d.district_name
    ORDER BY point_count DESC
  "

4. Pythonによる空間結合（GeoPandas）¶

4.1 GeoPandasのsjoin()の基本¶

import geopandas as gpd
import pandas as pd

# データの読み込み
gdf_points = gpd.read_file("data/raw/llm_log.geojson")
gdf_districts = gpd.read_file("data/raw/districts.geojson")

# CRSを確認（両方が同じCRSである必要がある）
print("ポイントのCRS:", gdf_points.crs)
print("ポリゴンのCRS:", gdf_districts.crs)

# CRSが異なる場合は合わせる
if gdf_points.crs != gdf_districts.crs:
    gdf_districts = gdf_districts.to_crs(gdf_points.crs)

# sjoin()による空間結合
# how="left": 左（ポイント）を全件保持、マッチしない場合はNaN
# predicate="within": ポイントがポリゴンの内側にある

gdf_joined = gpd.sjoin(
    left_df=gdf_points,       # ポイントデータ（左側）
    right_df=gdf_districts,   # ポリゴンデータ（右側）
    how="left",               # 結合方法（left / right / inner）
    predicate="within"        # 空間関係（within / intersects / contains）
)

print(f"結合前のポイント数: {len(gdf_points)}")
print(f"結合後の行数: {len(gdf_joined)}")
print(gdf_joined[["id", "category", "district_name"]].head(10))

4.2 sjoin()の結合オプション¶

# inner join: 空間的にマッチしたもののみを保持（区域外のポイントは除外）
gdf_inner = gpd.sjoin(gdf_points, gdf_districts, how="inner", predicate="within")

# left join: ポイントをすべて保持（ポリゴン外のポイントはdistrictがNaN）
gdf_left = gpd.sjoin(gdf_points, gdf_districts, how="left", predicate="within")

# intersects: より緩い条件（境界上の点も含む）
gdf_intersects = gpd.sjoin(gdf_points, gdf_districts, how="left", predicate="intersects")

4.3 区ごとの集計¶

import geopandas as gpd
import pandas as pd

# データ読み込み
gdf_points = gpd.read_file("data/raw/llm_log.geojson")
gdf_districts = gpd.read_file("data/raw/districts.geojson")

# CRS統一
gdf_districts = gdf_districts.to_crs(gdf_points.crs)

# 空間結合
gdf_joined = gpd.sjoin(gdf_points, gdf_districts, how="left", predicate="within")

# 区ごとのポイント数を集計
count_by_district = (
    gdf_joined
    .groupby("district_name")["id"]
    .count()
    .reset_index()
    .rename(columns={"id": "point_count"})
    .sort_values("point_count", ascending=False)
)

print(count_by_district)

# CSV出力（encoding="utf-8-sig" はExcelでも文字化けしない設定）
count_by_district.to_csv("results/count_by_district.csv", index=False, encoding="utf-8-sig")
print("results/count_by_district.csv に保存しました")

5. ハンズオン手順¶

所要時間の目安: 全手順で60〜90分程度です。授業時間内に終わらなければ持ち帰ってください。どこまで終わったかをメモしておくと再開しやすいです。

手順1: データの準備と確認¶

cd ~/bigdata2/pipeline_practice

# ポイントデータの確認
ogrinfo -al -so data/raw/llm_log.geojson

# ポリゴンデータの確認
ogrinfo -al -so data/raw/districts.geojson

確認すること：

両方のデータのCRS（EPSGコード）を確認する
CRSが同じかどうか確認する（違う場合は次の手順で統一する）

手順2: 両データのCRSを統一する¶

# ポリゴンをポイントのCRS（WGS84）に合わせる
ogr2ogr \
  -f "GeoJSON" \
  -t_srs EPSG:4326 \
  data/processed/districts_wgs84.geojson \
  data/raw/districts.geojson

# CRSが一致したことを確認
ogrinfo -al -so data/processed/districts_wgs84.geojson

うまくいかないとき: CRSが元から同じ場合はこの手順はスキップしてOKです。

手順3: ogr2ogrで空間結合を実行する¶

# 各ポイントにどの区かを付与する
ogr2ogr \
  -f "CSV" \
  data/processed/llm_log_with_district.csv \
  data/raw/llm_log.geojson \
  -dialect SQLite \
  -sql "
    SELECT
      p.id,
      p.timestamp,
      p.category,
      d.district_name
    FROM llm_log p
    LEFT JOIN 'data/processed/districts_wgs84.geojson'.districts d
      ON ST_Within(p.geometry, d.geometry)
  "

# 結果を確認
head -10 data/processed/llm_log_with_district.csv

うまくいかないとき（SQLエラー）: SQL query failed が出た場合、ファイルパスやレイヤー名が正しいか確認してください。ファイルパスの 'data/processed/...' は実際のパスに合わせてください。

手順4: GeoPandasで空間結合と集計を実行する¶

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

# src/spatial_join.py

import geopandas as gpd
import pandas as pd
import os

# ディレクトリ作成
os.makedirs("results", exist_ok=True)

# データ読み込み
print("データを読み込んでいます...")
gdf_points = gpd.read_file("data/raw/llm_log.geojson")
gdf_districts = gpd.read_file("data/raw/districts.geojson")

print(f"  ポイント数: {len(gdf_points)}")
print(f"  ポリゴン数: {len(gdf_districts)}")

# CRS統一
gdf_districts = gdf_districts.to_crs(gdf_points.crs)
print(f"  CRS: {gdf_points.crs}")

# 空間結合
print("空間結合を実行しています...")
gdf_joined = gpd.sjoin(gdf_points, gdf_districts, how="left", predicate="within")
print(f"  結合後の行数: {len(gdf_joined)}")

# 区ごとの集計
result = (
    gdf_joined
    .groupby("district_name", dropna=False)["id"]
    .count()
    .reset_index()
    .rename(columns={"id": "point_count"})
    .sort_values("point_count", ascending=False)
)

print("\n--- 区ごとのポイント数 ---")
print(result.to_string(index=False))

# CSV出力
output_path = "results/count_by_district.csv"
result.to_csv(output_path, index=False, encoding="utf-8-sig")
print(f"\n結果を保存しました: {output_path}")

python3 src/spatial_join.py

手順5: 結果の確認¶

cat results/count_by_district.csv

6. 演習課題（空間演習①）¶

この課題は評価対象です。

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

課題内容¶

学内配布の llm_log.geojson（ポイントデータ）と districts.geojson（行政区域ポリゴン）を使い、以下の処理を実装・実行してください。

要件:

GeoJSONのポイントデータと行政区域ポリゴンを空間結合する
区ごとのポイント数を集計する
結果を results/count_by_district.csv としてCSVファイルに出力する
CSVには最低限、以下の列を含めること：
district_name（区名）
point_count（ポイント数）
ポイント数の多い順に並べること

発展課題（任意・加点対象）:

カテゴリ別のクロス集計（results/cross_district_category.csv）も作成する
ポリゴンに集計結果を結合して、results/districts_with_count.geojson を作成する

# 発展課題のヒント：集計結果をポリゴンに戻す
gdf_districts_count = gdf_districts.merge(
    count_by_district,
    on="district_name",
    how="left"
)
gdf_districts_count["point_count"] = gdf_districts_count["point_count"].fillna(0).astype(int)
gdf_districts_count.to_file("results/districts_with_count.geojson", driver="GeoJSON")

提出物¶

以下を1つのフォルダにまとめてzipで提出してください。

src/spatial_join.py（実装したPythonスクリプト）
results/count_by_district.csv（出力されたCSVファイル）
README.md（手順の説明・躓いた点・気づいたこと）

README.mdには、完全に動かなかった場合でも、どこまでできたか・どこで詰まったかを書いてください。途中経過の提出も評価します。

まとめ¶

項目	内容
空間結合	位置関係（内包・交差など）をキーとしてデータを結合する
PiP（ポイント-in-ポリゴン）	地点がどのエリアに属するかを特定する
R-tree インデックス	空間インデックスの仕組み。MBRの階層で検索を高速化
ST_Within	ポイントがポリゴンの内側にある場合に真となるSQL関数
ogr2ogr -dialect SQLite	SQLiteのSpatialite拡張でGDALに空間SQL機能を追加
gpd.sjoin()	GeoPandasの空間結合関数。predicate で空間関係を指定
how="left"	左側（ポイント）の全件を保持する結合方式

次回（第04回）は、バッファ分析と最近傍分析により、地点間の距離・到達圏を計算します。

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

第03回 空間結合¶