{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.7.10","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nfrom pathlib import Path\n \nTRAIN_DIR  = Path(\"/kaggle/input/google-smartphone-decimeter-challenge/train\")\nBASELINE   = Path(\"/kaggle/input/google-smartphone-decimeter-challenge/baseline_locations_train.csv\")\nOUTPUT     = Path(\"/kaggle/working/geofusion_full_feature_dataset.csv\")\n \nGPS_EPOCH_OFFSET_MS = 315_964_800_000   # milliseconds between Unix epoch and GPS epoch","metadata":{"execution":{"iopub.status.busy":"2026-07-07T21:10:39.019978Z","iopub.execute_input":"2026-07-07T21:10:39.020541Z","iopub.status.idle":"2026-07-07T21:10:39.028446Z","shell.execute_reply.started":"2026-07-07T21:10:39.020491Z","shell.execute_reply":"2026-07-07T21:10:39.027251Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def llh_to_ecef(lat_deg, lng_deg, h_m):\n    \"\"\"Convert WGS-84 geodetic coordinates to ECEF (x, y, z) in metres.\"\"\"\n    a  = 6_378_137.0            # semi-major axis (m)\n    e2 = 6.6943799901414e-3     # first eccentricity squared\n    lat = np.radians(lat_deg)\n    lng = np.radians(lng_deg)\n    N   = a / np.sqrt(1 - e2 * np.sin(lat) ** 2)\n    x   = (N + h_m) * np.cos(lat) * np.cos(lng)\n    y   = (N + h_m) * np.cos(lat) * np.sin(lng)\n    z   = (N * (1 - e2) + h_m)   * np.sin(lat)\n    return x, y, z","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-07T21:10:39.030039Z","iopub.execute_input":"2026-07-07T21:10:39.030365Z","iopub.status.idle":"2026-07-07T21:10:39.041106Z","shell.execute_reply.started":"2026-07-07T21:10:39.030324Z","shell.execute_reply":"2026-07-07T21:10:39.039948Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def ecef_diff_to_north_east(lat_ref_deg, lng_ref_deg,\n                             dx, dy, dz):\n    \"\"\"\n    Rotate an ECEF difference vector (phone − NovAtel) into the local\n    North-East frame at the reference (NovAtel) position.\n    Returns (error_N_m, error_E_m).\n    \"\"\"\n    lat = np.radians(lat_ref_deg)\n    lng = np.radians(lng_ref_deg)\n    # ENU rotation matrix rows for North and East only\n    # North = [-sin(lat)cos(lng), -sin(lat)sin(lng),  cos(lat)]\n    # East  = [-sin(lng),          cos(lng),           0       ]\n    north = (-np.sin(lat) * np.cos(lng) * dx\n             - np.sin(lat) * np.sin(lng) * dy\n             + np.cos(lat) * dz)\n    east  = (-np.sin(lng) * dx\n             + np.cos(lng) * dy)\n    return north, east","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-07T21:10:39.043592Z","iopub.execute_input":"2026-07-07T21:10:39.043937Z","iopub.status.idle":"2026-07-07T21:10:39.061757Z","shell.execute_reply.started":"2026-07-07T21:10:39.043894Z","shell.execute_reply":"2026-07-07T21:10:39.060606Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def parse_gnss_status(gnss_log_path):\n    \"\"\"\n    Extract Status records from a GnssLog file.\n    Returns a DataFrame with one row per signal per epoch,\n    with millisSinceGpsEpoch as the epoch key.\n    \"\"\"\n    rows = []\n    with open(gnss_log_path, \"r\", errors=\"replace\") as f:\n        for line in f:\n            if not line.startswith(\"Status,\"):\n                continue\n            p = line.strip().split(\",\")\n            try:\n                unix_ms = int(p[1])\n                if unix_ms <= 0:\n                    continue\n                rows.append({\n                    \"millisSinceGpsEpoch\": unix_ms - GPS_EPOCH_OFFSET_MS,\n                    \"cn0DbHz\":             float(p[7]),\n                    \"azimuthDeg\":          float(p[8]),\n                    \"elevationDeg\":        float(p[9]),\n                    \"usedInFix\":           int(p[10]),\n                })\n            except (IndexError, ValueError):\n                continue\n    return pd.DataFrame(rows)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-07T21:10:39.063809Z","iopub.execute_input":"2026-07-07T21:10:39.064227Z","iopub.status.idle":"2026-07-07T21:10:39.076491Z","shell.execute_reply.started":"2026-07-07T21:10:39.064166Z","shell.execute_reply":"2026-07-07T21:10:39.075414Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"STATUS_EMPTY_COLS = [\n    \"millisSinceGpsEpoch\", \"n_sats_used\",\n    \"avg_cn0\", \"min_cn0\", \"max_cn0\",\n    \"avg_elevation\", \"min_elevation\", \"max_elevation\",\n    \"avg_azimuth\",\n]\n \ndef aggregate_status(status_df):\n    \"\"\"\n    Aggregate per-signal Status records to one row per epoch.\n    Returns a DataFrame with STATUS_EMPTY_COLS even when input is empty.\n    \"\"\"\n    if status_df.empty:\n        return pd.DataFrame(columns=STATUS_EMPTY_COLS)\n    return (\n        status_df\n        .groupby(\"millisSinceGpsEpoch\", sort=False)\n        .agg(\n            n_sats_used   = (\"usedInFix\",    \"sum\"),\n            avg_cn0       = (\"cn0DbHz\",      \"mean\"),\n            min_cn0       = (\"cn0DbHz\",      \"min\"),\n            max_cn0       = (\"cn0DbHz\",      \"max\"),\n            avg_elevation = (\"elevationDeg\", \"mean\"),\n            min_elevation = (\"elevationDeg\", \"min\"),\n            max_elevation = (\"elevationDeg\", \"max\"),\n            avg_azimuth   = (\"azimuthDeg\",   \"mean\"),\n        )\n        .reset_index()\n    )","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-07T21:10:39.077819Z","iopub.execute_input":"2026-07-07T21:10:39.078156Z","iopub.status.idle":"2026-07-07T21:10:39.093799Z","shell.execute_reply.started":"2026-07-07T21:10:39.078123Z","shell.execute_reply":"2026-07-07T21:10:39.092785Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"DERIVED_EMPTY_COLS = [\n    \"millisSinceGpsEpoch\",\n    \"avg_iono\", \"avg_tropo\", \"avg_rawPrUnc\", \"n_signals\",\n]\n \ndef aggregate_derived(derived_df):\n    \"\"\"\n    Aggregate per-signal derived measurements to one row per epoch.\n    Returns a DataFrame with DERIVED_EMPTY_COLS even when input is empty.\n    \"\"\"\n    if derived_df.empty:\n        return pd.DataFrame(columns=DERIVED_EMPTY_COLS)\n    return (\n        derived_df\n        .groupby(\"millisSinceGpsEpoch\", sort=False)\n        .agg(\n            avg_iono     = (\"ionoDelayM\",  \"mean\"),\n            avg_tropo    = (\"tropoDelayM\", \"mean\"),\n            avg_rawPrUnc = (\"rawPrUncM\",   \"mean\"),\n            n_signals    = (\"rawPrUncM\",   \"count\"),\n        )\n        .reset_index()\n    )","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-07T21:10:39.095811Z","iopub.execute_input":"2026-07-07T21:10:39.096312Z","iopub.status.idle":"2026-07-07T21:10:39.106584Z","shell.execute_reply.started":"2026-07-07T21:10:39.096274Z","shell.execute_reply":"2026-07-07T21:10:39.105517Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def compute_j_avg(derived_df, baseline_epoch_df):\n    \"\"\"\n    Compute J_avg (El Abbous satellite geometry cost) per epoch.\n \n    derived_df       : derived.csv for one phone, with satellite ECEF columns\n    baseline_epoch_df: baseline_locations_train rows for this phone,\n                       indexed by millisSinceGpsEpoch, with columns\n                       latDeg, lngDeg, heightAboveWgs84EllipsoidM\n \n    Returns a DataFrame with columns [millisSinceGpsEpoch, j_avg].\n    \"\"\"\n    # Merge receiver position (from baseline) into derived signals\n    der = derived_df.merge(\n        baseline_epoch_df[[\"millisSinceGpsEpoch\",\n                            \"latDeg\", \"lngDeg\",\n                            \"heightAboveWgs84EllipsoidM\"]],\n        on=\"millisSinceGpsEpoch\", how=\"inner\"\n    )\n    if der.empty:\n        return pd.DataFrame(columns=[\"millisSinceGpsEpoch\", \"j_avg\"]).astype(\n            {\"millisSinceGpsEpoch\": \"int64\", \"j_avg\": \"float64\"}\n        )\n \n    # Receiver ECEF from phone's reported position (no leakage)\n    rx_x, rx_y, rx_z = llh_to_ecef(\n        der[\"latDeg\"].values,\n        der[\"lngDeg\"].values,\n        der[\"heightAboveWgs84EllipsoidM\"].values\n    )\n \n    # LOS vector: satellite ECEF − receiver ECEF\n    dx = der[\"xSatPosM\"].values - rx_x\n    dy = der[\"ySatPosM\"].values - rx_y\n    dz = der[\"zSatPosM\"].values - rx_z\n    norm = np.sqrt(dx**2 + dy**2 + dz**2)\n    # Avoid division by zero\n    valid = norm > 0\n    dx[valid] /= norm[valid]\n    dy[valid] /= norm[valid]\n    dz[valid] /= norm[valid]\n \n    der[\"los_x\"] = dx\n    der[\"los_y\"] = dy\n    der[\"los_z\"] = dz\n \n    results = []\n    for epoch, grp in der.groupby(\"millisSinceGpsEpoch\", sort=False):\n        los = grp[[\"los_x\", \"los_y\", \"los_z\"]].values\n        n   = len(los)\n        if n < 2:\n            results.append({\"millisSinceGpsEpoch\": epoch, \"j_avg\": np.nan})\n            continue\n        # Pairwise cosines of angles between LOS vectors\n        cos_theta = np.clip(los @ los.T, -1.0, 1.0)        # (n, n)\n        J_matrix  = np.cos(2 * np.arccos(cos_theta))        # J_ij = cos(2θ)\n        np.fill_diagonal(J_matrix, 0)                        # exclude self-pairs\n        J_i = J_matrix.sum(axis=1)                           # per-satellite cost\n        results.append({\n            \"millisSinceGpsEpoch\": epoch,\n            \"j_avg\": J_i.mean()\n        })\n \n    return pd.DataFrame(results)\n ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-07T21:10:39.108116Z","iopub.execute_input":"2026-07-07T21:10:39.108474Z","iopub.status.idle":"2026-07-07T21:10:39.125926Z","shell.execute_reply.started":"2026-07-07T21:10:39.108441Z","shell.execute_reply":"2026-07-07T21:10:39.124714Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"Loading baseline locations...\")\nbaseline_all = pd.read_csv(BASELINE)\nprint(f\"  Baseline shape: {baseline_all.shape}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-07T21:10:39.128293Z","iopub.execute_input":"2026-07-07T21:10:39.128643Z","iopub.status.idle":"2026-07-07T21:10:39.332305Z","shell.execute_reply.started":"2026-07-07T21:10:39.128612Z","shell.execute_reply":"2026-07-07T21:10:39.331135Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"all_epochs = []\n \ndrive_dirs = sorted([d for d in TRAIN_DIR.iterdir() if d.is_dir()])\nprint(f\"Found {len(drive_dirs)} drives.\\n\")\n \nfor drive_dir in drive_dirs:\n    drive_id = drive_dir.name\n    phone_dirs = sorted([p for p in drive_dir.iterdir() if p.is_dir()])\n \n    for phone_dir in phone_dirs:\n        phone_name = phone_dir.name\n        print(f\"Processing: {drive_id} / {phone_name}\")\n \n        # ── Locate required files ──────────────────────────────────────────────\n        gt_path      = phone_dir / \"ground_truth.csv\"\n        derived_path = phone_dir / f\"{phone_name}_derived.csv\"\n        gnss_path    = phone_dir / f\"{phone_name}_GnssLog.txt\"\n \n        if not gt_path.exists():\n            print(f\"  [SKIP] ground_truth.csv not found\")\n            continue\n        if not derived_path.exists():\n            print(f\"  [SKIP] derived.csv not found\")\n            continue\n        if not gnss_path.exists():\n            print(f\"  [SKIP] GnssLog not found\")\n            continue\n \n        # ── Step 1: Load ground truth ──────────────────────────────────────────\n        gt = pd.read_csv(gt_path)\n \n        # ── Step 1: Compute error labels ───────────────────────────────────────\n        # Phone's reported position for this drive/phone\n        baseline = baseline_all[\n            (baseline_all[\"collectionName\"] == drive_id) &\n            (baseline_all[\"phoneName\"]      == phone_name)\n        ][[\"millisSinceGpsEpoch\", \"latDeg\", \"lngDeg\",\n           \"heightAboveWgs84EllipsoidM\"]].copy()\n \n        # Merge ground truth with baseline on epoch key\n        merged = gt.merge(baseline, on=\"millisSinceGpsEpoch\", how=\"inner\",\n                          suffixes=(\"_gt\", \"_phone\"))\n \n        if merged.empty:\n            print(f\"  [SKIP] No matching epochs between ground truth and baseline\")\n            continue\n \n        # Convert both positions to ECEF\n        gt_x, gt_y, gt_z = llh_to_ecef(\n            merged[\"latDeg_gt\"].values,\n            merged[\"lngDeg_gt\"].values,\n            merged[\"heightAboveWgs84EllipsoidM_gt\"].values\n        )\n        ph_x, ph_y, ph_z = llh_to_ecef(\n            merged[\"latDeg_phone\"].values,\n            merged[\"lngDeg_phone\"].values,\n            merged[\"heightAboveWgs84EllipsoidM_phone\"].values\n        )\n \n        # ECEF difference: phone − NovAtel\n        dx = ph_x - gt_x\n        dy = ph_y - gt_y\n        dz = ph_z - gt_z\n \n        # Rotate into local North-East frame\n        error_N, error_E = ecef_diff_to_north_east(\n            merged[\"latDeg_gt\"].values,\n            merged[\"lngDeg_gt\"].values,\n            dx, dy, dz\n        )\n        merged[\"error_N_m\"]     = error_N\n        merged[\"error_E_m\"]     = error_E\n        merged[\"radial_error_m\"] = np.sqrt(error_N**2 + error_E**2)\n \n        # ── Step 2: Vehicle state and DOP already in ground truth ──────────────\n        # (speedMps, courseDegree, hDop, vDop, timeSinceFirstFixSeconds)\n        # carried from gt via merge above\n \n        # ── Step 3: Aggregate GnssLog Status records ───────────────────────────\n        status_raw  = parse_gnss_status(gnss_path)\n        status_agg  = aggregate_status(status_raw)\n \n        # ── Step 4: Aggregate derived.csv ──────────────────────────────────────\n        derived     = pd.read_csv(derived_path)\n        derived_agg = aggregate_derived(derived)\n \n        # ── Step 5: Compute J_avg ──────────────────────────────────────────────\n        j_avg_df = compute_j_avg(derived, baseline)\n \n        # ── Step 6: Join all epoch-level features ──────────────────────────────\n        # Ensure all right-hand DataFrames have the merge key even if empty\n        if \"millisSinceGpsEpoch\" not in status_agg.columns:\n            status_agg  = pd.DataFrame(columns=STATUS_EMPTY_COLS)\n        if \"millisSinceGpsEpoch\" not in derived_agg.columns:\n            derived_agg = pd.DataFrame(columns=DERIVED_EMPTY_COLS)\n        if \"millisSinceGpsEpoch\" not in j_avg_df.columns:\n            j_avg_df    = pd.DataFrame(columns=[\"millisSinceGpsEpoch\", \"j_avg\"])\n \n        epoch_df = (merged[[\n                \"collectionName\", \"phoneName\", \"millisSinceGpsEpoch\",\n                \"latDeg_gt\", \"lngDeg_gt\",\n                \"latDeg_phone\", \"lngDeg_phone\",\n                \"radial_error_m\", \"error_N_m\", \"error_E_m\",\n                \"speedMps\", \"courseDegree\",\n                \"hDop\", \"vDop\", \"timeSinceFirstFixSeconds\"\n            ]]\n            .merge(status_agg,  on=\"millisSinceGpsEpoch\", how=\"left\")\n            .merge(derived_agg, on=\"millisSinceGpsEpoch\", how=\"left\")\n            .merge(j_avg_df,    on=\"millisSinceGpsEpoch\", how=\"left\")\n        )\n \n        # Convert millisSinceGpsEpoch to Unix seconds for readability\n        epoch_df[\"epoch_unix_s\"] = (\n            epoch_df[\"millisSinceGpsEpoch\"] // 1000 + GPS_EPOCH_OFFSET_MS // 1000\n        )\n        epoch_df.drop(columns=[\"millisSinceGpsEpoch\"], inplace=True)\n \n        print(f\"  → {len(epoch_df):,} epochs, {epoch_df.shape[1]} columns\")\n        all_epochs.append(epoch_df)\n ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-07T21:10:39.333843Z","iopub.execute_input":"2026-07-07T21:10:39.334146Z","iopub.status.idle":"2026-07-07T21:14:08.790863Z","shell.execute_reply.started":"2026-07-07T21:10:39.334108Z","shell.execute_reply":"2026-07-07T21:14:08.789458Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"\\nConcatenating all sessions...\")\nfinal_df = pd.concat(all_epochs, ignore_index=True)\nprint(f\"Raw dataset shape: {final_df.shape}\")\n\n# Reorder columns cleanly\ncol_order = [\n    \"collectionName\", \"phoneName\", \"epoch_unix_s\",\n    \"latDeg_gt\", \"lngDeg_gt\",\n    \"latDeg_phone\", \"lngDeg_phone\",\n    \"radial_error_m\", \"error_N_m\", \"error_E_m\",\n    \"speedMps\", \"courseDegree\",\n    \"hDop\", \"vDop\", \"timeSinceFirstFixSeconds\",\n    \"n_sats_used\",\n    \"avg_cn0\", \"min_cn0\", \"max_cn0\",\n    \"avg_elevation\", \"min_elevation\", \"max_elevation\",\n    \"avg_azimuth\",\n    \"avg_iono\", \"avg_tropo\",\n    \"avg_rawPrUnc\",\n    \"n_signals\",\n    \"j_avg\",\n]\nfinal_df = final_df[[c for c in col_order if c in final_df.columns]]\n\n# ── Cleaning Step 1: Drop columns with >50% nulls ──────────────────────────\n# Exception: retain the 42.1% group (n_signals, avg_iono, avg_rawPrUnc,\n# avg_tropo, j_avg) — these will be handled by row-level dropping below.\nnull_pct = final_df.isna().mean()\ndrop_cols = [c for c in null_pct[null_pct > 0.5].index]\nfinal_df = final_df.drop(columns=drop_cols)\nprint(f\"After dropping {len(drop_cols)} high-null columns (>40%): {final_df.shape}\")\nprint(f\"  Dropped: {drop_cols}\")\n\n# ── Cleaning Step 2: Drop rows with any remaining nulls ────────────────────\nbefore = len(final_df)\nfinal_df = final_df.dropna()\nprint(f\"After dropping null rows: {len(final_df):,}  (removed {before - len(final_df):,} rows)\")\n\n# ── Cleaning Step 3: Remove outliers beyond ±5 IQR in N and E error ───────\nfor col in ['error_N_m', 'error_E_m']:\n    q1, q3 = final_df[col].quantile(0.25), final_df[col].quantile(0.75)\n    iqr = q3 - q1\n    lo, hi = q1 - 5 * iqr, q3 + 5 * iqr\n    before = len(final_df)\n    final_df = final_df[(final_df[col] >= lo) & (final_df[col] <= hi)]\n    print(f\"Outlier removal {col}: kept {len(final_df):,}  \"\n          f\"(removed {before - len(final_df):,})  range=[{lo:.2f}, {hi:.2f}] m\")\n\nprint(f\"\\nFinal clean dataset shape: {final_df.shape}\")\nassert final_df.isna().sum().sum() == 0, \"Nulls remain after cleaning — check pipeline.\"\n\nprint(f\"\\nColumn summary:\")\nprint(final_df.dtypes)\nprint(f\"\\nSample statistics:\")\nprint(final_df.describe())\n\nfinal_df.to_csv(OUTPUT, index=False)\nprint(f\"\\nSaved to {OUTPUT}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-07T21:14:08.793161Z","iopub.execute_input":"2026-07-07T21:14:08.793868Z","iopub.status.idle":"2026-07-07T21:14:10.608900Z","shell.execute_reply.started":"2026-07-07T21:14:08.793818Z","shell.execute_reply":"2026-07-07T21:14:10.607507Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}