# Spatiotemporal Graph Feature Engineering

> Feature engineering patterns for spatiotemporal graph prediction tasks. Use when: (1) Working with time series on graph structures, (2) Predicting node-level outcomes over time, (3) Network dynamics with spatial dependencies, (4) Kaggle competitions involving traffic, flood, power grid, or social networks. Covers: neighbor aggregation, temporal windows, cross-domain features.

- Skill: `topprismdata/spatiotemporal-graph-feature-engineering` (Agent Skill)
- Install (CLI): `npx skillmds@latest add topprismdata/spatiotemporal-graph-feature-engineering`
- Raw SKILL.md: https://api.skillmd.com/api/skills/topprismdata/spatiotemporal-graph-feature-engineering/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: topprismdata (https://skillmd.com/u/topprismdata)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/topprismdata/spatiotemporal-graph-feature-engineering

---


# Spatiotemporal Graph Feature Engineering

## Problem
Predicting outcomes on spatiotemporal graphs (nodes + edges + time) requires
specialized features that capture both spatial relationships and temporal patterns.
Standard tabular or time-series features often miss critical graph structure.

## Context / Trigger Conditions
- Data has nodes, edges, and time steps (e.g., traffic sensors, hydraulic networks)
- Need to predict node-level values over time
- Graph structure influences predictions (neighbors matter)
- Common in: flood modelling, traffic prediction, power grid, social networks

## Solution

### Core Insight: Graph Features Trump Individual Features

**Critical discovery**:
```
Single-node features:
  elevation:     1179 importance
  degree:         345 importance

Neighbor aggregation:
  neigh_elev_mean: 1239 importance  ← MOST IMPORTANT!
```

**Lesson**: Aggregating information from neighbors consistently outperforms
single-node features in graph-structured prediction tasks.

### Feature Categories

#### 1. Graph Structure Features (Static)

**Node Degree**:
```python
from_counts = pd.concat([
    edges['from_node'],
    edges['to_node']
]).value_counts()

nodes['degree'] = nodes['node_idx'].map(from_counts).fillna(0).astype(int)
```

**Neighbor Statistics** (MOST IMPORTANT):
```python
from collections import defaultdict

# Build adjacency list
adj = defaultdict(list)
for _, row in edges.iterrows():
    adj[int(row['from_node'])].append(int(row['to_node']))

# Aggregate neighbor features
elev_map = dict(zip(nodes['node_idx'], nodes['elevation']))

neighbor_features = []
for n in nodes['node_idx']:
    neighbors = adj.get(n, [])
    if neighbors:
        # Multiple aggregation strategies
        neighbor_values = [elev_map.get(x, elev_map[n]) for x in neighbors]

        neighbor_features.append({
            'neigh_elev_mean': np.mean(neighbor_values),
            'neigh_elev_std': np.std(neighbor_values),
            'neigh_elev_min': np.min(neighbor_values),
            'neigh_elev_max': np.max(neighbor_values),
            'neigh_count': len(neighbors),
        })
    else:
        # Isolated nodes
        neighbor_features.append({
            'neigh_elev_mean': elev_map[n],
            'neigh_elev_std': 0,
            'neigh_elev_min': elev_map[n],
            'neigh_elev_max': elev_map[n],
            'neigh_count': 0,
        })
```

**Higher-Order Features**:
```python
# 2-hop neighbors
def get_k_hop_neighbors(adj, node, k=2):
    visited = set()
    current = {node}
    for _ in range(k):
        current = {n for curr in current for n in adj.get(curr, [])}
        visited.update(current)
    return visited

# Clustering coefficient
def clustering_coefficient(adj, node):
    neighbors = adj.get(node, [])
    if len(neighbors) < 2:
        return 0
    # Count edges between neighbors
    neighbor_pairs = 0
    for n1 in neighbors:
        for n2 in neighbors:
            if n2 in adj.get(n1, []):
                neighbor_pairs += 1
    return neighbor_pairs / (len(neighbors) * (len(neighbors) - 1))
```

#### 2. Temporal Features (Dynamic)

**Cumulative Features**:
```python
def create_cumulative_features(time_series_array):
    """Compute cumulative sums over time"""
    T, N = time_series_array.shape  # T timesteps, N nodes
    features = {}

    # Cumulative sum
    features['cumsum'] = np.cumsum(time_series_array, axis=0)

    # Cumulative max (peak so far)
    features['cummax'] = np.maximum.accumulate(time_series_array, axis=0)

    # Cumulative min
    features['cummin'] = np.minimum.accumulate(time_series_array, axis=0)

    return features
```

**Rolling Window Features**:
```python
def create_rolling_features(time_series_array, windows=[3, 6, 12]):
    """
    Create rolling window features.
    For 5-minute timesteps: windows=[3,6,12] = [15min, 30min, 1hour]
    """
    T, N = time_series_array.shape
    features = {}

    cumsum = np.cumsum(time_series_array, axis=0)

    for w in windows:
        # Current cumulative - w-step cumulative = sum of last w steps
        rolled = np.zeros_like(time_series_array)
        rolled[w:] = cumsum[:-w]
        features[f'roll_{w}'] = cumsum - rolled

    return features
```

**Rate of Change**:
```python
def create_rate_features(time_series_array):
    """Rate of change over time"""
    features = {}

    # First derivative
    features['delta_1'] = np.diff(time_series_array, axis=0, prepend=0)

    # Second derivative (acceleration)
    features['delta_2'] = np.diff(time_series_array, axis=0, n=2, prepend=0, append=0)

    # Percent change
    features['pct_change'] = np.diff(time_series_array, axis=0, prepend=time_series_array[:1])
    features['pct_change'] = features['pct_change'] / (time_series_array + 1e-6)

    return features
```

**Warm-up Features**:
```python
def create_warmup_features(time_series_array, warmup=10):
    """Statistics from initial warmup period"""
    T, N = time_series_array.shape
    features = {}

    # Compute stats from first 'warmup' timesteps
    for col in range(N):
        warmup_data = time_series_array[:warmup, col]

        # Broadcast same value to all timesteps
        features[f'warmup_mean_{col}'] = np.full(T, np.mean(warmup_data))
        features[f'warmup_std_{col}'] = np.full(T, np.std(warmup_data))
        features[f'warmup_min_{col}'] = np.full(T, np.min(warmup_data))
        features[f'warmup_max_{col}'] = np.full(T, np.max(warmup_data))

    return features
```

#### 3. Cross-Domain Features

**Multi-Modal Predictions**:
```python
# Train model on domain A, use predictions as features for domain B
model_a = train_model(X_a, y_a)

# Get predictions for domain B
preds_a_for_b = model_a.predict(X_b_features)

# Use as feature
X_b['pred_from_a'] = preds_a_for_b

# Example: 1D hydraulic network → 2D surface flooding
```

**Connection Mapping**:
```python
# Map nodes across domains
conn_map = defaultdict(list)
for _, row in connections.iterrows():
    node_a = row['node_idx_a']
    node_b = row['node_idx_b']
    conn_map[node_b].append(node_a)

# Aggregate predictions from connected nodes
def get_cross_domain_prediction(node_b, t):
    connected_nodes = conn_map.get(node_b, [])
    if not connected_nodes:
        return 0

    # Average predictions from connected nodes in domain A
    preds = [model_a.predict(get_features_a(n, t)) for n in connected_nodes]
    return np.mean(preds)
```

#### 4. Spatial-Temporal Interaction

**Spatiotemporal Lag**:
```python
def create_st_lag_features(node_data, adj, lags=[1, 2, 3]):
    """Features from neighbors at previous timesteps"""
    T, N = node_data.shape
    features = {}

    for lag in lags:
        neighbor_lag = np.zeros_like(node_data)
        for t in range(lag, T):
            for n in range(N):
                neighbors = adj.get(n, [])
                if neighbors:
                    # Average of neighbors at t-lag
                    neighbor_lag[t, n] = np.mean([
                        node_data[t-lag, neigh] for neigh in neighbors
                    ])
        features[f'neigh_lag_{lag}'] = neighbor_lag

    return features
```

**Diffusion Features**:
```python
def create_diffusion_features(node_data, adj, steps=2):
    """Simulate information diffusion across graph"""
    features = {}

    for step in range(1, steps + 1):
        diffused = node_data.copy()

        # Each step: spread to neighbors
        for _ in range(step):
            new_values = diffused.copy()
            for n in range(node_data.shape[1]):
                neighbors = adj.get(n, [])
                if neighbors:
                    # Average of self and neighbors
                    new_values[:, n] = np.mean(
                        np.column_stack([diffused[:, n]] + [diffused[:, neigh] for neigh in neighbors]),
                        axis=1
                    )
            diffused = new_values

        features[f'diffused_{step}'] = diffused

    return features
```

### Feature Selection Strategy

**1. Start Simple**:
```python
# Baseline
features = ['elevation', 'degree', 'rainfall']
```

**2. Add Neighbor Aggregations**:
```python
features += ['neigh_elev_mean', 'neigh_elev_std', 'neigh_count']
# Expected: +20-30% improvement
```

**3. Add Temporal Windows**:
```python
features += ['rain_cumsum', 'rain_roll_15', 'rain_roll_30', 'rain_roll_60']
# Expected: +10-15% improvement
```

**4. Add Cross-Domain**:
```python
features += ['pred_1d']
# Expected: +5-10% improvement
```

**5. Validate with Feature Importance**:
```python
importance = model.feature_importance()
feature_importance = sorted(zip(features, importance), key=lambda x: -x[1])

# Remove features with importance < 1% of max
max_imp = max(importance)
keep_features = [f for f, imp in feature_importance if imp > 0.01 * max_imp]
```

## Verification

**Expected Feature Importance Pattern**:
```
1. Neighbor aggregation (neigh_elev_mean)   ← Should be #1
2. Static node features (elevation, degree)
3. Temporal windows (rain_roll_15, rain_roll_30)
4. Cross-domain predictions (pred_1d)
5. Cumulative features (rain_cumsum)
```

**Red Flags**:
- ❌ Single-node features outrank neighbor features → Check graph construction
- ❌ Temporal features have near-zero importance → Wrong time windows
- ❌ Cross-domain feature has zero importance → Model A not useful for Model B

**Success Indicators**:
- ✅ Neighbor aggregation is top feature (importance > single-node)
- ✅ Adding features improves validation score monotonically
- ✅ Feature importance aligns with domain knowledge

## Example

**Task**: Urban flood water level prediction

**Data Structure**:
- 1D network: 17 nodes (pipes), features: inlet_flow
- 2D network: 3716 nodes (surface), features: rainfall
- Connections: 16 links between 1D and 2D

**Implementation**:
```python
# Step 1: Build graph features
def add_graph_features(nodes, edges):
    from_counts = pd.concat([edges['from_node'], edges['to_node']]).value_counts()
    nodes['degree'] = nodes['node_idx'].map(from_counts).fillna(0).astype(int)

    # Neighbor elevation stats
    elev_map = dict(zip(nodes['node_idx'], nodes['elevation']))
    adj = defaultdict(list)
    for _, row in edges.iterrows():
        adj[int(row['from_node'])].append(int(row['to_node']))

    neighbor_means = []
    for n in nodes['node_idx']:
        neighbors = adj.get(n, [])
        if neighbors:
            neighbor_elevs = [elev_map.get(x, elev_map[n]) for x in neighbors]
            neighbor_means.append(np.mean(neighbor_elevs))
        else:
            neighbor_means.append(elev_map[n])

    nodes['neigh_elev_mean'] = neighbor_means
    return nodes

nodes_2d = add_graph_features(nodes_2d, edges_2d)

# Step 2: Temporal features
rain_feats = create_rainfall_features(rainfall_2d.values)
# Returns: {'rcum': ..., 'r15': ..., 'r30': ..., 'r1h': ...}

# Step 3: Cross-domain (1D → 2D)
model_1d = train_1d_model(inlet_flow, water_level_1d)
for nid_2d in nodes_2d:
    connected_1d_nodes = get_connections(nid_2d)
    pred_1d = mean([model_1d.predict(n1d) for n1d in connected_1d_nodes])
    features_2d['pred_1d'] = pred_1d

# Step 4: Train final model
X = pd.concat([
    pd.DataFrame(rain_feats),
    nodes_2d[['elevation', 'degree', 'neigh_elev_mean']],
    pd.Series(pred_1d_for_2d.flatten())
], axis=1)

model = lgb.train(params, lgb.Dataset(X, y), num_boost_round=100)
```

**Results**:
```
Feature Importance:
  neigh_elev_mean    1239  ← TOP!
  elevation           1179
  degree               345
  rain_r15              77
  pred_1d               71
  rain_r30              52
  rain_rcum             25

RMSE: 1.0264 (vs ~2.5 without neighbor features)
Improvement: ~60% from neighbor aggregation alone
```

## Notes

**Computational Considerations**:
- Neighbor aggregation is O(E) where E = edges, efficient for sparse graphs
- For large graphs (>10K nodes), sample neighbors or use approximate aggregation
- Rolling windows can be pre-computed once per time series

**When to Use Each Feature Type**:
- **Neighbor stats**: Always useful in graph problems, highest ROI
- **Temporal windows**: When recent history matters more than distant
- **Cumulative**: When total accumulation matters (e.g., rainfall, traffic volume)
- **Cross-domain**: When multiple interacting systems (1D→2D, upstream→downstream)
- **Diffusion**: When propagation speed matters (epidemic, flood, cascade)

**Common Pitfalls**:
1. Using node IDs as features (they're not meaningful)
2. Ignoring graph structure and treating as tabular data
3. Not aggregating neighbor information (biggest missed opportunity)
4. Using wrong time windows (too short or too long)
5. Data leakage: including future information in features

**Domain-Specific Considerations**:
- **Hydraulic**: Elevation differences drive flow (use neighbor elevation)
- **Traffic**: Congestion propagates upstream (use reverse adjacency)
- **Power**: Flow follows path of least resistance (use impedance weights)
- **Social**: Influence decays with distance (use weighted aggregation)

## References

- Graph Neural Networks: `torch_geometric.nn.GCNConv`, `GATConv`
- Temporal Convolutional Networks for time series on graphs
- Urban Flood Modelling Kaggle Competition (0.3728 solution)