How to Assign Points to H3 Cells in Python

On this page

Problem statement

Assigning a point to an H3 cell is one function call: h3.latlng_to_cell(lat, lng, res). Doing it for a whole table goes wrong in three ways that the one-line examples do not show.

It can be slow for the wrong reason. Measured on one million GeoNames coordinates at resolution 9, DataFrame.apply took 5.81 s. A plain list comprehension doing identical work took 0.51 s, the integer API 0.38 s, and DuckDB's H3 extension 0.42 s including the file read (Example 4). The H3 call itself is cheap; the time goes on the per-row machinery around it.

It can store the result expensively. Kept as the usual hexadecimal strings, two million cells occupy 46.0 MB in pandas. The same cells as unsigned 64-bit integers take 16.0 MB.

It can be silently wrong. Swap latitude and longitude and H3 raises nothing. It returns a valid cell somewhere else on Earth. Over 100,000 GeoNames points with swapped columns, the median cell was 7,007 km from the right one. Only non-finite input such as NaN raises an error. This guide covers a fast, integer-based, validated assignment for tables of any size.

Quick answer

import h3
import pandas as pd

df = pd.read_parquet("cities500.parquet", columns=["name", "lat", "lon"])
df["h3_9"] = [h3.latlng_to_cell(lat, lng, 9)
              for lat, lng in zip(df["lat"].tolist(), df["lon"].tolist())]
print(df.head(3))
      name       lat      lon             h3_9
0     Vila  42.53176  1.56654  89396223157ffff
1   Soldeu  42.57688  1.66769  893962228bbffff
2  Sispony  42.53368  1.51613  89396223687ffff

That is every row of cities500 with h3 4.5. Latitude goes first. For large tables, use h3.api.basic_int to get integers directly, and check the columns before trusting the output: nothing else will warn you that they are swapped.

Time to assign one million GeoNames points to H3 resolution 9 by five methods. 1,000,000 points β†’ resolution 9 cellsDataFrame.apply5.81 snp.vectorize0.57 slist comprehension0.51 sbasic_int API0.38 sDuckDB h3 (incl. read)0.42 s07 sapply builds a pandas Series for every row; the H3 call itself is the cheap part spatialworkflow.io
Measured with h3 4.5 on GeoNames coordinates; all four Python routes returned identical cells.

Check degrees and order, pick a resolution, store integers

1. Check that the columns are degrees, in the right order

H3 wants WGS84 latitude and longitude in degrees, latitude first. Before indexing anything, look at the ranges:

def coordinate_report(df, lat="lat", lon="lon"):
    la, lo = df[lat], df[lon]
    print(f"{lat}: {la.min():.4f} .. {la.max():.4f}   {lon}: {lo.min():.4f} .. {lo.max():.4f}")
    print(f"missing: {la.isna().sum() + lo.isna().sum():,}")
    print(f"|{lat}| > 90 (impossible): {(la.abs() > 90).sum():,}")
    print(f"|{lon}| > 90 (impossible if swapped): {(lo.abs() > 90).sum():,}")


geonames = pd.read_parquet("geonames_all.parquet")
coordinate_report(geonames)
print()
coordinate_report(geonames.rename(columns={"lat": "lon", "lon": "lat"}))      # deliberately swapped
lat: -90.0000 .. 90.0000   lon: -180.0000 .. 180.0000
missing: 0
|lat| > 90 (impossible): 0
|lon| > 90 (impossible if swapped): 4,659,217

lat: -180.0000 .. 180.0000   lon: -90.0000 .. 90.0000
missing: 0
|lat| > 90 (impossible): 4,659,217
|lon| > 90 (impossible if swapped): 0

That is the whole GeoNames dump, with the columns right and then deliberately swapped.

A "latitude" column with values beyond Β±90 is proof of a swap. A dataset confined to longitudes between βˆ’90 and 90, such as most of Europe and Africa, gives no such proof, so there you must check a known location by eye. Values in the hundreds of thousands mean projected metres, which have to be transformed to EPSG:4326 first.

2. Choose the resolution for the question

The resolution sets how far each point is moved to its cell centre and how many points share a cell. At resolution 9, the 235,735 cities500 places were a median 128 m and at most 219 m from their cell centres. Choosing an H3 resolution covers how to choose. Put the level in the column name (h3_9) so tables at different levels are never joined by mistake.

3. Loop with a comprehension, never with apply

The same one million coordinates, timed five ways:

import h3
import h3.api.basic_int as hi
import numpy as np

la, lo = df["lat"].to_numpy(), df["lon"].to_numpy()   # df: 1,000,000 GeoNames rows; timings from Example 4

a = [h3.latlng_to_cell(y, x, 9) for y, x in zip(la, lo)]                      # 0.57 s
b = [h3.latlng_to_cell(y, x, 9) for y, x in zip(la.tolist(), lo.tolist())]    # 0.53 s
c = [hi.latlng_to_cell(y, x, 9) for y, x in zip(la.tolist(), lo.tolist())]    # 0.41 s
d = np.vectorize(h3.latlng_to_cell, otypes=[object])(la, lo, 9)                # 0.57 s
e = df.apply(lambda r: h3.latlng_to_cell(r.lat, r.lon, 9), axis=1)             # 6.11 s

All five returned identical cells. np.vectorize is not vectorised: it is a Python loop with extra setup. apply(axis=1) builds a pandas Series for every row, and that is where 12 times the time goes. The integer API is fastest because it skips formatting each 64-bit index as a string.

4. Store cells as integers

import h3.api.basic_int as hi

df["h3_9"] = pd.Series([hi.latlng_to_cell(lat, lng, 9)
                        for lat, lng in zip(df["lat"].tolist(), df["lon"].tolist())],
                       dtype="uint64", index=df.index)

The string and integer forms are the same 64-bit number. h3.str_to_int("89194ad30d3ffff") is 617438095495921663, and h3.int_to_str converts back. Measured on cities500, the string column took 5.43 MB in pandas 3 (Arrow-backed str) and the integer column 1.89 MB. The largest cell identifier in GeoNames, 621,280,605,650,812,927, is below 2⁢³, so a signed int64 column works too, which matters for engines without unsigned types.

5. Deal with missing and impossible coordinates explicitly

NaN and infinity raise. Everything else that is finite is accepted:

import h3

try:
    h3.latlng_to_cell(float("nan"), -0.1278, 9)
except h3.H3LatLngDomainError as exc:
    print(type(exc).__name__, repr(str(exc)))

print(h3.latlng_to_cell(-0.1278, 51.5074, 9))   # swapped: no error
H3LatLngDomainError ''
897b81723b3ffff

Uncaught, the traceback ends in a bare h3._cy.error_system.H3LatLngDomainError, because the error carries no message. It derives from ValueError (H3LatLngDomainError β†’ H3ValueError β†’ H3BaseException β†’ ValueError), so one NaN in ten million rows stops a plain loop at that row. Mask unusable rows before the loop rather than catching exceptions inside it. Example 1 does this.

6. Parallelise, or move the work into DuckDB

A single process indexed all 13,472,318 GeoNames points in 5.1 s. Four worker processes took 1.8 s (Example 2). DuckDB's h3 community extension adds the column while reading Parquet, and it wrote the whole table with an h3_9 column in 1.4 s on four threads (Example 3). A 300,000-row sample of DuckDB's cells matched the Python library's exactly.

7. Verify on a sample

Measure the distance from each point to its cell centre. It should never exceed about one edge length at your resolution:

import numpy as np

import h3.api.basic_int as h3i

cc = np.array([h3i.cell_to_latlng(int(c)) for c in df["h3_9"]])      # the column holds integers
R = 6371008.8
p1, p2 = np.radians(df["lat"].to_numpy()), np.radians(cc[:, 0])
dl = np.radians(cc[:, 1] - df["lon"].to_numpy())
d = 2 * R * np.arcsin(np.sqrt(np.sin((p2 - p1) / 2) ** 2
                              + np.cos(p1) * np.cos(p2) * np.sin(dl / 2) ** 2))
print(f"distance to cell centre: median {np.median(d):.0f} m, max {d.max():.0f} m")
distance to cell centre: median 128 m, max 219 m

This check catches arithmetic mistakes, not a consistent swap: swapped coordinates are close to their own, wrong, cells. For a swap, map a few known places. Fixing H3 cells in the wrong place covers that diagnosis.

The same H3 cells stored as hexadecimal strings and as unsigned 64-bit integers. one identifier, two encodingsmeasurestringuint64example89194ad30d3ffff617438095495921663pandas, 2M cells46.0 MB16.0 MBParquet zstd, 13.47M57.7 MB41.4 MBlargest valueβ€”below 2^63: fits int64h3.str_to_int and h3.int_to_str convert losslessly either wayCompression narrows the gap on disk; in memory the integers are almost three times smaller spatialworkflow.io
Integers also join and group faster, and every engine has a 64-bit integer type.

A validated, parallel assignment of GeoNames points

The data used on this page

Every feature in the GeoNames allCountries dump as a point; GeoNames' cities500, every populated place of 500 people or more. 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)


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], names=["lat", "lon"],
                    quoting=3).to_parquet("geonames_all.parquet")
if not os.path.exists("cities500.parquet"):
    # GeoNames cities500 as saved on 8 October 2026 (rebuilt nightly); live: https://download.geonames.org/export/dump/cities500.zip
    fetch("https://data.spatialworkflow.io/geonames-cities500-2026-10-08.zip", "cities500.zip")
    with zipfile.ZipFile("cities500.zip") as z, z.open("cities500.txt") as f:
        pd.read_csv(f, sep="\t", header=None, usecols=[1, 4, 5, 8, 14], quoting=3,
                    names=["name", "lat", "lon", "country", "population"],
                    keep_default_na=False, na_values={"population": [""]}).to_parquet("cities500.parquet")
print(len(pd.read_parquet("geonames_all.parquet", columns=["lat"])), len(pd.read_parquet("cities500.parquet", columns=["name"])))
13472318 236017

Example 1 β€” a validated assignment that returns nullable integers

import h3
import h3.api.basic_int as h3i
import numpy as np
import pandas as pd


def assign_h3(df, res, lat="lat", lon="lon", as_int=True):
    """Return a Series of H3 cells; rows with unusable coordinates get <NA>."""
    la = df[lat].to_numpy(dtype="float64", na_value=np.nan)
    lo = df[lon].to_numpy(dtype="float64", na_value=np.nan)
    impossible = np.abs(la) > 90
    if impossible.mean() > 0.01:
        raise ValueError(f"{impossible.sum():,} rows have |{lat}| > 90 β€” are {lat} and {lon} swapped?")
    usable = np.isfinite(la) & np.isfinite(lo) & ~impossible & (np.abs(lo) <= 180)

    fn = h3i.latlng_to_cell if as_int else h3.latlng_to_cell
    values = [fn(a, b, res) for a, b in zip(la[usable].tolist(), lo[usable].tolist())]

    out = pd.Series(pd.NA, index=df.index, dtype="UInt64" if as_int else "string",
                    name=f"h3_{res}")
    out[usable] = values
    skipped = int((~usable).sum())
    if skipped:
        print(f"assign_h3: {skipped:,} of {len(df):,} rows had no usable coordinates")
    return out


import time

cities = pd.read_parquet("cities500.parquet", columns=["name", "lat", "lon"])
test = pd.concat([cities, pd.DataFrame({"name": ["no coordinates", "impossible"],
                                        "lat": [np.nan, 91.5], "lon": [0.0, 0.0]})], ignore_index=True)
start = time.perf_counter()
print(assign_h3(test, 9).tail(3))
print(f"{time.perf_counter() - start:.2f} s")
try:
    assign_h3(test.rename(columns={"lat": "lon", "lon": "lat"}), 9)
except ValueError as exc:
    print("ValueError:", exc)
assign_h3: 2 of 236,019 rows had no usable coordinates
236016    619634483258195967
236017                  <NA>
236018                  <NA>
Name: h3_9, dtype: UInt64
0.09 s
ValueError: 71,286 rows have |lat| > 90 β€” are lat and lon swapped?

That is cities500 with a NaN row and a latitude of 91.5 appended, then with the columns swapped.

It took 0.09 s. With the columns swapped, the 1% threshold stops it before any cell is computed: ValueError: 71,286 rows have |lat| > 90 β€” are lat and lon swapped? A handful of bad rows passes through as missing values, while a systematic problem stops the run.

Example 2 β€” the same assignment across worker processes

from multiprocessing import get_context

import h3.api.basic_int as h3i
import numpy as np


def _cells(args):
    lat, lng, res = args
    return np.fromiter((h3i.latlng_to_cell(a, b, res) for a, b in zip(lat.tolist(), lng.tolist())),
                       dtype=np.uint64, count=len(lat))


def cells_parallel(lat, lng, res, processes=4, chunk=500_000):
    parts = [(lat[i:i + chunk], lng[i:i + chunk], res) for i in range(0, len(lat), chunk)]
    with get_context("fork").Pool(processes) as pool:
        return np.concatenate(pool.map(_cells, parts))


import time

import pandas as pd

points = pd.read_parquet("geonames_all.parquet").dropna()
lat, lng = points["lat"].to_numpy(), points["lon"].to_numpy()
for run in range(2):
    start = time.perf_counter()
    single = _cells((lat, lng, 9))
    one = time.perf_counter() - start
    start = time.perf_counter()
    many = cells_parallel(lat, lng, 9)
    four = time.perf_counter() - start
    print(f"run {run}: {len(lat):,} points  single {one:.1f}s  Pool(4) {four:.1f}s  "
          f"identical {np.array_equal(single, many)}")
run 0: 13,472,318 points  single 5.1s  Pool(4) 1.8s  identical True
run 1: 13,472,318 points  single 5.1s  Pool(4) 1.7s  identical True

Four processes gave about a 3Γ— speed-up, not 4Γ—. Each chunk's arrays and results are pickled across process boundaries, and that cost does not shrink as workers are added. The fork context works on Linux. On macOS and Windows, where spawn is the default, _cells must live in an importable module; see the multiprocessing fix.

Example 3 β€” add the cell column while rewriting a Parquet file in DuckDB

import os
import time

import duckdb


def cells_duckdb(src, dst, res=9, threads=4):
    con = duckdb.connect()
    con.execute("install h3 from community")
    con.execute("load h3")
    con.execute(f"set threads = {threads}")
    con.execute("set enable_progress_bar = false")
    start = time.perf_counter()
    con.execute(f"""
        copy (select *, h3_latlng_to_cell(lat, lon, {res}) as h3_{res}
              from read_parquet('{src}')
              where lat is not null and lon is not null)
        to '{dst}' (format parquet, compression zstd)
    """)
    rows, kind = con.execute(f"select count(*), any_value(typeof(h3_{res})) from read_parquet('{dst}')").fetchone()
    print(f"{rows:,} rows, h3_{res} as {kind}, {time.perf_counter() - start:.1f}s, "
          f"{os.path.getsize(dst) / 1e6:,.0f} MB")


cells_duckdb("geonames_all.parquet", "geonames_h3.parquet")
13,472,318 rows, h3_9 as UBIGINT, 1.4s, 149 MB

The column comes out as UBIGINT, DuckDB's unsigned 64-bit integer, which is the same number Python's integer API returns. None of the table passes through Python objects. How to use H3 in DuckDB continues from here.

Example 4 β€” the measurements quoted on this page

import os
import time

import duckdb
import h3
import h3.api.basic_int as h3i
import numpy as np
import pandas as pd
from pyproj import Geod

points = pd.read_parquet("geonames_all.parquet").dropna()
million = points.sample(1_000_000, random_state=0)
lats, lons = million["lat"].tolist(), million["lon"].tolist()
la, lo = million["lat"].to_numpy(), million["lon"].to_numpy()
million.to_parquet("million.parquet")
con = duckdb.connect()
con.execute("install h3 from community; load h3")
for label, run in [
        ("DataFrame.apply", lambda: million.apply(lambda r: h3.latlng_to_cell(r.lat, r.lon, 9), axis=1)),
        ("np.vectorize", lambda: np.vectorize(h3.latlng_to_cell, otypes=[object])(la, lo, 9)),
        ("over NumPy scalars", lambda: [h3.latlng_to_cell(a, b, 9) for a, b in zip(la, lo)]),
        ("list comprehension", lambda: [h3.latlng_to_cell(a, b, 9) for a, b in zip(lats, lons)]),
        ("integer API", lambda: [h3i.latlng_to_cell(a, b, 9) for a, b in zip(lats, lons)]),
        ("DuckDB, with the read", lambda: con.execute("select h3_latlng_to_cell(lat, lon, 9) from 'million.parquet'").fetchall())]:
    start = time.perf_counter()
    run()
    print(f"{label:22} {time.perf_counter() - start:5.2f} s per million points")

geod = Geod(ellps="WGS84")
sample = points.sample(100_000, random_state=1)
right = np.array([h3.cell_to_latlng(h3.latlng_to_cell(a, b, 9)) for a, b in zip(sample.lat, sample.lon)])
wrong = np.array([h3.cell_to_latlng(h3.latlng_to_cell(b, a, 9)) for a, b in zip(sample.lat, sample.lon)])
moved = np.array(geod.inv(right[:, 1], right[:, 0], wrong[:, 1], wrong[:, 0])[2]) / 1000
print(f"swapped columns move the median cell {np.median(moved):,.0f} km; "
      f"London {geod.inv(-0.1278, 51.5074, 51.5074, -0.1278)[2] / 1000:,.0f} km")

two = points.head(2_000_000)
as_str = pd.Series([h3.latlng_to_cell(a, b, 9) for a, b in zip(two.lat, two.lon)], dtype="str")
as_int = pd.Series([h3i.latlng_to_cell(a, b, 9) for a, b in zip(two.lat, two.lon)], dtype="uint64")
print(f"two million cells in memory: str {as_str.memory_usage(deep=True) / 1e6:.1f} MB, "
      f"uint64 {as_int.memory_usage(deep=True) / 1e6:.1f} MB")

for name, expr in [("int", "h3_latlng_to_cell(lat, lon, 9)"), ("str", "h3_h3_to_string(h3_latlng_to_cell(lat, lon, 9))")]:
    con.execute(f"copy (select {expr} as h3_9 from 'geonames_all.parquet' where lat is not null) "
                f"to 'cells_{name}.parquet' (format parquet, compression zstd)")
    print(f"every cell as {name}: {os.path.getsize(f'cells_{name}.parquet') / 1e6:.1f} MB of Parquet")
busiest = con.execute("""select count(*) from 'geonames_all.parquet' where lat is not null
                         group by h3_latlng_to_cell(lat, lon, 10) order by 1 desc limit 1""").fetchone()[0]
print(f"busiest resolution-10 cell: {busiest} features in about {h3.average_hexagon_area(10, 'km^2'):.3f} kmΒ²")
DataFrame.apply         5.81 s per million points
np.vectorize            0.57 s per million points
over NumPy scalars      0.57 s per million points
list comprehension      0.51 s per million points
integer API             0.38 s per million points
DuckDB, with the read   0.42 s per million points
swapped columns move the median cell 7,007 km; London 7,484 km
two million cells in memory: str 46.0 MB, uint64 16.0 MB
every cell as int: 41.4 MB of Parquet
every cell as str: 57.7 MB of Parquet
busiest resolution-10 cell: 360 features in about 0.015 kmΒ²

Comprehensions over apply, and swaps that never raise

Why apply is twelve times slower than a comprehension

latlng_to_cell is a thin wrapper around a C function, and at 0.38 s per million, the integer API costs about 380 ns per point including the Python loop. DataFrame.apply(axis=1) adds a pandas Series for every row: an index, a dtype check and attribute lookups for r.lat and r.lon. That wrapping costs about ten times more than the geometry.

.tolist() helps a little because iterating a Python list of floats avoids creating a NumPy scalar per element. The larger gain is dropping string formatting with the integer API, which brings it to 0.38 s.

Why swapped coordinates do not raise

H3 converts the two numbers to radians and then to a point on the unit sphere with sines and cosines, and every finite pair of angles is some point on a sphere. Latitude 95 is simply 5Β° past the pole. Measured, latlng_to_cell(95.0, 10.0, 9) returned a cell centred at 84.999Β° N, 169.999Β° W, which is the same point reached by going over the pole to the opposite meridian.

Swapped coordinates are therefore different positions, not invalid ones. London's coordinates reversed land in the Indian Ocean, east of the Somali coast, 7,484 km away. Across 100,000 GeoNames points the median displacement was 7,007 km. H3 cannot detect this; only the data's ranges can, and only when they happen to exceed Β±90.

Why integers are the right storage

An H3 index is a 64-bit integer: one reserved high bit, 4 bits of mode, 3 more reserved bits, 4 bits of resolution, 7 bits of base cell and 3 bits for each of 15 resolution digits. The hexadecimal string is a display format of that integer. Storing strings costs 15 characters plus offsets per row, and every join and group-by then compares strings instead of machine words.

Measured, in memory the gap is nearly threefold: 46.0 MB against 16.0 MB for two million cells. On disk, Parquet's dictionary encoding and zstd narrow it to 57.7 MB against 41.4 MB for all 13.47 million cells. Keep strings for display, logs and URLs, and integers everywhere else.

Why parallelism gives less than it promises

What a process pool actually does per chunk when indexing points into H3 cells. why four processes do not give four timeschunk500,000 pointspickle outto a workerindexC, no GILpickle backthe result arrayconcatenatein the parent5.9 s became 2.3 s on four processes; the rest is serialisation and the parent’s own work spatialworkflow.io
DuckDB skips the two orange steps entirely by running the same C library over columnar batches.

Indexing is embarrassingly parallel, since no point depends on another, but the process pool is not free. Each 500,000-point chunk is pickled to a worker, and each result array is pickled back and concatenated. On 13.47 million points, four processes turned 5.1 s into 1.8 s. The remaining time is serialisation and the parent's own work, and more processes only add to it.

DuckDB avoids the problem by running the same C library on its own threads over columnar batches, with no Python objects in between. When the data is already in Parquet, that is the better route.

What h3.latlng_to_cell does with five kinds of bad coordinate input. only some bad input raisesinputwhat H3 doeslongitude and latitude swappedvalid cell, median 7,007 km away β€”silentlatitude 95valid cell near 85Β° N β€” silentprojected metresvalid cell somewhere β€” silentNaN or infinityH3LatLngDomainErrorNoneTypeError: must be real numberThe three silent cases are the dangerous ones: check the columns before indexing spatialworkflow.io
H3 treats any finite pair of numbers as a position on the sphere, so it cannot tell you the order was wrong.

Latitude first, projected input and domain errors

  • Latitude first. latlng_to_cell(lat, lng, res) is the reverse of Shapely's (x, y) and of GeoJSON. A swapped call raises nothing.
  • Projected coordinates are accepted too. latlng_to_cell(180000.0, 530000.0, 9) returned a valid cell. Transform to EPSG:4326 first.
  • NaN and infinity raise H3LatLngDomainError with an empty message; None raises TypeError: must be real number, not NoneType. Mask them before the loop.
  • Resolution 16 raises H3ResDomainError. Valid levels are 0 to 15.
  • The v4 API renamed everything. geo_to_h3 is gone and latlng_to_cell replaces it; see the v3 to v4 rename.
  • Nullable UInt64 keeps missing rows as <NA>, while plain uint64 cannot hold them. Use the nullable type if any row can be unusable.
  • Duplicate coordinates are common in real data. The busiest resolution-10 cell in GeoNames, about 0.015 kmΒ², holds 360 features. That is a feature of the source, not a bug in the index.

FAQ

How do I convert latitude and longitude to an H3 cell in Python?

Call h3.latlng_to_cell(lat, lng, res) with latitude first, inside a list comprehension over the two columns. With h3 4.5 that indexed 235,735 points in 0.23 s.

Why is DataFrame.apply so slow for H3?

It constructs a pandas Series for every row before H3 is even called. On one million points it took 5.81 s, against 0.51 s for a comprehension doing the same work.

Should I store H3 cells as strings or integers?

Integers. They are the same 64-bit value. Two million cells took 16.0 MB as integers against 46.0 MB as strings, and integer joins and group-bys compare machine words rather than text.

Does H3 raise an error if latitude and longitude are swapped?

No. Any finite pair of numbers is a valid position on the sphere, so you get a real cell in the wrong place, a median of 7,007 km away in a GeoNames test. Only NaN and infinity raise.

How do I index tens of millions of points quickly?

Use DuckDB's h3 extension on the Parquet file, or split the arrays across worker processes. All 13.47 million GeoNames points took 1.4 s to rewrite with a cell column in DuckDB and 1.8 s with four Python processes.

What happens to rows with missing coordinates?

latlng_to_cell raises H3LatLngDomainError on NaN, which stops a loop. Mask those rows first and store the result in a nullable integer column so they stay as missing values.