How to Use H3 Neighbours for k-Ring Smoothing and Buffers

On this page

Problem statement

Counts per hexagon at a fine resolution are mostly noise. Binning the 43,712 GeoNames populated places in the United Kingdom into res-7 cells (about 5 km² each) left 54.1% of the 53,100 land cells empty, and the correlation between a cell and its neighbours was only 0.406 — a map of that surface is a scatter of speckle. The usual remedies are a coarser resolution, which throws away detail everywhere, or kernel density, which needs a projected CRS and a raster. H3 offers a third: every cell knows its neighbours by index arithmetic, so a moving-window average is a dictionary lookup per neighbour.

Measured on the same cells, a k=1 neighbourhood mean cut the empty share to 16.1%, halved the peak from 31 to 15.14, and kept the total within 0.2%. The same neighbour functions give you an approximate buffer — with an error that has to be handled, because a disk of hexagons is not a circle.

Quick answer

grid_disk(cell, k) returns the cell and every cell within k steps. Average over it:

import h3


def smooth(counts, k=1):
    """Mean of each cell and its k-ring neighbours; missing neighbours count as zero."""
    return {
        cell: sum(counts.get(n, 0) for n in h3.grid_disk(cell, k)) / (1 + 3 * k * (k + 1))
        for cell in counts
    }


counts = {h3.latlng_to_cell(51.5074, -0.1278, 8): 12}
print(smooth(counts, k=1))

{'88195da49bfffff': 1.7142857142857142}

Twelve points spread over the seven cells of the disk give 12 / 7. That version divides by the nominal disk size and only visits cells that already have counts; Example 1 fixes both, handling pentagons and the edges of a study area, and it is the one to use on real data.

Measured speed: 100,000 disks at k=1 in 0.26 s, at k=2 in 0.48 s, at k=5 in 1.80 s.

The disk is index arithmetic: no geometry is built to find these nineteen cells.

Disks, rings, grid distance and a full study area

1. Get a disk or a ring

origin = h3.latlng_to_cell(51.5074, -0.1278, 8)       # Charing Cross
for k in range(4):
    print(k, len(h3.grid_disk(origin, k)), len(h3.grid_ring(origin, k)))

grid_disk is filled — everything within k steps, 1 + 3k(k + 1) cells. grid_ring is hollow — exactly k steps away, 6k cells. Use the disk for windows and buffers, the ring when neighbours at different distances need different weights.

In 2,000 random disks the origin came first and ring 1 came next every time, but the h3-py documentation does not promise an order. Ask grid_ring for distances rather than slicing a disk.

2. Expect pentagons to return fewer cells

pentagon = h3.get_pentagons(8)[0]
print(len(h3.grid_disk(pentagon, 1)), len(h3.grid_disk(pentagon, 2)), len(h3.grid_ring(pentagon, 1)))

6 16 5

A pentagon has five neighbours, so its disks hold 6 and 16 cells rather than 7 and 19. Dividing by 1 + 3k(k + 1), as the quick answer does, under-weights those cells. There are only twelve pentagons per resolution and all twelve centres lie at sea, but a global pipeline will meet them.

3. Measure grid distance, and expect it to fail

london = h3.latlng_to_cell(51.5, -0.1, 5)
paris = h3.latlng_to_cell(48.86, 2.35, 5)
new_york = h3.latlng_to_cell(40.7, -74.0, 5)
print(h3.grid_distance(london, paris))
try:
    h3.grid_distance(london, new_york)
except h3.H3FailedError as exc:
    print(type(exc).__name__)

21
H3FailedError

grid_distance counts steps between two cells at the same resolution. It is computed by unfolding the icosahedron faces into flat local IJ coordinates, which cannot be done consistently across a pentagon, where five faces meet, or over long distances — so rather than return a wrong number it raises for cells far apart. It also fails close to a pentagon: in the 3-disk around one pentagon, 360 of 961 cell pairs raised.

4. Build the whole study area, zeros included

A dictionary of occupied cells cannot be smoothed into the gaps — the empty cells are not in it. Fill the region first:

import duckdb
import geopandas as gpd

con = duckdb.connect()
con.execute("set threads = 4")
con.execute("install h3 from community; load h3")
occupied = con.execute("""
    select h3_h3_to_string(h3_latlng_to_cell(lat, lon, 7)) as cell, count(*) as n
    from read_parquet('geonames_all.parquet')
    where country_code = 'GB' and feature_class = 'P'
    group by 1
""").df()

countries = gpd.read_file("zip://ne_10m_admin_0_countries.zip")
uk = countries.loc[countries.ADMIN == "United Kingdom", "geometry"].iloc[0]
study = h3.geo_to_cells(uk, 7)

base = dict.fromkeys(study, 0)
base.update({c: n for c, n in zip(occupied.cell, occupied.n) if c in base})
print(f"{len(study):,} cells in the study area, {sum(v > 0 for v in base.values()):,} occupied")

53,100 cells in the study area, 24,359 occupied

The GeoNames points occupied 25,174 cells in total. The other 815 are cells whose centres fall outside the Natural Earth coastline, which a centre-based polygon fill excludes even though places lie inside them.

5. Smooth with a k-ring mean

With Example 1 applied to base, measured on the 53,100 cells:

raw           zero  54.1%  sd 1.212  max  31.00  sum   42,577  neighbour r 0.406
k=1           zero  16.1%  sd 0.832  max  15.14  sum   42,497  neighbour r 0.892
k=2           zero   6.1%  sd 0.742  max  10.95  sum   42,421  neighbour r 0.961
k=1 weighted  zero  16.1%  sd 0.844  max  16.38  sum   42,517  neighbour r 0.888

All three surfaces took 1.0 s together. The weighted version counts the centre cell twice as heavily as its neighbours, and it is only marginally sharper than the plain k=1 mean.

6. Decide what happens at the edge of the study area

9.2% of the UK cells were missing at least one neighbour — coast, mostly. Averaging over the neighbours that exist (outside="ignore") moved the sum from 42,577 to 42,497; treating outside cells as zero (outside="zero") dropped it to 41,307, because coastal counts leak into the sea.

For counts of things that cannot be offshore, ignore. For a sampled surface that really continues beyond the boundary, fetch the neighbouring region instead of pretending it is empty.

7. Treat a disk as a buffer only with a distance check

A disk is a cheap stand-in for "everything within about this distance", and it is not a circle. At res 8 a k=1 disk reached 969 m in its narrowest direction and 1,402 m in its widest (Example 2). Tested against true great-circle distance for a 1 km radius at res 9 around 300 UK places, k=3 still missed 14 of 287 true neighbours, and k=4 missed none — but 275 of its 562 candidates were beyond 1 km (Example 2). Use the disk to find candidates, then refine. Choosing k for a given distance, with a rule verified on 140,000 adversarial pairs, is covered in How to join two point datasets on an H3 index.

The neighbour correlation of a smoothed surface is mostly the smoother, so test autocorrelation on the raw counts.

Weighted smoothing and how round a disk is

The data used on this page

Country outlines from Natural Earth's 1:10m admin-0 layer; every feature in the GeoNames allCountries dump as a point. GeoNames is edited daily, so counts drift a little between downloads.

import os
import zipfile

import pandas as pd
import requests


def fetch(url, path):
    if not os.path.exists(path):
        with requests.get(url, stream=True, timeout=600) as reply:
            reply.raise_for_status()
            with open(path, "wb") as out:
                for chunk in reply.iter_content(1 << 20):
                    out.write(chunk)


fetch("https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_0_countries.zip",
      "ne_10m_admin_0_countries.zip")
if not os.path.exists("geonames_all.parquet"):          # every GeoNames feature, about 400 MB zipped
    fetch("https://download.geonames.org/export/dump/allCountries.zip", "allCountries.zip")
    with zipfile.ZipFile("allCountries.zip") as z, z.open("allCountries.txt") as f:
        pd.read_csv(f, sep="\t", header=None, usecols=[4, 5, 6, 8], quoting=3,
                    names=["lat", "lon", "feature_class", "country_code"],
                    dtype={"feature_class": str, "country_code": str}).to_parquet("geonames_all.parquet")
print(len(pd.read_parquet("geonames_all.parquet", columns=["lat"])))

13472317

Example 1 — smoothing a cell surface with weights and edge handling

import h3
import numpy as np
import pandas as pd


def kring_smooth(values, k=1, weights=None, outside="ignore"):
    """Smooth a {cell: value} surface over H3 rings.

    weights: one weight per ring distance 0..k (default: all 1).
    outside: 'ignore' averages over neighbours that exist in `values`;
             'zero' treats neighbours outside the study area as zero.
    """
    weights = weights or [1.0] * (k + 1)
    out = {}
    for cell in values:
        num = den = 0.0
        for d in range(k + 1):
            ring = [cell] if d == 0 else h3.grid_ring(cell, d)
            for n in ring:
                if n in values:
                    num += weights[d] * values[n]
                    den += weights[d]
                elif outside == "zero":
                    den += weights[d]
        out[cell] = num / den
    return pd.Series(out)


def neighbour_correlation(series):
    """Correlation between each cell's value and each of its ring-1 neighbours'."""
    a, b = [], []
    for cell, v in series.items():
        for n in h3.grid_ring(cell, 1):
            if n in series.index:
                a.append(v)
                b.append(series[n])
    return np.corrcoef(a, b)[0, 1]


import time

raw = pd.Series(base)
start = time.perf_counter()
surfaces = {"raw": raw, "k=1": kring_smooth(base, 1), "k=2": kring_smooth(base, 2),
            "k=1 weighted": kring_smooth(base, 1, weights=[2.0, 1.0])}
took = time.perf_counter() - start
for name, s in surfaces.items():
    print(f"{name:13} zero {(s == 0).mean():5.1%}  sd {s.std():.3f}  max {s.max():6.2f}  sum {s.sum():8,.0f}  "
          f"neighbour r {neighbour_correlation(s):.3f}")
edge = sum(any(n not in base for n in h3.grid_ring(c, 1)) for c in base) / len(base)
print(f"three surfaces in {took:.1f} s; {edge:.1%} of cells miss a neighbour; "
      f"sum with outside='zero' {kring_smooth(base, 1, outside='zero').sum():,.0f}")
occupied_all = set(occupied.cell)
print(f"{len(occupied_all):,} occupied cells in all, {len(occupied_all - set(base)):,} outside the centre fill")

raw           zero 54.1%  sd 1.212  max  31.00  sum   42,577  neighbour r 0.406
k=1           zero 16.1%  sd 0.832  max  15.14  sum   42,497  neighbour r 0.892
k=2           zero  6.1%  sd 0.742  max  10.95  sum   42,421  neighbour r 0.961
k=1 weighted  zero 16.1%  sd 0.844  max  16.38  sum   42,517  neighbour r 0.888
three surfaces in 0.8 s; 9.2% of cells miss a neighbour; sum with outside='zero' 41,307
25,174 occupied cells in all, 815 outside the centre fill

The function divides by the weights that actually took part, so a pentagon's five neighbours and a coastal cell's surviving neighbours are averaged correctly. That is the table in step 5, and neighbour_correlation produced its last column.

Example 2 — how round is a disk?

from pyproj import Transformer
from shapely import union_all
from shapely.geometry import Point, Polygon


def disk_radii(lat, lng, res, k):
    """Inner and outer radius, in metres, of the k-disk around the cell containing a point."""
    cell = h3.latlng_to_cell(lat, lng, res)
    clat, clng = h3.cell_to_latlng(cell)
    to_m = Transformer.from_crs(
        "EPSG:4326", f"+proj=aeqd +lat_0={clat} +lon_0={clng} +units=m", always_xy=True
    )
    hexes = []
    for n in h3.grid_disk(cell, k):
        lngs, lats = zip(*[(b, a) for a, b in h3.cell_to_boundary(n)])
        hexes.append(Polygon(zip(*to_m.transform(lngs, lats))))
    outline = union_all(hexes).exterior
    inner = outline.distance(Point(0, 0))
    outer = max(Point(xy).distance(Point(0, 0)) for xy in outline.coords)
    return inner, outer


for res, k in [(8, 1), (8, 3), (8, 6), (9, 5)]:
    inner, outer = disk_radii(51.5074, -0.1278, res, k)
    print(f"res {res} k={k}: inner {inner:6,.0f} m  outer {outer:6,.0f} m  ratio {outer / inner:.3f}")

res 8 k=1: inner    969 m  outer  1,402 m  ratio 1.447
res 8 k=3: inner  2,421 m  outer  3,214 m  ratio 1.328
res 8 k=6: inner  4,598 m  outer  5,947 m  ratio 1.293
res 9 k=5: inner  1,453 m  outer  1,912 m  ratio 1.316

The inner radius is the guarantee: everything closer than that to the cell's centre is inside the disk. The outer radius is the worst over-reach.

And as a buffer on real points — every UK place within 1 km of one of 300 sampled UK places, against the res-9 disk around the sampled place's cell:

from scipy.spatial import cKDTree

places = con.execute("""select lat, lon from read_parquet('geonames_all.parquet')
                        where country_code = 'GB' and feature_class = 'P'""").df()
xyz = np.column_stack([np.cos(np.radians(places.lat)) * np.cos(np.radians(places.lon)),
                       np.cos(np.radians(places.lat)) * np.sin(np.radians(places.lon)),
                       np.sin(np.radians(places.lat))])
cells9 = [h3.latlng_to_cell(a, b, 9) for a, b in zip(places.lat, places.lon)]
origins = np.random.default_rng(0).choice(len(places), 300, replace=False)
tree = cKDTree(xyz)
true_pairs = [(o, j) for o in origins for j in tree.query_ball_point(xyz[o], 1.0 / 6371.0088) if j != o]
for k in (3, 4):
    missed = sum(cells9[j] not in set(h3.grid_disk(cells9[o], k)) for o, j in true_pairs)
    print(f"k={k}: missed {missed} of {len(true_pairs):,} places truly within 1 km")
by_cell = {}
for j, c in enumerate(cells9):
    by_cell.setdefault(c, []).append(j)
candidates = [(o, j) for o in origins for c in h3.grid_disk(cells9[o], 4) for j in by_cell.get(c, []) if j != o]
far = sum(h3.great_circle_distance((places.lat[o], places.lon[o]), (places.lat[j], places.lon[j]), unit="km") > 1
          for o, j in candidates)
print(f"k=4 candidates: {len(candidates):,}, of which {far:,} are beyond 1 km")

k=3: missed 14 of 287 places truly within 1 km
k=4: missed 0 of 287 places truly within 1 km
k=4 candidates: 562, of which 275 are beyond 1 km

Example 3 — a distance that does not raise

def cell_distance_km(a, b):
    """Grid steps when H3 can compute them, great-circle distance otherwise."""
    try:
        steps = h3.grid_distance(a, b)
    except h3.H3FailedError:
        steps = None
    km = h3.great_circle_distance(h3.cell_to_latlng(a), h3.cell_to_latlng(b), "km")
    return steps, round(km, 1)


print(cell_distance_km(london, paris))
print(cell_distance_km(london, new_york))

(21, 340.8)
(None, 5582.3)

Grid steps are useful for "within k cells" logic; kilometres are what most questions actually ask. Returning both, and never raising, keeps a batch job alive when one pair straddles a pentagon.

Example 4 — the remaining measurements

import math
import random
import time

for k in (1, 2, 5):
    cells = [h3.latlng_to_cell(random.uniform(-60, 60), random.uniform(-180, 180), 7) for _ in range(100_000)]
    start = time.perf_counter()
    for c in cells:
        h3.grid_disk(c, k)
    print(f"100,000 disks at k={k}: {time.perf_counter() - start:.2f} s")

random.seed(1)
ordered = 0
for _ in range(2000):
    c = h3.latlng_to_cell(random.uniform(-60, 60), random.uniform(-180, 180), 7)
    disk = h3.grid_disk(c, 2)
    ordered += disk[0] == c and set(disk[1:7]) == set(h3.grid_ring(c, 1))
print(f"origin first and ring 1 next in {ordered} of 2,000 disks")

london = h3.latlng_to_cell(51.5074, -0.1278, 7)
d = [h3.great_circle_distance(h3.cell_to_latlng(london), h3.cell_to_latlng(n), unit="km") for n in h3.grid_ring(london, 1)]
print(f"res 7 neighbour centres in London: {min(d):.3f} to {max(d):.3f} km ({max(d) / min(d) - 1:.1%} apart)")
dlat, dlng = 180 / 2 ** 15, 360 / 2 ** 15                 # a precision-6 geohash cell
km = 6371.0088 * math.pi / 180
edge_ns, edge_ew = dlat * km, dlng * km * math.cos(math.radians(51.5074))
print(f"geohash precision 6 in London: neighbours {min(edge_ns, edge_ew):.3f} to {math.hypot(edge_ns, edge_ew):.3f} km "
      f"(ratio {math.hypot(edge_ns, edge_ew) / min(edge_ns, edge_ew):.3f})")

pentagon = next(p for p in h3.get_pentagons(7))
near = h3.grid_disk(pentagon, 3)
failed = 0
for a in near:
    for b in near:
        try:
            h3.grid_distance(a, b)
        except h3.H3BaseException:
            failed += 1
print(f"grid_distance near a pentagon: {failed} of {len(near) ** 2:,} pairs raised")
try:
    h3.grid_distance(h3.latlng_to_cell(51.5074, -0.1278, 5), h3.latlng_to_cell(40.7128, -74.006, 5))
except h3.H3BaseException as exc:
    print("London to New York at res 5:", type(exc).__name__)

100,000 disks at k=1: 0.26 s
100,000 disks at k=2: 0.49 s
100,000 disks at k=5: 1.88 s
origin first and ring 1 next in 2000 of 2,000 disks
res 7 neighbour centres in London: 2.243 to 2.426 km (8.2% apart)
geohash precision 6 in London: neighbours 0.611 to 0.975 km (ratio 1.597)
grid_distance near a pentagon: 360 of 961 pairs raised
London to New York at res 5: H3FailedError

Good smoothing windows and what they cost

Why hexagon neighbours make a good smoothing window

Every H3 neighbour shares an edge, and neighbour centres sit at nearly the same distance. Measured at res 7 in London, the six centre distances ranged from 2.243 to 2.426 km — within 8.2%. The eight neighbours of a geohash cell at precision 6 ranged from 0.611 to 0.975 km, a ratio of 1.597, because corner neighbours are further away than edge neighbours.

A square window therefore weights diagonal directions differently from straight ones. A hexagonal ring does not need a correction.

Why smoothing trades location for stability

Averaging seven cells divides the sampling noise, which is why the standard deviation fell from 1.212 to 0.832 at k=1 and the empty share from 54.1% to 16.1%. The price is location: the single busiest raw cell (31 places, in the City of London) was diluted to a mean of its neighbourhood, and the k=1 maximum moved to a cell about 2 km west.

At k=2 each value averages nineteen res-7 cells, about 98 km² or nearly three res-6 cells. The advantage over simply choosing a coarser resolution is that the window slides: every small hexagon gets its own neighbourhood, so the map keeps its fine grain.

Why the neighbour correlation is not a finding

Smoothing makes neighbours share inputs. Two adjacent k=1 cells average overlapping windows, so their values are correlated by construction: 0.406 raw, 0.892 at k=1, 0.961 at k=2.

That correlation belongs to the smoother, not the data. Any statistic about spatial dependence — Moran's I, a hot-spot test, a regression residual check — must run on the raw counts, with the neighbour structure supplied as weights.

Why a disk is not a circle

A disk of hexagons is itself a rough hexagon with a stepped edge. Its inner radius is set by the notches between outer cells and its outer radius by their corners, so the ratio is large for small k — 1.447 at k=1 — and shrinks towards the shape of a hexagon as k grows: 1.328 at k=3, 1.293 at k=6.

On top of that, a point is not at its cell's centre, so the disk around the point's cell is offset from the point by up to about one edge length. Both effects are why a disk buffer needs a margin of extra rings and a true-distance refinement.

Near a pentagon, 360 of 961 grid_distance calls raised in one measured neighbourhood.

Gaps, the antimeridian and quadratic cost

Smoothing a dictionary of occupied cells never fills gaps. Build the study area with a polygon fill first, zeros included.
Disks cross the antimeridian without special handling. A res-7 disk near Fiji returned its seven cells with centres on both sides of ±180°; only drawing them needs care.
Cost grows with the square of k. A k=5 disk is 91 cells; 100,000 of them took 1.80 s against 0.26 s at k=1.
The total drifts at the edges. 42,577 became 42,497 averaging over existing neighbours, and 41,307 when outside cells counted as zero.
Centre weighting changes little. Doubling the centre weight kept the maximum at 16.38 against 15.14 unweighted.
All cells must share one resolution. Neighbour functions work within a resolution; mixed-resolution sets need uncompacting first.
Smoothed surfaces are for display. Run spatial statistics on the raw counts.

Internal links

H3 explained: how a hexagonal grid indexes the whole Earth — pentagons, faces and why the grid folds
How to join two point datasets on an H3 index — choosing k for a distance, verified
How to aggregate points into H3 hexagons and map them — producing the counts being smoothed
How to fill a polygon with H3 cells — building the study area with its zeros
Fixing an H3 map that misleads — empty cells and unequal areas on the finished map
Kernel density explained — the continuous alternative to k-ring smoothing
How to calculate Moran's I in Python — testing dependence on the raw counts

FAQ

What is the difference between a disk and a ring?

grid_disk returns every cell within k steps, 1 + 3k(k + 1) cells including the origin. grid_ring returns only the cells exactly k steps away, 6k of them.

Why does H3 grid distance raise H3FailedError?

It works in local grid coordinates that cannot be unfolded across pentagons or over long distances. London to New York at res 5 raised, and so did 360 of 961 pairs near one pentagon. Fall back to great-circle distance.

Should I smooth with k=1 or k=2?

k=1 cut empty UK cells from 54.1% to 16.1% while keeping the peak at 15.14; k=2 reached 6.1% empty but flattened the peak to 10.95. Start with k=1 and move up only if the map is still speckled.

Can a k-ring disk replace a buffer?

Only as a candidate filter. A k=1 disk at res 8 reached 969 m one way and 1,402 m the other, so refine the candidates by true distance.

Is a smoothed H3 surface safe for statistics?

No. Smoothing made neighbouring UK cells correlate at 0.892 against 0.406 raw, entirely by construction. Use the raw counts with explicit neighbour weights.