Optimizing Geospatial Queries: Finding the K Nearest Weather Stations in Go

We start with a practical scenario: locating the K nearest national weather observation stations to the 768 Creative Industry Park in Haidian, Beijing. The simplest method is to iterate over all candidate stations, compute the distance to each, sort by distance, and pick the top K. Although straightforward to implement, this approach is slow.

How to find the K nearest national observation stations to the 768 Creative Industry Park?

The performance bottleneck arises because distance calculation cannot rely on simple planar formulas; the Earth’s curvature matters for larger separations. The most accurate method is Vincenty's formulae, which involves many trigonometric operations and is computationally expensive.

Vincenty's formulae (see [1]) is shown in the following diagram:

An alternative is the haversine (half‑chord) formula, which is slightly less accurate but sufficient for meteorological data. It still contains multiple trigonometric calls, and even with CPU‑level instruction optimizations, these operations remain slower than basic arithmetic, especially when processing tens of thousands of candidate stations.

Our first optimization adds a coarse bounding‑box filter: only stations whose latitude and longitude fall within a predefined box are subjected to the haversine calculation; others are skipped. This reduces typical query times from around 2,000,000 ns to 4,000,000 ns for slower cases.

Profiling with pprof revealed that the dominant CPU cost stems from memory allocations caused by iterating over a slice of length in the ten‑thousands. To avoid full scans, we need a spatial index.

We evaluated community R‑Tree and Quad‑Tree implementations but found performance issues in some scenarios, such as query times around 1,000,000 ns and lack of built‑in sorting for K‑nearest‑neighbor results, which still required manual distance computation and sorting.

Inspired by map tiling, we devised a custom tile‑based index. Each station is assigned to a tile at a specific zoom level (x, y, z). By examining the nine tiles surrounding the target tile (the central tile plus its eight neighbors), we can limit distance calculations to a small subset.

┌───────────┬───────────┬───────────┐
│           │           │           │
│           │           │           │
│ x-1,y-1,z │ x+0,y-1,z │ x+1,y-1,z │
│           │           │           │
│           │           │           │
├───────────┼───────────┼───────────┤
│           │           │           │
│           │           │           │
│ x-1,y+0,z │ x+0,y+0,z │ x+1,y+0,z │
│           │           │           │
│           │           │           │
├───────────┼───────────┼───────────┤
│           │           │           │
│           │           │           │
│ x-1,y+1,z │ x+0,y+1,z │ x+1,y+1,z │
│           │           │           │
│           │           │           │
└───────────┴───────────┴───────────┘

Because the tile index uses only integer arithmetic after the initial trigonometric conversion, subsequent lookups involve simple addition and subtraction.

Internally we store tile keys in reverse order z,y,x within a map, allowing us to retrieve the nine relevant keys without scanning the entire dataset.

Observational data shows uneven spatial distribution: U.S. stations are relatively uniform, while Chinese stations are denser in the east than the west, and Western Europe is denser than Eastern Europe. To handle this, we employ a two‑level pre‑index: if a high‑resolution tile yields no stations, we fall back to a coarser tile level.

After applying these optimizations, query latency dropped to roughly 900 ns–3,000 ns, and CPU consumption decreased by about 10 %. The design also leaves room for optional R‑Tree indexing for sparse datasets with variable spatial resolution.

References

[1] Vincenty's formulae: https://en.wikipedia.org/wiki/Vincenty%27s_formulae