Python 64.6%
TeX 23.3%
HTML 11.6%
Makefile 0.5%
1"""Geospatial feature engineering for hedonic pricing analysis."""23from __future__ import annotations45import logging6from typing import Any78import numpy as np9import pandas as pd10from sklearn.cluster import KMeans11from sklearn.neighbors import BallTree1213logger = logging.getLogger(__name__)1415EARTH_RADIUS_M = 6_371_000.01617# ---------------------------------------------------------------------------18# Default city centres (major Quebec cities)19# ---------------------------------------------------------------------------2021DEFAULT_CITY_CENTERS: dict[str, tuple[float, float]] = {22 "montréal": (45.5017, -73.5673),23 "montreal": (45.5017, -73.5673),24 "québec": (46.8139, -71.2080),25 "quebec": (46.8139, -71.2080),26 "gatineau": (45.4765, -75.7013),27 "sherbrooke": (45.4042, -71.8929),28 "saguenay": (48.4280, -71.0686),29 "laval": (45.6066, -73.7124),30 "trois-rivières": (46.3432, -72.5432),31 "trois-rivieres": (46.3432, -72.5432),32 "lévis": (46.8032, -71.1780),33 "levis": (46.8032, -71.1780),34 "terrebonne": (45.6960, -73.6473),35 "saint-jean-sur-richelieu": (45.3073, -73.2628),36 "repentigny": (45.7402, -73.4598),37 "drummondville": (45.8838, -72.4843),38 "granby": (45.4000, -72.7333),39 "saint-hyacinthe": (45.6307, -72.9571),40 "rimouski": (48.4490, -68.5233),41 "victoriaville": (46.0500, -71.9667),42 "saint-georges": (46.1167, -70.6667),43 "matane": (48.8500, -67.5333),44 "mont-tremblant": (46.2103, -74.5847),45}4647# ---------------------------------------------------------------------------48# Haversine helpers49# ---------------------------------------------------------------------------5051def _to_radians(lat: np.ndarray, lon: np.ndarray) -> np.ndarray:52 """Convert lat/lon arrays to radians, shape (n, 2)."""53 return np.column_stack([np.radians(lat), np.radians(lon)])545556def _haversine_vector(lat1: np.ndarray, lon1: np.ndarray, lat2: float, lon2: float) -> np.ndarray:57 """Vectorised haversine distance in metres."""58 lat1_r, lon1_r = np.radians(lat1), np.radians(lon1)59 lat2_r, lon2_r = np.radians(lat2), np.radians(lon2)60 dlat = lat2_r - lat1_r61 dlon = lon2_r - lon1_r62 a = np.sin(dlat / 2) ** 2 + np.cos(lat1_r) * np.cos(lat2_r) * np.sin(dlon / 2) ** 263 return 2 * EARTH_RADIUS_M * np.arcsin(np.sqrt(a))646566# ---------------------------------------------------------------------------67# BallTree construction68# ---------------------------------------------------------------------------6970def _build_balltree(lat: np.ndarray, lon: np.ndarray) -> BallTree:71 """Build a BallTree from latitude/longitude in radians using haversine metric."""72 coords_rad = _to_radians(lat, lon)73 return BallTree(coords_rad, metric="haversine")747576# ---------------------------------------------------------------------------77# Distance to city centre78# ---------------------------------------------------------------------------7980def compute_distance_to_center(81 df: pd.DataFrame,82 city_centers_dict: dict[str, tuple[float, float]] | None = None,83) -> pd.Series:84 """Compute haversine distance (m) from each listing to its city centre."""85 centers = city_centers_dict or DEFAULT_CITY_CENTERS8687 distances = pd.Series(np.nan, index=df.index, dtype=float)88 city_col = df["city"].str.lower().str.strip()8990 for city_key, (clat, clon) in centers.items():91 mask = city_col == city_key92 if mask.any():93 d = _haversine_vector(94 df.loc[mask, "lat"].values,95 df.loc[mask, "long"].values,96 clat,97 clon,98 )99 distances.loc[mask] = d100101 # For cities not in the dictionary, use the median of their own coordinates as centre102 missing_mask = distances.isna()103 if missing_mask.any():104 for city_name, grp in df.loc[missing_mask].groupby(city_col[missing_mask]):105 clat = grp["lat"].median()106 clon = grp["long"].median()107 d = _haversine_vector(grp["lat"].values, grp["long"].values, clat, clon)108 distances.loc[grp.index] = d109110 return distances111112113# ---------------------------------------------------------------------------114# Local density115# ---------------------------------------------------------------------------116117def compute_local_density(df: pd.DataFrame, radius_m: float = 1000.0) -> pd.Series:118 """Count of other listings within radius_m metres using BallTree."""119 valid = df[["lat", "long"]].notna().all(axis=1)120 sub = df.loc[valid]121 lat = sub["lat"].values122 lon = sub["long"].values123 tree = _build_balltree(lat, lon)124125 radius_rad = radius_m / EARTH_RADIUS_M126 coords_rad = _to_radians(lat, lon)127 counts = tree.query_radius(coords_rad, r=radius_rad, count_only=True)128129 density = np.maximum(counts - 1, 0)130 logger.info("Local density (r=%dm): mean=%.1f, max=%d", radius_m, density.mean(), density.max())131 result = pd.Series(np.nan, index=df.index, name=f"local_density_{int(radius_m)}m")132 result.loc[valid] = density133 return result134135136# ---------------------------------------------------------------------------137# Local average price138# ---------------------------------------------------------------------------139140def compute_local_avg_price(df: pd.DataFrame, radius_m: float = 1000.0) -> pd.Series:141 """Average price of other listings within radius_m metres."""142 coord_valid = df[["lat", "long"]].notna().all(axis=1)143 sub = df.loc[coord_valid]144 lat = sub["lat"].values145 lon = sub["long"].values146 prices = sub["price_numeric"].values.copy() if "price_numeric" in sub.columns else np.full(len(sub), np.nan)147 tree = _build_balltree(lat, lon)148149 radius_rad = radius_m / EARTH_RADIUS_M150 coords_rad = _to_radians(lat, lon)151 indices_list = tree.query_radius(coords_rad, r=radius_rad)152153 avg_prices = np.full(len(sub), np.nan)154 for i, neighbours in enumerate(indices_list):155 others = neighbours[neighbours != i]156 neighbour_prices = prices[others]157 valid = neighbour_prices[~np.isnan(neighbour_prices)]158 if len(valid) > 0:159 avg_prices[i] = valid.mean()160161 result = pd.Series(np.nan, index=df.index, name=f"local_avg_price_{int(radius_m)}m")162 result.loc[coord_valid] = avg_prices163 logger.info("Local avg price (r=%dm): computed for %d listings", radius_m, np.sum(~np.isnan(avg_prices)))164 return result165166167# ---------------------------------------------------------------------------168# Spatial lag (KNN)169# ---------------------------------------------------------------------------170171def compute_spatial_lag(df: pd.DataFrame, col: str, k: int = 5) -> pd.Series:172 """K-nearest-neighbour spatial lag of any numeric column."""173 coord_valid = df[["lat", "long"]].notna().all(axis=1)174 sub = df.loc[coord_valid]175 lat = sub["lat"].values176 lon = sub["long"].values177 values = sub[col].values.astype(float)178 tree = _build_balltree(lat, lon)179180 coords_rad = _to_radians(lat, lon)181 dists, indices = tree.query(coords_rad, k=k + 1)182183 lag = np.full(len(sub), np.nan)184 for i in range(len(sub)):185 nbr_idx = indices[i, 1:]186 nbr_vals = values[nbr_idx]187 valid = nbr_vals[~np.isnan(nbr_vals)]188 if len(valid) > 0:189 lag[i] = valid.mean()190191 result = pd.Series(np.nan, index=df.index, name=f"spatial_lag_{col}_k{k}")192 result.loc[coord_valid] = lag193 logger.info("Spatial lag of '%s' (k=%d): computed for %d listings", col, k, np.sum(~np.isnan(lag)))194 return result195196197# ---------------------------------------------------------------------------198# Spatial clusters199# ---------------------------------------------------------------------------200201def create_spatial_clusters(df: pd.DataFrame, n_clusters: int = 20) -> pd.Series:202 """K-means clustering on lat/long coordinates."""203 coords = df[["lat", "long"]].values204 valid_mask = ~np.isnan(coords).any(axis=1)205206 labels = pd.Series(-1, index=df.index, dtype=int, name="spatial_cluster")207 if valid_mask.sum() < n_clusters:208 logger.warning("Too few valid coordinates (%d) for %d clusters", valid_mask.sum(), n_clusters)209 return labels210211 kmeans = KMeans(n_clusters=n_clusters, random_state=42, n_init=10)212 kmeans.fit(coords[valid_mask])213 labels.iloc[np.where(valid_mask)[0]] = kmeans.labels_214215 logger.info("Spatial clusters: %d clusters from %d valid points", n_clusters, valid_mask.sum())216 return labels217218219# ---------------------------------------------------------------------------220# Combined distance features221# ---------------------------------------------------------------------------222223def compute_distance_features(224 df: pd.DataFrame,225 city_centers_dict: dict[str, tuple[float, float]] | None = None,226) -> pd.DataFrame:227 """Compute the standard suite of distance-based features."""228 out = df.copy()229230 out["distance_to_center"] = compute_distance_to_center(out, city_centers_dict)231 out["log_distance_to_center"] = np.log1p(out["distance_to_center"])232233 out["local_density_500m"] = compute_local_density(out, radius_m=500.0)234 out["local_density_1km"] = compute_local_density(out, radius_m=1000.0)235236 out["local_avg_price_500m"] = compute_local_avg_price(out, radius_m=500.0)237 out["local_avg_price_1km"] = compute_local_avg_price(out, radius_m=1000.0)238239 out["spatial_cluster"] = create_spatial_clusters(out)240241 logger.info("Distance features complete: %d new columns", out.shape[1] - df.shape[1])242 return out243244245# ---------------------------------------------------------------------------246# Spatial weights matrix247# ---------------------------------------------------------------------------248249def build_spatial_weights(df: pd.DataFrame, k: int = 5) -> np.ndarray:250 """Build a row-standardised KNN spatial weights matrix (n x n)."""251 lat = df["lat"].values252 lon = df["long"].values253 tree = _build_balltree(lat, lon)254 coords_rad = _to_radians(lat, lon)255256 n = len(df)257 W = np.zeros((n, n), dtype=float)258259 _, indices = tree.query(coords_rad, k=k + 1)260 for i in range(n):261 nbrs = indices[i, 1:] # exclude self262 W[i, nbrs] = 1.0263264 # Row-standardise265 row_sums = W.sum(axis=1, keepdims=True)266 row_sums[row_sums == 0] = 1.0267 W = W / row_sums268269 logger.info("Spatial weights matrix built: %dx%d, k=%d", n, n, k)270 return W271272273# ---------------------------------------------------------------------------274# Moran's I275# ---------------------------------------------------------------------------276277def compute_morans_i(278 df: pd.DataFrame,279 col: str,280 w: np.ndarray | None = None,281 k: int = 5,282) -> dict[str, float]:283 """Compute Moran's I statistic for spatial autocorrelation."""284 y = df[col].values.astype(float)285 valid_mask = ~np.isnan(y)286287 if valid_mask.sum() < 10:288 logger.warning("Too few valid observations (%d) for Moran's I", valid_mask.sum())289 return {"morans_i": np.nan, "expected_i": np.nan, "z_score": np.nan}290291 if w is None:292 w = build_spatial_weights(df, k=k)293294 # Restrict to valid observations295 idx = np.where(valid_mask)[0]296 y_sub = y[idx]297 w_sub = w[np.ix_(idx, idx)]298299 # Re-row-standardise the sub-matrix300 row_sums = w_sub.sum(axis=1, keepdims=True)301 row_sums[row_sums == 0] = 1.0302 w_sub = w_sub / row_sums303304 n = len(y_sub)305 y_mean = y_sub.mean()306 z = y_sub - y_mean307308 numerator = float(n * z @ w_sub @ z)309 denominator = float(w_sub.sum() * (z ** 2).sum())310311 if denominator == 0:312 return {"morans_i": np.nan, "expected_i": np.nan, "z_score": np.nan}313314 I = numerator / denominator315 expected_I = -1.0 / (n - 1)316317 # Variance under normality assumption318 S0 = w_sub.sum()319 S1 = 0.5 * ((w_sub + w_sub.T) ** 2).sum()320 S2 = ((w_sub.sum(axis=1) + w_sub.sum(axis=0)) ** 2).sum()321 k2 = (n * (z ** 4).sum()) / ((z ** 2).sum() ** 2)322323 A = n * ((n ** 2 - 3 * n + 3) * S1 - n * S2 + 3 * S0 ** 2)324 B = k2 * ((n ** 2 - n) * S1 - 2 * n * S2 + 6 * S0 ** 2)325 C = (n - 1) * (n - 2) * (n - 3) * S0 ** 2326327 if C == 0:328 z_score = np.nan329 else:330 var_I = (A - B) / C - expected_I ** 2331 if var_I > 0:332 z_score = float((I - expected_I) / np.sqrt(var_I))333 else:334 z_score = np.nan335336 result = {"morans_i": float(I), "expected_i": float(expected_I), "z_score": z_score}337 logger.info("Moran's I for '%s': I=%.4f, E[I]=%.4f, z=%.2f", col, I, expected_I, z_score if not np.isnan(z_score) else 0.0)338 return result339