How We Made This Map

Rows are spaced 0.0045 degrees of latitude apart, about 500 meters. A degree of longitude shrinks with latitude by a factor of cos φ, so columns are spaced 0.0045 / cos(φ) × 2/√3 degrees to make the triangles near-equilateral on the ground, with every other row offset by half a column.

Deciding which points are in water uses the same shoreline geometry the map draws: lake and bay polygons from OpenStreetMap, plus the Sound assembled from its mapped coastline. Each point is tested with the standard ray-casting trick: draw a line from the point out to the edge of the map and count how many shoreline borders it crosses. An odd count means the point is in water. 883 points survive on land.

The mesh behind the X-ray button connects each surviving point to its neighbors within about 1.6 grid steps. If a connecting line touches water, it is deleted, with one exception: lines that pass within 350 meters of a named bridge are kept, because there the crossing exists. Each bridge keeps only its two shortest crossings so the deck reads as a deck and not a fan of diagonals.

The street network becomes a directed graph with about 136,000 nodes and 248,000 edges. An edge's weight is the time it takes to drive it: the segment's length divided by its speed limit. Length comes from the haversine formula. The speed limit comes from the road's maxspeed tag in OpenStreetMap, and when a road has no tag, a default for its class fills in, 25 mph for a residential street, 60 for a motorway. A one-way street contributes an edge in one direction only.

Street data is full of stranded fragments: parking aisles, gated stubs, driveways. Routing between them would report infinite times, so the router first finds the largest set of mutually reachable nodes (Kosaraju's algorithm) and stays inside it. That set holds 98% of the network. Each sample point then snaps to its nearest node in the set, and one run of Dijkstra's algorithm gives the time from that point to every other node. 883 runs fill the whole table in a few seconds.

Traffic enters as a per-edge multiplier. TomTom reports, for the road nearest a sample point, how long it currently takes compared to free flow. A reading of 1.4 means 40% slower than an open road. Each edge of the graph averages the four nearest readings within 1.5 km, weighting closer readings more, and the multiplier is capped at 4× so one stalled sensor cannot distort a whole corridor. Where no reading is close, the road keeps its free-flow time. Nothing else in the pipeline changes with traffic; the only difference downstream is the numbers in the table.

Classical MDS is the one-step start. Square the time matrix entrywise to get D⁽²⁾, then double-center it: B = -½ J D⁽²⁾ J with J = I - (1/n) 1 1ᵀ. If the times were distances between points in some space, B would be exactly the matrix of dot products of those points, and its top two eigenvectors, scaled by the square roots of their eigenvalues, would be the best 2D coordinates. The eigenvectors come from power iteration: multiply a random vector by B repeatedly until it stops turning, then deflate and repeat for the second.

Drive times are not distances in any flat space, so that answer is only a good start. SMACOF (stress majorization) refines it. The target is raw stress, Σ (dᵢⱼ(X) - δᵢⱼ)²: drawn distance minus measured time, squared, summed over pairs. Each round applies the Guttman transform: point i moves to the average, over partners j, of the position that would make the pair's drawn distance exactly δᵢⱼ, which is xⱼ + δᵢⱼ (xᵢ - xⱼ) / dᵢⱼ. Majorization theory guarantees stress never increases; 150 rounds are enough for this data.

The transit map changes this step in two ways. First, each pair's error is weighted by 1/δ², so a 10-minute pair drawn 5 minutes wrong costs as much as an 80-minute pair drawn 40 minutes wrong. Transit times run from minutes to nearly two hours, and without the weighting the longest pairs dominate and bend the neighborhoods out of shape. Second, the starting arrangement is geography itself, scaled into time units, instead of classical MDS. A network built around a few spines leaves whole districts free to swivel without changing any pairwise distance, and the settling keeps whatever swivel its starting point had. Starting from geography, there is none to keep.

This is orthogonal Procrustes with uniform scaling. Center the layout X and the geographic reference G on their means, then accumulate the four cross-covariance sums sₓₓ = Σ gₓ xₓ, sₓₔ = Σ gₓ xₔ, and so on. In 2D the rotation that best matches X to G has the closed form θ = atan2(sₔₓ - sₓₔ, sₓₓ + sₔₔ). When the determinant of the cross-covariance is negative the best fit includes a mirror flip, so the layout is reflected and the angle re-derived, keeping whichever version correlates better with geography.

The scale is the least-squares ratio Σ ⟨R xᵢ, gᵢ⟩ / Σ |xᵢ|². Every transform used here (rotation, reflection, translation, uniform scale) changes all distances by the same factor. The map exists to show relative distortion, and relative distortion is exactly what these transforms leave alone.

The reported number is Kruskal's stress-1: √(Σ(dᵢⱼ - δᵢⱼ)² / Σδᵢⱼ²), the residual as a fraction of the total. The free-flow map fits at 0.124, the Friday evening map at 0.111. Congestion makes the times more uniform, and more uniform distances flatten into two dimensions with less error. The transit map fits at 0.377 by this measure. That is what it should say: a network of spines and transfers is much further from flat than a street grid. The number reports the unweighted residual even though the transit settling optimizes the weighted one. Each anchor's dot radius scales with its share of the residual.

Here is the bridge example with exact numbers: three points on a circle of radius 10 are pairwise farthest when they form an equilateral triangle, with side 10√3 ≈ 17.3. The measured 20s exceed that bound, so no flat drawing satisfies all six numbers; the excess spreads across the neighborhood as stress. Drive-time data is full of such conflicts, which is also why the embedding targets 2D and not, say, 3D: the goal is a readable map, and the residual is shown rather than hidden.