xjtulyc/mobility-trajectories

Urban Mobility Trajectory Analysis

One-line summary: Process GPS mobility traces to detect activity stops, build origin-destination matrices, measure individual mobility entropy, and model urban travel demand patterns.

---

When to Use This Skill

• When processing GPS/CDR traces to identify stops and trips
• When building origin-destination (OD) matrices from mobility data
• When measuring individual mobility entropy and radius of gyration
• When detecting home/work locations from trajectory data
• When analyzing commuting patterns and modal split
• When computing flow maps and desire lines

Trigger keywords: GPS trajectories, mobility data, origin-destination matrix, stop detection, radius of gyration, mobility entropy, CDR data, travel behavior, commuting, flow map, home detection, individual mobility, urban mobility, transportation

---

Background & Key Concepts

Radius of Gyration

$$ r_g = \sqrt{\frac{1}{N}\sum_{i=1}^N (r_i - r_{cm})^2} $$

where $r_{cm}$ is the center of mass of all visited locations. Measures spatial extent of individual mobility.

Mobility Entropy (Diversity of Locations)

$$ S = -\sum_{i} p_i \ln p_i $$

where $p_i$ = fraction of time spent at location $i$. High S = diverse, low S = routine/concentrated.

Stop Detection

A stop = sequence of GPS points where speed < threshold AND dwell time > minimum duration. Between stops = trips.

---

Environment Setup

Install Dependencies

pip install pandas>=2.0 numpy>=1.24 scipy>=1.11 scikit-learn>=1.3 \
            matplotlib>=3.7 geopandas>=0.14

Verify Installation

import pandas as pd
import numpy as np
from scipy.spatial.distance import cdist

# Haversine distance test
def haversine(lat1, lon1, lat2, lon2):
    R = 6371  # Earth radius km
    lat1, lon1, lat2, lon2 = map(np.radians, [lat1, lon1, lat2, lon2])
    dlat = lat2 - lat1; dlon = lon2 - lon1
    a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2
    return 2 * R * np.arcsin(np.sqrt(a))

dist = haversine(40.7128, -74.0060, 40.7589, -73.9851)  # NYC points
print(f"Distance: {dist:.3f} km  (expected ~5.7 km)")

---

Core Workflow

Step 1: GPS Trajectory Simulation and Stop Detection

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from datetime import datetime, timedelta

# ------------------------------------------------------------------ #
# Simulate GPS traces for 10 individuals over one week
# ------------------------------------------------------------------ #

def haversine(lat1, lon1, lat2, lon2):
    """Vectorized haversine distance in km."""
    R = 6371.0
    lat1, lon1, lat2, lon2 = (np.radians(x) for x in [lat1, lon1, lat2, lon2])
    dlat = lat2 - lat1; dlon = lon2 - lon1
    a = np.sin(dlat/2)**2 + np.cos(lat1)*np.cos(lat2)*np.sin(dlon/2)**2
    return 2 * R * np.arcsin(np.sqrt(np.clip(a, 0, 1)))

np.random.seed(42)
city_center = (40.7128, -74.0060)  # NYC

def generate_individual_trace(uid, n_days=7, city_lat=40.7128, city_lon=-74.0060):
    """
    Generate synthetic GPS trace for one individual.
    Activity pattern: home(night) → commute → work(day) → commute → home(evening)
    """
    records = []
    t = datetime(2023, 3, 6, 7, 0, 0)  # Start Monday

    # Fixed locations for this individual
    home_lat = city_lat + np.random.uniform(-0.05, 0.05)
    home_lon = city_lon + np.random.uniform(-0.05, 0.05)
    work_lat = city_lat + np.random.uniform(-0.03, 0.03)
    work_lon = city_lon + np.random.uniform(-0.03, 0.03)
    # Occasional visit locations
    pois = [(city_lat + np.random.uniform(-0.04, 0.04),
             city_lon + np.random.uniform(-0.04, 0.04)) for _ in range(5)]

    for day in range(n_days):
        # Morning: at home, GPS ping every 10min
        for h in range(0, 7):
            for m in range(0, 60, 10):
                records.append({'uid': uid, 'ts': t,
                                'lat': home_lat + np.random.normal(0, 0.0005),
                                'lon': home_lon + np.random.normal(0, 0.0005)})
                t += timedelta(minutes=10)

        # Commute (30 min)
        for step in range(6):
            frac = step / 5
            records.append({'uid': uid, 'ts': t,
                            'lat': home_lat + frac*(work_lat-home_lat) + np.random.normal(0, 0.001),
                            'lon': home_lon + frac*(work_lon-home_lon) + np.random.normal(0, 0.001)})
            t += timedelta(minutes=5)

        # At work (8 hours, weekdays only; walk around on weekends)
        if day < 5:
            for h in range(48):
                records.append({'uid': uid, 'ts': t,
                                'lat': work_lat + np.random.normal(0, 0.0005),
                                'lon': work_lon + np.random.normal(0, 0.0005)})
                t += timedelta(minutes=10)
        else:
            # Weekend: visit 1-2 POIs
            poi = pois[np.random.randint(len(pois))]
            for h in range(48):
                records.append({'uid': uid, 'ts': t,
                                'lat': poi[0] + np.random.normal(0, 0.001),
                                'lon': poi[1] + np.random.normal(0, 0.001)})
                t += timedelta(minutes=10)

        # Commute back
        for step in range(6):
            frac = step / 5
            records.append({'uid': uid, 'ts': t,
                            'lat': work_lat + frac*(home_lat-work_lat) + np.random.normal(0, 0.001),
                            'lon': work_lon + frac*(home_lon-work_lon) + np.random.normal(0, 0.001)})
            t += timedelta(minutes=5)

        # Evening at home
        for h in range(24):
            records.append({'uid': uid, 'ts': t,
                            'lat': home_lat + np.random.normal(0, 0.0005),
                            'lon': home_lon + np.random.normal(0, 0.0005)})
            t += timedelta(minutes=10)

    return records

# Generate traces for 10 individuals
all_records = []
for uid in range(10):
    all_records.extend(generate_individual_trace(uid))

df_gps = pd.DataFrame(all_records)
df_gps['ts'] = pd.to_datetime(df_gps['ts'])
df_gps = df_gps.sort_values(['uid', 'ts']).reset_index(drop=True)
print(f"GPS dataset: {len(df_gps):,} records, {df_gps['uid'].nunique()} individuals")

# ---- Stop Detection -------------------------------------------- #
def detect_stops(trace, speed_threshold_kmh=1.5, min_duration_min=15, radius_m=50):
    """
    Detect activity stops from GPS trace.

    Parameters
    ----------
    trace : DataFrame — sorted by 'ts', with 'lat', 'lon' columns
    speed_threshold_kmh : float — max speed to be considered stationary
    min_duration_min : float — minimum stop duration

    Returns
    -------
    DataFrame with stop records: lat, lon, arrive_ts, leave_ts, duration_min
    """
    trace = trace.copy().reset_index(drop=True)

    # Compute speed between consecutive points
    speeds = [0.0]
    for i in range(1, len(trace)):
        dt = (trace.loc[i,'ts'] - trace.loc[i-1,'ts']).total_seconds() / 3600  # hours
        d  = haversine(trace.loc[i-1,'lat'], trace.loc[i-1,'lon'],
                       trace.loc[i,'lat'],   trace.loc[i,'lon'])
        speeds.append(d/dt if dt > 0 else 0)
    trace['speed_kmh'] = speeds

    # Mark stationary points
    trace['stationary'] = trace['speed_kmh'] < speed_threshold_kmh

    # Cluster consecutive stationary points
    stops = []
    in_stop = False
    stop_pts = []

    for _, row in trace.iterrows():
        if row['stationary']:
            stop_pts.append(row)
            in_stop = True
        else:
            if in_stop and len(stop_pts) > 0:
                arrive = stop_pts[0]['ts']
                leave  = stop_pts[-1]['ts']
                dur_min = (leave - arrive).total_seconds() / 60
                if dur_min >= min_duration_min:
                    stops.append({
                        'lat': np.mean([p['lat'] for p in stop_pts]),
                        'lon': np.mean([p['lon'] for p in stop_pts]),
                        'arrive_ts': arrive,
                        'leave_ts': leave,
                        'duration_min': dur_min,
                        'n_points': len(stop_pts),
                    })
            stop_pts = []
            in_stop = False

    return pd.DataFrame(stops)

# Apply to all individuals
stops_all = []
for uid, group in df_gps.groupby('uid'):
    group_sorted = group.sort_values('ts')
    stops = detect_stops(group_sorted)
    if len(stops) > 0:
        stops['uid'] = uid
        stops_all.append(stops)

df_stops = pd.concat(stops_all, ignore_index=True)
print(f"\nDetected {len(df_stops)} stops across {df_stops['uid'].nunique()} individuals")
print(f"Mean stop duration: {df_stops['duration_min'].mean():.1f} min")
print(f"Median stops per person: {df_stops.groupby('uid').size().median():.0f}")

Step 2: Origin-Destination Matrix

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.cluster.hierarchy import linkage, dendrogram, fcluster

# ------------------------------------------------------------------ #
# Build OD matrix by clustering stop locations into zones
# ------------------------------------------------------------------ #

# Cluster stops into geographic zones (k-means on lat/lon)
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler

# Use all unique stop locations
stop_coords = df_stops[['lat', 'lon']].values
scaler = StandardScaler()
coords_scaled = scaler.fit_transform(stop_coords)

n_zones = 8  # Number of geographic zones
kmeans = KMeans(n_clusters=n_zones, random_state=42, n_init=10)
df_stops['zone'] = kmeans.fit_predict(coords_scaled)

zone_centers = scaler.inverse_transform(kmeans.cluster_centers_)
zone_centers_df = pd.DataFrame(zone_centers, columns=['lat', 'lon'])

# Assign zone names (simplified)
zone_names = [f"Zone {i+1}" for i in range(n_zones)]

# ---- Build OD matrix ------------------------------------------- #
def build_od_matrix(stops_df, n_zones):
    """
    Build OD matrix from consecutive stops (same individual).
    Each trip = from departure zone to arrival zone.
    """
    od = np.zeros((n_zones, n_zones), dtype=int)
    for uid, group in stops_df.groupby('uid'):
        group_sorted = group.sort_values('arrive_ts').reset_index(drop=True)
        for i in range(len(group_sorted) - 1):
            origin = int(group_sorted.loc[i, 'zone'])
            dest   = int(group_sorted.loc[i+1, 'zone'])
            if origin != dest:
                od[origin, dest] += 1
    return od

od_matrix = build_od_matrix(df_stops, n_zones)

# Print OD matrix
od_df = pd.DataFrame(od_matrix, index=zone_names, columns=zone_names)
print("\nOrigin-Destination Matrix (trips):")
print(od_df)

# ---- Visualize OD matrix and flow map -------------------------- #
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# OD matrix heatmap
im = axes[0].imshow(od_matrix, cmap='YlOrRd', aspect='auto')
plt.colorbar(im, ax=axes[0], label='Number of trips')
axes[0].set_xticks(range(n_zones)); axes[0].set_xticklabels(zone_names, rotation=30, ha='right', fontsize=8)
axes[0].set_yticks(range(n_zones)); axes[0].set_yticklabels(zone_names, fontsize=8)
axes[0].set_xlabel("Destination"); axes[0].set_ylabel("Origin")
axes[0].set_title("Origin-Destination Matrix")
for i in range(n_zones):
    for j in range(n_zones):
        if od_matrix[i,j] > 0:
            axes[0].text(j, i, str(od_matrix[i,j]), ha='center', va='center',
                          fontsize=8, color='white' if od_matrix[i,j] > od_matrix.max()/2 else 'black')

# Flow map (desire lines)
axes[1].scatter(zone_centers_df['lon'], zone_centers_df['lat'],
                c='red', s=200, zorder=5, edgecolors='black', linewidths=0.7)
for i, name in enumerate(zone_names):
    axes[1].annotate(name, (zone_centers_df.loc[i,'lon'], zone_centers_df.loc[i,'lat']),
                      fontsize=8, xytext=(5,5), textcoords='offset points')

max_flow = od_matrix.max()
for i in range(n_zones):
    for j in range(i+1, n_zones):
        flow = od_matrix[i,j] + od_matrix[j,i]
        if flow > 0:
            x = [zone_centers_df.loc[i,'lon'], zone_centers_df.loc[j,'lon']]
            y = [zone_centers_df.loc[i,'lat'],  zone_centers_df.loc[j,'lat']]
            lw = flow / max_flow * 5 + 0.5
            axes[1].plot(x, y, 'b-', linewidth=lw, alpha=0.5)

axes[1].set_xlabel("Longitude"); axes[1].set_ylabel("Latitude")
axes[1].set_title("Mobility Flow Map (Desire Lines)")
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("od_matrix.png", dpi=150)
plt.show()

Step 3: Mobility Metrics (Entropy, Radius of Gyration)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# ------------------------------------------------------------------ #
# Compute individual mobility metrics
# ------------------------------------------------------------------ #

def radius_of_gyration(stops_df, uid):
    """Radius of gyration for one individual (in km)."""
    user_stops = stops_df[stops_df['uid'] == uid][['lat','lon','duration_min']]
    if len(user_stops) == 0:
        return np.nan
    weights = user_stops['duration_min'].values
    weights = weights / weights.sum()
    lat_cm = np.average(user_stops['lat'], weights=weights)
    lon_cm = np.average(user_stops['lon'], weights=weights)
    dists = haversine(user_stops['lat'].values, user_stops['lon'].values,
                      np.full(len(user_stops), lat_cm),
                      np.full(len(user_stops), lon_cm))
    return np.sqrt(np.average(dists**2, weights=weights))

def mobility_entropy(stops_df, uid, zone_col='zone'):
    """Shannon entropy over zone visit time fractions."""
    user_stops = stops_df[stops_df['uid'] == uid]
    if len(user_stops) == 0:
        return np.nan
    zone_time = user_stops.groupby(zone_col)['duration_min'].sum()
    p = zone_time / zone_time.sum()
    p = p[p > 0]
    return -np.sum(p * np.log(p))

def detect_home(stops_df, uid, night_hours=(22, 6)):
    """Detect home location: most frequent stop during night hours."""
    user_stops = stops_df[stops_df['uid'] == uid].copy()
    hour = user_stops['arrive_ts'].dt.hour
    night_mask = (hour >= night_hours[0]) | (hour < night_hours[1])
    night_stops = user_stops[night_mask]
    if len(night_stops) == 0:
        return np.nan, np.nan
    home_zone = night_stops.groupby('zone')['duration_min'].sum().idxmax()
    home_loc = night_stops[night_stops['zone'] == home_zone][['lat','lon']].mean()
    return home_loc['lat'], home_loc['lon']

# Compute metrics for all individuals
mobility_metrics = []
for uid in df_stops['uid'].unique():
    rg  = radius_of_gyration(df_stops, uid)
    S   = mobility_entropy(df_stops, uid)
    h_lat, h_lon = detect_home(df_stops, uid)
    n_stops_per_day = df_stops[df_stops['uid']==uid]['arrive_ts'].apply(lambda x: x.date()).nunique()
    mobility_metrics.append({
        'uid': uid,
        'n_stops': len(df_stops[df_stops['uid']==uid]),
        'radius_of_gyration_km': rg,
        'mobility_entropy': S,
        'n_zones_visited': df_stops[df_stops['uid']==uid]['zone'].nunique(),
        'home_lat': h_lat, 'home_lon': h_lon,
    })

df_mobility = pd.DataFrame(mobility_metrics)
print("Individual Mobility Metrics:")
print(df_mobility[['uid','n_stops','radius_of_gyration_km','mobility_entropy','n_zones_visited']].round(3).to_string(index=False))

# ---- Visualization --------------------------------------------- #
fig, axes = plt.subplots(2, 2, figsize=(12, 9))

# Radius of gyration distribution
axes[0][0].hist(df_mobility['radius_of_gyration_km'].dropna(), bins=8,
                color='steelblue', edgecolor='black', linewidth=0.7)
axes[0][0].set_xlabel("Radius of Gyration (km)"); axes[0][0].set_ylabel("Count")
axes[0][0].set_title("Radius of Gyration Distribution"); axes[0][0].grid(axis='y', alpha=0.3)

# Mobility entropy distribution
axes[0][1].hist(df_mobility['mobility_entropy'].dropna(), bins=8,
                color='coral', edgecolor='black', linewidth=0.7)
axes[0][1].set_xlabel("Mobility Entropy (nats)"); axes[0][1].set_ylabel("Count")
axes[0][1].set_title("Mobility Entropy Distribution"); axes[0][1].grid(axis='y', alpha=0.3)

# Rg vs. entropy scatter
axes[1][0].scatter(df_mobility['radius_of_gyration_km'],
                    df_mobility['mobility_entropy'],
                    c='purple', s=80, edgecolors='black', linewidths=0.5)
for _, row in df_mobility.iterrows():
    axes[1][0].annotate(f"U{row['uid']}", (row['radius_of_gyration_km'], row['mobility_entropy']),
                         fontsize=8, xytext=(3,3), textcoords='offset points')
axes[1][0].set_xlabel("Radius of Gyration (km)"); axes[1][0].set_ylabel("Mobility Entropy")
axes[1][0].set_title("Mobility Extent vs. Diversity"); axes[1][0].grid(True, alpha=0.3)

# Stops per day by zone (heatmap)
zone_visit_matrix = df_stops.groupby(['uid','zone'])['duration_min'].sum().unstack(fill_value=0)
im = axes[1][1].imshow(zone_visit_matrix.values, cmap='YlOrRd', aspect='auto')
plt.colorbar(im, ax=axes[1][1], label='Total time (min)')
axes[1][1].set_xticks(range(n_zones)); axes[1][1].set_xticklabels(zone_names, rotation=30, ha='right', fontsize=8)
axes[1][1].set_yticks(range(len(zone_visit_matrix))); axes[1][1].set_yticklabels([f"U{uid}" for uid in zone_visit_matrix.index], fontsize=8)
axes[1][1].set_xlabel("Zone"); axes[1][1].set_ylabel("Individual")
axes[1][1].set_title("Time Spent per Zone per Individual")

plt.tight_layout()
plt.savefig("mobility_metrics.png", dpi=150)
plt.show()

---

Advanced Usage

Commuting Pattern Analysis

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# ---- Hourly trip departure distribution ------------------------ #
df_stops['hour_of_day'] = df_stops['arrive_ts'].dt.hour
trip_hours = []
for uid, group in df_stops.groupby('uid'):
    g = group.sort_values('arrive_ts')
    for i in range(len(g)-1):
        o_zone = g.iloc[i]['zone']
        d_zone = g.iloc[i+1]['zone']
        if o_zone != d_zone:
            trip_hours.append(g.iloc[i]['leave_ts'].hour if 'leave_ts' in g.columns else g.iloc[i]['arrive_ts'].hour)

fig, ax = plt.subplots(figsize=(10, 4))
ax.hist(trip_hours, bins=24, range=(0,24), color='steelblue', edgecolor='black', linewidth=0.5)
ax.set_xlabel("Hour of day"); ax.set_ylabel("Number of trips")
ax.set_title("Trip Departure Distribution by Hour"); ax.grid(axis='y', alpha=0.3)
ax.set_xticks(range(0, 25, 2)); ax.axvspan(7, 9, alpha=0.1, color='red', label='AM peak')
ax.axvspan(17, 19, alpha=0.1, color='orange', label='PM peak')
ax.legend(); plt.tight_layout(); plt.savefig("trip_hours.png", dpi=150); plt.show()

---

Troubleshooting

Stop detection misses short stops

# Reduce minimum duration
stops = detect_stops(trace, min_duration_min=5)  # 5 min instead of 15

GPS traces with gaps (network outage)

# Filter out inter-point intervals > threshold
df_gps['time_gap'] = df_gps.groupby('uid')['ts'].diff().dt.total_seconds() / 60
df_gps = df_gps[df_gps['time_gap'] < 60]  # Remove gaps > 60 min

Haversine distance overflow for antipodal points

# Clip arcsin argument to [-1, 1]
return 2 * R * np.arcsin(np.sqrt(np.clip(a, 0, 1)))

Version Compatibility

| Package | Tested versions | Notes | |:--------|:----------------|:------| | pandas | 2.0, 2.1 | dt.hour, groupby stable | | scikit-learn | 1.3, 1.4 | KMeans API stable | | geopandas | 0.14 | For spatial joining with zones |

---

External Resources

Official Documentation

• pandas datetime

Key Papers

• González, M.C., Hidalgo, C.A. & Barabási, A.-L. (2008). Understanding individual human mobility patterns. Nature.
• Song, C. et al. (2010). Limits of predictability in human mobility. Science.

---

Examples

Example 1: Daily Mobility Profile

import numpy as np
import matplotlib.pyplot as plt

# Average time at home vs. outside by hour
for uid in range(3):
    user_trace = df_gps[df_gps['uid'] == uid].copy()
    user_trace['hour'] = user_trace['ts'].dt.hour
    # Simplified: count records near home
    home_lat = user_trace.groupby('uid').apply(lambda x: x.iloc[0]['lat']).values[0] if len(user_trace) > 0 else 40.71
    home_lon = user_trace.groupby('uid').apply(lambda x: x.iloc[0]['lon']).values[0] if len(user_trace) > 0 else -74.01

print("Daily profile analysis requires full trace data — see Step 1 output.")

Example 2: Mobility Predictability (Maximum Entropy)

import numpy as np

# Maximum entropy = log(n_zones) if uniformly distributed
n_zones_visited = df_mobility['n_zones_visited'].values
max_entropy = np.log(n_zones_visited)
actual_entropy = df_mobility['mobility_entropy'].values
predictability = 1 - actual_entropy / (max_entropy + 1e-8)

print("Individual mobility predictability (1 = perfectly predictable):")
for uid, pred in zip(df_mobility['uid'], predictability):
    print(f"  User {uid}: {pred:.3f}")
print(f"\nMean predictability: {predictability.mean():.3f}")

---

Last updated: 2026-03-17 | Maintainer: @xjtulyc Issues: GitHub Issues