Exploring Visitor Census Data in U.S. Cities: Insights from Large-Scale Mobile Positioning Data

import sys, subprocess, importlib.util, os, warnings
from datetime import datetime
from typing import List, Optional

for pkg in ["geopandas", "matplotlib", "pytz"]:
    if importlib.util.find_spec(pkg) is None:
        subprocess.check_call([sys.executable, "-m", "pip", "install", pkg])

import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter
import seaborn as sns
import numpy as np
import pytz

warnings.filterwarnings("ignore")

def home_visit(timezone: str, start_date: int, end_date: int, demo_path, cbg_path,) -> pd.DataFrame:
    if start_date > end_date:
        raise ValueError("start_date must be ≤ end_date")

    # Load data
    demo_raw = pd.read_csv(demo_path, usecols=["device_key", "unix_time", "Target_CBG", "Visitor_CBG"])
    msa_cbg = gpd.read_file(cbg_path)[["GISJOIN"]]

    # Convert timestamps *vectorised and select time range
    true_tz = pytz.timezone(timezone)
    demo_raw['real_time'] = demo_raw['unix_time'].apply(lambda x: datetime.fromtimestamp(x,true_tz))
    demo_raw["hour"] = demo_raw["real_time"].dt.floor("H")
    demo_raw["day"]  = demo_raw["real_time"].dt.day
    demo_raw = demo_raw.loc[demo_raw["day"].between(start_date, end_date)]

    if demo_raw.empty:
        return pd.DataFrame(
            columns=["Time", "Target_CBG", "Visitor_CBG","Num_visitor", "Num_visit"])

    # Aggregate visits 
    hourly_counts = (demo_raw
        .groupby(["Target_CBG", "Visitor_CBG", "hour"], as_index=False)
        .agg(Num_visit=("device_key", "size"), Num_visitor=("device_key", "nunique")))

    # Attach missing CBGs (left join)
    result = (msa_cbg.merge(hourly_counts, left_on="GISJOIN", right_on="Target_CBG",how="left")
              .dropna(subset=["Num_visit", "Num_visitor"]))

    # ISO-style timestamp (yyyy-mm-dd HH)
    result["Time"] = result["hour"].dt.strftime("%Y-%m-%d %H")

    # reorder & clean
    result = result[["Time", "Target_CBG", "Visitor_CBG", "Num_visitor", "Num_visit"]].sort_values(["Time", "Target_CBG", "Visitor_CBG"]).reset_index(drop=True)
    return result

# TCMA (Central time), Juy 19-25 (all records in demo raw)
demo_path = "./visitorcensus_data/demo_data/demo_raw.csv"
cbg_path = "./visitorcensus_data/demo_data/TCMA_cbg.geojson"
demo_homevisit = home_visit('US/Central', 19, 25, demo_path, cbg_path)

demo_homevisit.info()

# optional: save the generated data
demo_homevisit.to_csv('./visitorcensus_data/demo_data/demo_homevisit.csv')

def visitor_census_aggregation(agg_type: str, homevisit_path, census_path, cbg_path,
    attrs: Optional[List[str]] = None,) -> gpd.GeoDataFrame:

    # Sanity checks
    if agg_type not in {"visitor", "visit"}:
        raise ValueError("agg_type must be 'visitor' or 'visit'")

    # Fall-back list of census columns
    if attrs is None:
        attrs = ["Income", "Wp", "Pbb", "Upp"]

    weight_col = "Num_visitor" if agg_type == "visitor" else "Num_visit"

    # Load data
    home_visit_tab = pd.read_csv(homevisit_path).drop(columns="Unnamed: 0", errors="ignore")
    attribute_tab  = pd.read_csv(census_path)   .drop(columns="Unnamed: 0", errors="ignore")

    # Link visitor attributes onto each visit record
    link = (home_visit_tab.merge(attribute_tab, left_on="Visitor_CBG", right_on="CBG_ID", how="left")
        .dropna(subset = attrs + [weight_col]))  # keep rows with data

    if link.empty:
        raise ValueError("No rows left after merging & dropping NaNs; check inputs.")

    # Group & compute weighted means
    def _wavg(group: pd.DataFrame) -> pd.Series:
        w = group[weight_col].to_numpy()
        out = {}
        for col in attrs:
            out[f"{col}_v"] = np.average(group[col].to_numpy(), weights=w)
        return pd.Series(out)

    visitor_attribute = (link.groupby("Target_CBG", sort=False)   # preserve input order
            .apply(_wavg).reset_index())

    # merge back onto full CBG
    msa_cbg = gpd.read_file(cbg_path)[["GISJOIN"]]
    visitor_census = (msa_cbg.merge(visitor_attribute, left_on="GISJOIN", right_on="Target_CBG", how="left")
               .drop(columns="Target_CBG")
               .rename(columns={"GISJOIN": "CBG_ID_v"}))
    visitor_census = visitor_census[["CBG_ID_v", "Income_v", "Wp_v", "Pbb_v", "Upp_v"]]
    return visitor_census

homevisit_path = "./visitorcensus_data/demo_data/demo_homevisit.csv"
census_path = "./visitorcensus_data/demo_data/census_attribute.csv"
cbg_path = "./visitorcensus_data/demo_data/TCMA_cbg.geojson"
visitor_census_1 = visitor_census_aggregation('visitor', homevisit_path, census_path, cbg_path)

visitor_census_1.info()

# optional
visitor_census_1.to_csv('./visitorcensus_data/demo_data/demo_visitor.csv')

homevisit_path = "./visitorcensus_data/demo_data/demo_homevisit.csv"
census_path = "./visitorcensus_data/demo_data/census_attribute.csv"
cbg_path = "./visitorcensus_data/demo_data/TCMA_cbg.geojson"
visitor_census_2 = visitor_census_aggregation('visit', homevisit_path, census_path, cbg_path)

visitor_census_2.info()

# optional
visitor_census_2.to_csv('./visitorcensus_data/demo_data/demo_visit.csv')

def plot_census_violin_grid(
    base_dir="./visitorcensus_data", census_type="visitor", month="July2021",
    msa_names=("LA", "Houston", "Atlanta", "Miami", "Seattle", "Denver", "TCMA"),
    metric_map=None, figsize=(11, 8), bg_color="white"
):
    if metric_map is None:
        metric_map = [
            ("Income ($)", "Income_v", "r_income", "#7ea6cf"),
            ("Wp (%)", "Wp_v", "r_wr", "#b7b3ea"),
            ("Pbb (%)", "Pbb_v", "r_pbb", "#83d0a5"),
            ("Upp (%)", "Upp_v", "r_upr", "#e7d3a3")
        ]

    data_dict = {}
    for msa in msa_names:
        other = pd.read_csv(os.path.join(base_dir, msa, f"{msa}_{census_type}_{month}.csv"))
        resident = pd.read_csv(os.path.join(base_dir, msa, f"{msa}_resident_{month}.csv"))
        data_dict[msa] = {"other": other, "resident": resident}

    xlims = {}
    for title, other_col, resident_col, color in metric_map:
        vals = []
        for msa in msa_names:
            vals.extend(data_dict[msa]["resident"][resident_col].dropna())
            vals.extend(data_dict[msa]["other"][other_col].dropna())
        vmin, vmax = np.nanmin(vals), np.nanmax(vals)
        pad = (vmax - vmin) * 0.08 if vmax > vmin else 1
        xlims[title] = (max(0, vmin - pad), vmax + pad)

    xlims["Upp (%)"] = (0, 0.2)

    fig, axes = plt.subplots(len(msa_names), len(metric_map), figsize=figsize, sharex="col")
    if len(msa_names) == 1:
        axes = np.array([axes])

    other_label = census_type.capitalize() + "s"

    for i, msa in enumerate(msa_names):
        resident, other = data_dict[msa]["resident"], data_dict[msa]["other"]

        for j, (title, other_col, resident_col, color) in enumerate(metric_map):
            ax = axes[i, j]
            x1, x2 = resident[resident_col].dropna(), other[other_col].dropna()

            data = pd.DataFrame({
                "Value": pd.concat([x1, x2], ignore_index=True),
                "Category": ["Residents"] * len(x1) + [other_label] * len(x2),
                "Group": [""] * (len(x1) + len(x2))
            })

            sns.violinplot(
                data=data, x="Value", y="Group", hue="Category",
                hue_order=["Residents", other_label], split=True,
                inner="quart", palette=[color, color], cut=0,
                linewidth=1.5, density_norm="width", ax=ax
            )

            if ax.get_legend() is not None:
                ax.get_legend().remove()

            m1, m2 = np.median(x1), np.median(x2)
            xmin, xmax = xlims[title]
            dx = (xmax - xmin) * 0.03
            fmt = "{:,.0f}" if "Income" in title else "{:.1f}"

            ax.axvline(m1, ymin=0.5, ymax=0.9, color="firebrick", lw=1.6)
            ax.axvline(m2, ymin=0.1, ymax=0.5, color="firebrick", lw=1.6)
            ax.text(m1 + dx, -0.12, fmt.format(m1), va="center", fontsize=10)
            ax.text(m2 + dx, 0.12, fmt.format(m2), va="center", fontsize=10)

            ax.set_xlim(xlims[title])
            ax.set_facecolor(bg_color)
            ax.set_xlabel("")
            ax.set_ylabel("")
            ax.set_yticks([])

            if j == 0:
                ax.text(-0.07, 0.72, "Residents", transform=ax.transAxes, ha="right", va="center", fontsize=12)
                ax.text(-0.07, 0.28, other_label, transform=ax.transAxes, ha="right", va="center", fontsize=12)

            for spine in ax.spines.values():
                spine.set_color("gray")
                spine.set_linewidth(0.8)

            if i == 0:
                ax.set_title(title, fontsize=14)

            if i < len(msa_names) - 1:
                ax.tick_params(labelbottom=False)

            if "Income" in title:
                ax.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: f"{int(x/1000)}k"))

        axes[i, 0].text(
            -0.6, 0.5, msa,
            transform=axes[i, 0].transAxes,
            ha="right", va="center",
            fontsize=12, fontweight="bold"
        )

    plt.tight_layout()
    return fig, axes

plot_census_violin_grid(census_type="visitor")

plot_census_violin_grid(census_type="visit")

demo_homevisit.head(5)

def plot_od_flows(
    flow_df,
    cbg_gdf,
    time_mode="all",
    hour=None,
    timestamp=None,
    start_time=None,
    end_time=None,
    weight_col="Num_visitor",
    min_flow=None,
    top_n=None,
    origin_col="Visitor_CBG",
    destination_col="Target_CBG",
    cbg_id_col="GISJOIN",
    figsize=(12, 12),
    alpha=0.35,
    linewidth_range=(0.2, 4)
):
    df = flow_df.copy()
    df["Time"] = pd.to_datetime(df["Time"])

    if weight_col not in ["Num_visitor", "Num_visit"]:
        raise ValueError("weight_col must be 'Num_visitor' or 'Num_visit'.")

    if start_time is not None:
        df = df[df["Time"] >= pd.to_datetime(start_time)]
    if end_time is not None:
        df = df[df["Time"] <= pd.to_datetime(end_time)]

    if time_mode == "hour":
        if hour is None:
            raise ValueError("Provide hour=0–23.")
        df = df[df["Time"].dt.hour == hour]
        time_label = f"{hour:02d}:00"
    elif time_mode == "timestamp":
        if timestamp is None:
            raise ValueError("Provide timestamp.")
        timestamp = pd.to_datetime(timestamp)
        df = df[df["Time"] == timestamp]
        time_label = timestamp.strftime("%Y-%m-%d %H:%M")
    elif time_mode == "all":
        time_label = "All selected times"
    else:
        raise ValueError("time_mode must be 'all', 'hour', or 'timestamp'.")

    od = (
        df.groupby([origin_col, destination_col], as_index=False)[weight_col]
        .sum()
        .rename(columns={weight_col: "flow"})
    )

    if min_flow is not None:
        od = od[od["flow"] >= min_flow]
    if top_n is not None:
        od = od.nlargest(top_n, "flow")

    cbg = cbg_gdf.copy()
    centroids = cbg.geometry.centroid
    coords = pd.DataFrame({
        cbg_id_col: cbg[cbg_id_col],
        "x": centroids.x,
        "y": centroids.y
    })

    od = od.merge(
        coords.rename(columns={cbg_id_col: origin_col, "x": "ox", "y": "oy"}),
        on=origin_col,
        how="left"
    ).merge(
        coords.rename(columns={cbg_id_col: destination_col, "x": "dx", "y": "dy"}),
        on=destination_col,
        how="left"
    ).dropna(subset=["ox", "oy", "dx", "dy"])

    fmin, fmax = od["flow"].min(), od["flow"].max()
    lw_min, lw_max = linewidth_range
    if fmax == fmin:
        od["linewidth"] = (lw_min + lw_max) / 2
    else:
        s = (od["flow"]**0.5 - fmin**0.5) / (fmax**0.5 - fmin**0.5)
        od["linewidth"] = lw_min + s * (lw_max - lw_min)

    fig, ax = plt.subplots(figsize=figsize)
    cbg.plot(ax=ax, facecolor="none", edgecolor="0.75", linewidth=0.25)

    for r in od.itertuples():
        ax.plot([r.ox, r.dx], [r.oy, r.dy], linewidth=r.linewidth, alpha=alpha)

    ax.set_title(f"TCMA OD Flows | {time_label} | {weight_col}")
    ax.set_axis_off()
    plt.tight_layout()

    return od, fig, ax

flows = demo_homevisit.copy()
cbg = gpd.read_file("./visitorcensus_data/demo_data/TCMA_cbg.geojson")

# Merge all times
od, fig, ax = plot_od_flows(
    flows, cbg,
    time_mode="all",
    weight_col="Num_visitor",
    min_flow=3
)

# Merge the same hour across all days
od, fig, ax = plot_od_flows(
    flows, cbg,
    time_mode="hour",
    hour=8,
    weight_col="Num_visit"
)