diff --git a/copenhagen/scripts/render.py b/copenhagen/scripts/render.py index ed3794a..ce63d53 100644 --- a/copenhagen/scripts/render.py +++ b/copenhagen/scripts/render.py @@ -169,8 +169,16 @@ def classify_station(r): return None +DISTANCE_CRS = "EPSG:25832" # UTM 32N — accurate meters for Copenhagen + + def cluster_stops(coords, threshold_m): - """Cluster (x, y) coordinates within threshold_m; return list of centroids.""" + """Cluster (x, y) coordinates in EPSG:3857 within threshold_m (real meters). + + Reprojects to UTM 32N for accurate distance computation, clusters, then + returns centroids in EPSG:3857 for plotting. + """ + from pyproj import Transformer from scipy.cluster.hierarchy import fcluster, linkage from scipy.spatial.distance import pdist @@ -178,14 +186,21 @@ def cluster_stops(coords, threshold_m): if len(coords) <= 1: return coords.tolist() if len(coords) else [] - dists = pdist(coords) + # reproject to UTM for accurate distances + to_utm = Transformer.from_crs(TARGET_CRS, DISTANCE_CRS, always_xy=True) + back = Transformer.from_crs(DISTANCE_CRS, TARGET_CRS, always_xy=True) + utm = np.array([to_utm.transform(x, y) for x, y in coords]) + + dists = pdist(utm) links = linkage(dists, method="single") labels = fcluster(links, t=threshold_m, criterion="distance") centroids = [] for label in set(labels): - members = coords[labels == label] - centroids.append(members.mean(axis=0).tolist()) + members = utm[labels == label] + ux, uy = members.mean(axis=0) + cx, cy = back.transform(ux, uy) + centroids.append((cx, cy)) return centroids