Spatial Embeddings for Location Similarity

“Find me somewhere like this” is a nearest-neighbour query over a vector per place. Building that vector is almost entirely a preprocessing problem, and the two steps everyone skips — the transform and the standardisation — change the answer more than the choice of distance metric does.

Choosing the unit of location

Before any feature exists you have to decide what a location is. The three usual choices behave differently and are not interchangeable.

• A fixed grid cell. H3 at resolution 8 or 9, or a projected square grid. Every unit has the same area, which means counts are directly comparable and you can skip an area normalisation. It cuts through real boundaries — half a shopping centre in one cell, half in the next.
• An administrative polygon. Census tracts, output areas, postcode sectors. They align with the demographic data you want, which is why they are tempting, and their areas vary by orders of magnitude between a city centre and a rural district, so every count must become a density first.
• A point with a radius. One vector per store or per candidate site, built from everything within 500 m. This is what you want if the question is site selection rather than mapping, and it produces overlapping units, so treating them as independent samples in any downstream model is wrong.

The grid is usually the right default because equal area removes an entire class of error, and H3’s hexagons add the property that all neighbours are equidistant, which matters as soon as you smooth anything across cells.

The three feature families

A location vector is a concatenation. Each block has its own units and its own pathologies.

Points of interest. Counts per category within the unit or within a radius: cafes, schools, offices, transit stops, parking. From OpenStreetMap this is a tag query; from a commercial POI dataset it is a category taxonomy with a few hundred leaves. Category counts are sparse, heavy-tailed and strongly correlated with each other — cafes and restaurants and bars move together — so the effective dimensionality is far below the nominal one. Rolling a taxonomy up to twenty or thirty categories usually loses nothing that the distance metric was using.

Demographics. Population density, age distribution, household income, tenure, education. These arrive on administrative polygons and have to be reallocated onto your unit, which is areal interpolation and is where a quiet error enters: distributing tract population uniformly by area is wrong wherever the tract contains a park or a reservoir. Weighting by building footprints or by night-lights is the standard fix and it is the same problem as estimating population inside an arbitrary polygon.

Mobility. Visit counts by hour of day and day of week, origin-destination flows, dwell time. This is the block that distinguishes places the other two call identical: an office district and a residential district can have similar POI counts and similar resident demographics while having exactly opposite hourly profiles. Represent it as a normalised 24- or 168-dimensional profile — the shape, not the volume — and keep total volume as a separate scalar, because otherwise a busy place and a quiet place with the same rhythm look different when you wanted them to look alike.

Why the scaling decides the answer

Concatenate raw features and the vector is dominated by whichever column has the largest numbers. Population per cell might run to thousands, income to tens of thousands, cafe count to single digits. Euclidean distance on that is a distance in income with rounding noise from everything else.

Two transforms, in this order, do most of the work. First log1p on every count and every density, because these distributions are heavy-tailed and one central business district otherwise sets the scale for the whole dataset. Using log1p rather than log keeps zeros finite, and zeros are everywhere in POI data. Second, standardise each column to zero mean and unit variance across the whole study area, so that a one-standard- deviation difference in cafe count counts the same as a one-standard-deviation difference in income.

Then choose the metric deliberately. Cosine similarity compares composition and ignores magnitude, so a small town centre and a large city centre with the same mix come out similar. Euclidean distance on standardised features keeps magnitude, so they do not. Neither is correct in general; they answer different questions, and the way to decide is to state which pairs you want called similar before you look at any output. The tradeoff is the same one described in vector similarity metrics.

A worked similarity

Four features, three cells. Raw values first, then the same values after log1p and standardisation. The standardisation constants here are stated as assumptions — in a real run they come from the full study area, not from three rows.

raw                     cafes  offices  pop_density  evening_share
  A  high-street           38      12         6,400          0.41
  B  office district        9      74         1,100          0.06
  C  suburban centre       21       8         4,900          0.37

log1p, then z-scored using study-area mean/sd (assumed):
  A   [ 1.20,  -0.35,   0.88,   0.75 ]
  B   [-0.90,   1.42,  -1.05,  -1.30 ]
  C   [ 0.55,  -0.62,   0.60,   0.61 ]

Euclidean A-C:
  d^2 = (1.20-0.55)^2 + (-0.35+0.62)^2 + (0.88-0.60)^2 + (0.75-0.61)^2
      = 0.4225 + 0.0729 + 0.0784 + 0.0196 = 0.5934
  d   = 0.770

Euclidean A-B:
  d^2 = (1.20+0.90)^2 + (-0.35-1.42)^2 + (0.88+1.05)^2 + (0.75+1.30)^2
      = 4.410 + 3.133 + 3.725 + 4.203 = 15.471
  d   = 3.933

A is five times closer to C than to B, which is the answer you wanted. Now run the same comparison on the raw numbers and pop_density alone contributes (6400 − 4900)² = 2.25 million to the A–C distance while evening_share contributes 0.0016. The evening share — the one feature that actually separates a high street from an office district — is invisible. That is not a subtle degradation; it is the whole model reduced to one column, and it happens silently.

The same reasoning governs how many features to keep. Correlated blocks inflate the dimensions without adding information and the distances concentrate, which is the standard high-dimensional behaviour described in embedding dimensions. Thirty well-chosen, decorrelated features usually beat three hundred raw category counts.

Spatial autocorrelation breaks your split

Nearby places resemble each other. That is Tobler’s first law of geography and it is the reason a random train/test split over spatial units reports a score you cannot reproduce anywhere new.

Split 10,000 cells at random and the test cell is usually adjacent to several training cells that share its shopping centre, its transit stop and its census tract. The model does not need to have learned anything transferable to score well; it needs only to have memorised the neighbourhood. Held-out error is then an estimate of interpolation skill, and the question you actually asked was extrapolation skill: will this work in a city where I have no data?

The fix is spatial blocking. Partition the study area into contiguous blocks much larger than the correlation range, then assign whole blocks to folds. A useful sanity check is to hold out an entire city and see how far the score falls; the gap between the random split and the city-level split is the size of the leakage, and it is often large. Report the blocked number. Anything else is a claim about places you already have.