第10回 時空間データキューブ¶
科目：ビッグデータ応用Ⅱ〔特徴量エンジニアリング〕 担当：小野原彩香 対象：3年生（後期）
学習目標¶
この回を終えると、以下のことができるようになります。
- 時空間データキューブの3次元構造（x, y, time）を説明できる
- Polarsを使ってLLMログデータを時間ビンごとに集計できる
- xarrayでデータキューブを構築し、時系列スライスを取り出せる
- 折れ線グラフとヒートマップカレンダーで変化を可視化できる
- DuckDBの時空間SQLクエリを記述できる
1. 時空間データキューブとは¶
1.1 基本概念¶
通常の地理データは「空間（x, y）」の2次元です。 これに「時間（t）」を加えた3次元の構造を時空間データキューブと呼びます。
time →
t1 t2 t3 t4
y↑ ┌───────────────┐
│ v v v v │ ← ある(x,y)セルの時系列値
│ v v v v │
│ v v v v │
y↓ └───────────────┘
x →
各セル（x, y, t）には、そのグリッドセル・その時刻における集計値が入ります。
例：地理空間LLMログデータの場合 - x軸（経度）：東西方向のグリッドセルID - y軸（緯度）：南北方向のグリッドセルID - time軸：時間ビン（1時間単位、1日単位、1週単位など） - 値：クエリ数、エラー数、平均レイテンシ、ユニークユーザー数 など
1.2 なぜ時空間キューブが有用か¶
時空間キューブを使うと、以下のような問いに答えられます。
- 「ある地域のクエリ数は、曜日によって変化するか？」
- 「エラーが多発している時間帯と地域はどこか？」
- 「特定のイベント（祭り、災害）後に利用パターンはどう変わったか？」
これらは空間データだけでも時系列データだけでも答えられません。両方を組み合わせることで初めて分析できます。
1.3 時間ビンの種類¶
| ビン単位 | 用途 | DuckDB式 |
|---|---|---|
| 1時間 | 日中パターンの把握 | date_trunc('hour', ts) |
| 1日 | 曜日パターンの把握 | date_trunc('day', ts) |
| 1週 | 週次トレンド | date_trunc('week', ts) |
| 1月 | 季節性 | date_trunc('month', ts) |
2. データキューブの構築¶
2.1 LLMログデータの構造¶
data/llm_geo_logs.parquet
├── timestamp : タイムスタンプ（datetime）
├── longitude : 経度（float64）
├── latitude : 緯度（float64）
├── query_text : クエリテキスト（string）
├── response_time_ms: 応答時間（float64, ミリ秒）
├── is_error : エラーフラグ（bool）
├── error_type : エラー種別（string）
├── user_id : 匿名化ユーザーID（string）
└── grid_id : グリッドセルID（string, 例: "H36_V22"）
2.2 グリッドセルの定義¶
緯度経度を均等なグリッドに変換します。
import polars as pl
import numpy as np
# グリッドサイズ（度）
GRID_SIZE = 0.01 # ≈ 1km
def assign_grid(df: pl.DataFrame, lon_col: str, lat_col: str,
grid_size: float = 0.01) -> pl.DataFrame:
"""
緯度経度を指定グリッドサイズのセルIDに変換する
"""
return df.with_columns([
(pl.col(lon_col) / grid_size).floor().cast(pl.Int32).alias('grid_x'),
(pl.col(lat_col) / grid_size).floor().cast(pl.Int32).alias('grid_y'),
]).with_columns(
(pl.col('grid_x').cast(pl.Utf8) + '_' +
pl.col('grid_y').cast(pl.Utf8)).alias('grid_id')
)
3. 時系列集計と変化検出¶
3.1 Polarsによる時系列集計¶
Polarsはpandasより高速な列指向データフレームライブラリです。 特に大量のログデータ（数百万行）の集計に適しています。
import polars as pl
from datetime import datetime
# データ読み込み
df = pl.read_parquet('data/llm_geo_logs.parquet')
print(f"レコード数: {len(df):,}")
print(df.head())
# タイムスタンプを時間ビンに丸める
df = df.with_columns([
pl.col('timestamp').dt.truncate('1h').alias('hour_bin'),
pl.col('timestamp').dt.truncate('1d').alias('day_bin'),
pl.col('timestamp').dt.weekday().alias('weekday'), # 0=月曜
pl.col('timestamp').dt.hour().alias('hour_of_day'),
])
# グリッドセルIDを付与
df = assign_grid(df, 'longitude', 'latitude', GRID_SIZE)
3.2 時間×グリッドの集計¶
# 日×グリッドで集計（データキューブの断面を作る）
cube_daily = (
df
.group_by(['day_bin', 'grid_id', 'grid_x', 'grid_y'])
.agg([
pl.len().alias('query_count'),
pl.col('is_error').sum().alias('error_count'),
pl.col('is_error').mean().alias('error_rate'),
pl.col('response_time_ms').mean().alias('mean_response_ms'),
pl.col('response_time_ms').quantile(0.95).alias('p95_response_ms'),
pl.col('user_id').n_unique().alias('unique_users'),
])
.sort(['day_bin', 'grid_id'])
)
print("日次集計結果:")
print(cube_daily.head(10))
print(f"\n形状: {cube_daily.shape}")
3.3 変化検出：前日比の計算¶
# 特定のグリッドセルの時系列を取得
target_grid = 'H36_V22'
ts_target = (
cube_daily
.filter(pl.col('grid_id') == target_grid)
.sort('day_bin')
.with_columns([
pl.col('query_count').shift(1).alias('query_count_prev'),
pl.col('error_rate').shift(1).alias('error_rate_prev'),
])
.with_columns([
((pl.col('query_count') - pl.col('query_count_prev')) /
pl.col('query_count_prev') * 100).alias('query_change_pct'),
(pl.col('error_rate') - pl.col('error_rate_prev')).alias('error_rate_change'),
])
)
print("変化検出結果（最初の10日）:")
print(ts_target.select(['day_bin', 'query_count', 'query_change_pct',
'error_rate', 'error_rate_change']).head(10))
3.4 移動平均による平滑化¶
# 7日移動平均（ノイズ除去）
ts_smoothed = (
ts_target
.with_columns([
pl.col('query_count')
.rolling_mean(window_size=7, min_periods=3)
.alias('query_count_7d_ma'),
pl.col('error_rate')
.rolling_mean(window_size=7, min_periods=3)
.alias('error_rate_7d_ma'),
])
)
4. 可視化¶
4.1 折れ線グラフ（時系列トレンド）¶
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import pandas as pd
plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
# Polars → Pandasに変換して可視化
ts_pd = ts_smoothed.to_pandas()
fig, axes = plt.subplots(2, 1, figsize=(14, 8), sharex=True)
# 上段：クエリ数
axes[0].fill_between(ts_pd['day_bin'], ts_pd['query_count'],
alpha=0.3, color='steelblue', label='実測値')
axes[0].plot(ts_pd['day_bin'], ts_pd['query_count_7d_ma'],
color='navy', lw=2, label='7日移動平均')
axes[0].set_ylabel('クエリ数')
axes[0].set_title(f'グリッド {target_grid} の時系列変化')
axes[0].legend()
# 下段：エラー率
axes[1].fill_between(ts_pd['day_bin'], ts_pd['error_rate'] * 100,
alpha=0.3, color='crimson', label='エラー率（実測）')
axes[1].plot(ts_pd['day_bin'], ts_pd['error_rate_7d_ma'] * 100,
color='darkred', lw=2, label='7日移動平均')
axes[1].set_ylabel('エラー率（%）')
axes[1].set_xlabel('日付')
axes[1].legend()
# x軸の日付フォーマット設定
for ax in axes:
ax.xaxis.set_major_formatter(mdates.DateFormatter('%m/%d'))
ax.xaxis.set_major_locator(mdates.WeekdayLocator(interval=2))
plt.tight_layout()
plt.savefig('output/timeseries_target_grid.png', dpi=150, bbox_inches='tight')
plt.show()
4.2 ヒートマップカレンダー¶
# 時刻×曜日のヒートマップ
heatmap_data = (
df
.group_by(['weekday', 'hour_of_day'])
.agg([
pl.len().alias('query_count'),
pl.col('is_error').mean().alias('error_rate'),
])
.sort(['weekday', 'hour_of_day'])
.to_pandas()
)
# ピボットテーブル化
heatmap_pivot = heatmap_data.pivot(
index='weekday', columns='hour_of_day', values='query_count'
)
weekday_labels = ['月', '火', '水', '木', '金', '土', '日']
import seaborn as sns
fig, axes = plt.subplots(2, 1, figsize=(16, 8))
# クエリ数ヒートマップ
heatmap_error = heatmap_data.pivot(
index='weekday', columns='hour_of_day', values='error_rate'
)
sns.heatmap(heatmap_pivot, ax=axes[0], cmap='YlOrRd',
yticklabels=weekday_labels, linewidths=0.5,
cbar_kws={'label': 'クエリ数'})
axes[0].set_title('曜日×時間帯別 クエリ数ヒートマップ（全グリッド合計）')
axes[0].set_xlabel('時刻（0〜23時）')
axes[0].set_ylabel('曜日')
sns.heatmap(heatmap_error * 100, ax=axes[1], cmap='coolwarm',
yticklabels=weekday_labels, linewidths=0.5,
fmt='.1f', annot=True,
cbar_kws={'label': 'エラー率（%）'})
axes[1].set_title('曜日×時間帯別 エラー率ヒートマップ')
axes[1].set_xlabel('時刻（0〜23時）')
axes[1].set_ylabel('曜日')
plt.tight_layout()
plt.savefig('output/heatmap_calendar.png', dpi=150, bbox_inches='tight')
plt.show()
4.3 xarrayによるデータキューブ構築と空間スライス¶
import xarray as xr
# Polars → Pandas経由でxarrayに変換
cube_pd = cube_daily.to_pandas()
# グリッドの範囲を特定
grid_x_vals = sorted(cube_pd['grid_x'].unique())
grid_y_vals = sorted(cube_pd['grid_y'].unique())
day_vals = sorted(cube_pd['day_bin'].unique())
# xarrayデータセット構築
query_cube = cube_pd.pivot_table(
index='day_bin', columns=['grid_y', 'grid_x'],
values='query_count', fill_value=0
)
# xarray DataArray として整形
import itertools
ds = xr.Dataset()
for var in ['query_count', 'error_rate', 'mean_response_ms']:
pivot = cube_pd.pivot_table(
index='day_bin',
columns=['grid_y', 'grid_x'],
values=var,
fill_value=np.nan
)
# MultiIndex を展開
data_array = xr.DataArray(
pivot.values.reshape(
len(day_vals), len(grid_y_vals), len(grid_x_vals)
),
dims=['time', 'y', 'x'],
coords={
'time': day_vals,
'y': grid_y_vals,
'x': grid_x_vals,
},
name=var
)
ds[var] = data_array
print("データキューブの構造:")
print(ds)
# 特定の日のスナップショット取得
snapshot = ds['query_count'].sel(time='2024-03-15')
print(f"\n2024-03-15のスナップショット形状: {snapshot.shape}")
# 特定グリッドの時系列取得
ts_query = ds['query_count'].sel(x=136, y=35, method='nearest')
print(f"\n時系列の長さ: {len(ts_query)}")
5. ハンズオン手順¶
5.1 環境確認¶
# ライブラリのバージョン確認
python -c "import polars; print('Polars:', polars.__version__)"
python -c "import xarray; print('xarray:', xarray.__version__)"
python -c "import duckdb; print('DuckDB:', duckdb.__version__)"
5.2 DuckDBを使った時空間SQL集計¶
import duckdb
# DuckDBでインメモリ分析
conn = duckdb.connect()
# Parquetファイルを直接読み込む（コピー不要）
conn.execute("INSTALL spatial; LOAD spatial;")
# 時間帯別×グリッド別集計クエリ
result = conn.execute("""
SELECT
date_trunc('hour', timestamp) AS hour_bin,
floor(longitude / 0.01) AS grid_x,
floor(latitude / 0.01) AS grid_y,
COUNT(*) AS query_count,
SUM(is_error::INTEGER) AS error_count,
AVG(is_error::FLOAT) AS error_rate,
AVG(response_time_ms) AS mean_response_ms,
PERCENTILE_CONT(0.95) WITHIN GROUP
(ORDER BY response_time_ms) AS p95_response_ms
FROM read_parquet('data/llm_geo_logs.parquet')
WHERE timestamp >= '2024-01-01'
AND timestamp < '2024-04-01'
GROUP BY ALL
ORDER BY hour_bin, grid_x, grid_y
""").pl() # .pl() でPolars DataFrameとして取得
print("DuckDB集計結果:")
print(result.head(10))
print(f"\nレコード数: {len(result):,}")
5.3 空間フィルタリングクエリ¶
# 札幌都心部のみ抽出（ST_Contains使用）
result_sapporo = conn.execute("""
SELECT
date_trunc('day', timestamp) AS day_bin,
COUNT(*) AS query_count,
AVG(is_error::FLOAT) AS error_rate
FROM read_parquet('data/llm_geo_logs.parquet')
WHERE ST_Contains(
ST_GeomFromText('POLYGON((141.30 43.02, 141.42 43.02,
141.42 43.12, 141.30 43.12,
141.30 43.02))'),
ST_Point(longitude, latitude)
)
GROUP BY day_bin
ORDER BY day_bin
""").pl()
print("札幌都心部の日次集計:")
print(result_sapporo.head())
5.4 Parquetへの保存¶
# 集計結果をParquetで保存（次回以降の分析で再利用）
cube_daily.write_parquet('data/cube_daily.parquet')
result.write_parquet('data/cube_hourly.parquet')
print("保存完了:")
import os
for f in ['data/cube_daily.parquet', 'data/cube_hourly.parquet']:
size_mb = os.path.getsize(f) / 1024 / 1024
print(f" {f}: {size_mb:.1f} MB")
6. 演習課題¶
課題1（基本）¶
data/llm_geo_logs.parquetを読み込み、以下を実行してください。
- Polarsで「週×グリッドセル」の集計テーブルを作成する
- 最もクエリ数が多かったグリッドセルTOP10を表示する
課題2（応用）¶
- 曜日×時間帯のヒートマップを「エラー率」で作成する
- エラー率が特に高い時間帯・曜日を3つ特定し、その理由を推測してください（コメントとして記述）
課題3（発展）¶
- DuckDBの空間クエリ（
ST_Distance）を使って、「各グリッドセルから最寄りの鉄道駅までの距離」を計算する - 鉄道駅データ：
data/train_stations.parquet（longitude, latitude列あり） - 距離帯（0-500m, 500m-1km, 1km-5km, 5km以上）別のエラー率を比較するグラフを描く
提出形式：Jupyterノートブック、ファイル名：10_spatiotemporal_cube_[学籍番号].ipynb
まとめ¶
| 項目 | ポイント |
|---|---|
| 時空間キューブ | (x, y, time)の3次元構造で時空間パターンを表現 |
| Polarsの活用 | 遅延評価・列指向処理で大規模ログの高速集計 |
| 時間ビン | 分析目的に応じて時間・日・週単位を使い分ける |
| xarray | 多次元配列として扱い、スライス・インデックスが容易 |
| DuckDBの強み | Parquetを直接SQLで操作、空間関数も利用可能 |
| 可視化の選択 | 折れ線グラフ（トレンド）・ヒートマップ（パターン）を使い分ける |
| 次回への接続 | 今回作成した集計特徴量を第11回のk-meansクラスタリングに活用する |
次回予告： 第11回では、空間特徴量に基づいてグリッドセルを類型化（ゾーニング）するk-meansクラスタリングを実施します。
作成：小野原彩香 / ビッグデータ応用Ⅱ〔特徴量エンジニアリング〕