How to Aggregate Millions of Points into a Grid with DuckDB

Problem statement

You have millions of points and you need a summary: counts per cell, a mean per hexagon, a density surface, a table of totals per region.

GeoPandas can do it, and its ceiling arrives sooner than expected. Joining 13,472,195 real points to 4,596 provinces and counting per group measured:

                        time      peak memory
DuckDB                 50.9 s        290 MB
GeoPandas             121.8 s     23,813 MB

The memory is what bites: 23.8 GB for a job whose answer is a few thousand rows. (The DuckDB figures are Example 4's; the GeoPandas ones are from the benchmark guide, which runs the same join in both engines.) On an 8 GB laptop it is the difference between an answer and a crash, and on a shared machine it is the difference between finishing and being killed.

DuckDB streams. The aggregate is computed as the rows go past, so peak memory is a property of the output, not of the input.

Quick answer

Aggregate in SQL, and let the grid come from arithmetic rather than from a geometry join:

select floor(lon / 0.1) * 0.1 as cell_lon,
       floor(lat / 0.1) * 0.1 as cell_lat,
       count(*)               as n,
       avg(population)        as mean_population
from read_parquet('places.parquet')
where lat is not null
group by 1, 2
order by n desc;

Measured on 13,472,195 rows: 0.11 s from Parquet, against 0.62 s from the equivalent 1.79 GB text file. The busiest 0.1° cell held 3,610 points (Example 4).

No geometry, no join, no index — a grid cell is a rounded coordinate, and rounding is the fastest spatial operation there is.

13,472,195 GeoNames points; places.parquet 306 MB, allCountries.txt 1.79 GB

Geometry or arithmetic, cell size and CRS

1. Decide whether you need geometry at all

Three ways to aggregate points, in increasing order of cost:

If the output is a heat map or a density surface, the rectangular grid is almost always enough and it is two orders of magnitude cheaper.

2. Choose the cell size deliberately

The cell size is an analytical choice, not a rendering detail. Too large and the pattern disappears; too small and every cell holds one point.

select 0.1 as cell_degrees, count(*) as cells, max(n) as busiest, avg(n) as mean_n
from (select floor(lon / 0.1), floor(lat / 0.1), count(*) n
      from read_parquet('places.parquet') group by 1, 2);

Run it at three sizes and look at the distribution. A useful target is a mean of 10–50 points per occupied cell; much below that and the map shows sampling noise.

3. Work in a projected CRS if the cells must be equal-area

A grid in degrees is not a grid on the ground: a 0.1° cell is 11.1 km × 11.1 km at the equator and 11.1 km × 5.6 km at 60° north. For a national analysis that is usually fine; for anything comparing densities across latitudes it is a bias, not a rounding error.

The fix is to project once and grid in metres:

create table projected as
select st_transform(st_point(lon, lat), 'EPSG:4326', 'EPSG:3035',
                    always_xy := true) as geom, *
from read_parquet('places.parquet');

select floor(st_x(geom) / 10000) * 10000 as cell_x,
       floor(st_y(geom) / 10000) * 10000 as cell_y,
       count(*) as n
from projected group by 1, 2;

4. Aggregate into real areas with a spatial join when you need names

When the output has to be per district or per country, the grid does not help and a join is the operation:

select a.name, count(*) as n
from read_parquet('places.parquet') p
join st_read('provinces.shp') a
  on st_intersects(st_point(p.lon, p.lat), a.geom)
group by 1 order by n desc;

Measured on the full dataset: 13,472,195 points against 4,596 polygons in 50.9 s at 290 MB, matching 12,949,690 points in 249 countries. The 522,505 unmatched points — 3.9% — fall outside every Natural Earth province: offshore and small-island features (Example 4).

5. Keep the output small

The point of aggregating in the engine is that only the summary crosses into Python. A 13-million-row input producing 251 rows means 251 rows are materialised.

Resist select * before the aggregate. Every column named is a column read, and on a columnar file that is the difference between reading 0.06% of the file and 26% of it (measured here).

6. Turn the result into geometry only at the end

select cell_lon, cell_lat, n,
       st_makeenvelope(cell_lon, cell_lat,
                       cell_lon + 0.1, cell_lat + 0.1) as geom
from grid;

Building the 23,460 one-degree cell polygons is instant. Building 13 million point geometries to throw them away is not.

Gridding every GeoNames feature

The data used on this page

The points are GeoNames' full allCountries dump — every feature, about 13.5 million, rebuilt daily at download.geonames.org, so counts drift a little between runs — kept both as the 1.79 GB text file and as places.parquet (also linked as geonames.parquet) with name, latitude, longitude, feature class, country and population. The polygons are Natural Earth's 1:10m admin-1 states and provinces, 4,596 of them, as provinces.shp.

import os
import zipfile

import duckdb
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://download.geonames.org/export/dump/allCountries.zip", "allCountries.zip")   # about 420 MB
if not os.path.exists("allCountries.txt"):
    zipfile.ZipFile("allCountries.zip").extract("allCountries.txt")
if not os.path.exists("places.parquet"):
    duckdb.sql("""
        copy (select column01 as name, column04 as lat, column05 as lon, column06 as feature_class,
                     column08 as country, column14 as population
              from read_csv('allCountries.txt', delim = '\\t', header = false, quote = '',
                            all_varchar = false, sample_size = -1))
        to 'places.parquet' (format parquet)""")
if not os.path.exists("geonames.parquet"):
    os.link("places.parquet", "geonames.parquet")
fetch("https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_1_states_provinces.zip", "provinces.zip")
if not os.path.exists("provinces.shp"):
    with zipfile.ZipFile("provinces.zip") as z:
        for member in z.namelist():
            if member.startswith("ne_10m_admin_1_states_provinces."):
                open("provinces" + os.path.splitext(member)[1], "wb").write(z.read(member))
n = duckdb.sql("select count(*) from read_parquet('places.parquet')").fetchone()[0]
print(f"{n:,} GeoNames points; places.parquet {os.path.getsize('places.parquet') / 1e6:.0f} MB, "
      f"allCountries.txt {os.path.getsize('allCountries.txt') / 1e9:.2f} GB")

13,472,195 GeoNames points; places.parquet 306 MB, allCountries.txt 1.79 GB

Example 1 — a grid aggregation with the cell size as a parameter

import duckdb


def grid_counts(source, cell_size=0.1, lon="lon", lat="lat",
                extra_aggregates=None, con=None):
    """Count points per square cell. No geometry, no join, no index."""
    con = con or duckdb.connect()
    con.execute("set enable_progress_bar = false")

    extra = ""
    for alias, expression in (extra_aggregates or {}).items():
        extra += f", {expression} as {alias}"

    frame = con.execute(f"""
        select floor({lon} / {cell_size}) * {cell_size} as cell_lon,
               floor({lat} / {cell_size}) * {cell_size} as cell_lat,
               count(*) as n{extra}
        from {source}
        where {lat} is not null and {lon} is not null
        group by 1, 2
    """).df()

    print(f"{len(frame):,} occupied cells at {cell_size}°  "
          f"(mean {frame['n'].mean():.1f}, max {frame['n'].max():,} per cell)")
    return frame

grid = grid_counts("read_parquet('geonames.parquet')", cell_size=1.0,
                   extra_aggregates={"mean_pop": "avg(population)"})

23,460 occupied cells at 1.0°  (mean 574.3, max 33,850 per cell)

Example 2 — choosing the cell size from the data

def cell_size_sweep(con, source, sizes=(0.02, 0.05, 0.1, 0.25, 0.5, 1.0),
                    lon="lon", lat="lat"):
    """The distribution at each size, so the choice is informed."""
    print(f"{'size':>6} {'cells':>10} {'mean/cell':>10} {'median':>8} {'max':>10}")
    for size in sizes:
        row = con.execute(f"""
            with grid as (
              select floor({lon} / {size}) gx, floor({lat} / {size}) gy, count(*) n
              from {source} where {lat} is not null group by 1, 2)
            select count(*), avg(n), median(n), max(n) from grid
        """).fetchone()
        print(f"{size:6.2f} {row[0]:10,} {row[1]:10.1f} {row[2]:8.0f} {row[3]:10,}")

cell_size_sweep(duckdb.connect(), "read_parquet('geonames.parquet')")

  size      cells  mean/cell   median        max
  0.02  6,416,842        2.1        1        632
  0.05  2,513,819        5.4        2      1,699
  0.10    961,927       14.0        5      3,610
  0.25    224,457       60.0       14      6,722
  0.50     71,990      187.1       33     14,678
  1.00     23,460      574.3       64     33,850

A median of 1 — the 0.02° row — means the grid is finer than the data supports and the map will show sampling noise. A median in the tens is usually where a density surface becomes readable.

Example 3 — a streaming aggregate under a hard memory limit

def aggregate_within_memory(source, group_expression, memory_limit="500MB",
                            temp_dir="/tmp/duckdb_spill"):
    """Prove the aggregate streams: cap the memory and watch it still finish."""
    import duckdb
    import resource
    import time

    con = duckdb.connect()
    con.execute("set enable_progress_bar = false")
    con.execute(f"set memory_limit = '{memory_limit}'")
    con.execute(f"set temp_directory = '{temp_dir}'")

    started = time.perf_counter()
    frame = con.execute(f"""
        select {group_expression} as bucket, count(*) as n
        from {source} group by 1 order by n desc
    """).df()
    elapsed = time.perf_counter() - started
    peak = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024

    print(f"{len(frame):,} groups in {elapsed:.1f}s, peak RSS {peak:,.0f} MB "
          f"(limit was {memory_limit})")
    return frame


US_JOIN = """(select a.name as bucket_name from read_parquet('places.parquet') p
              join (select name, geom from st_read('provinces.shp') where adm0_a3 = 'USA') a
                on st_intersects(st_point(p.lon, p.lat), a.geom)
              where p.country = 'US')"""

Run in a fresh process each time, so the peak memory belongs to the query alone:

import subprocess
import sys

CHILD = """
import duckdb, resource, time
con = duckdb.connect()
con.execute("set enable_progress_bar = false")
con.execute("install spatial; load spatial")
LIMIT
t = time.perf_counter()
n = con.execute(\"\"\"select count(*) from (select a.name, count(*) from read_parquet('places.parquet') p
    join (select name, geom from st_read('provinces.shp') where adm0_a3 = 'USA') a
      on st_intersects(st_point(p.lon, p.lat), a.geom)
    where p.country = 'US' group by 1)\"\"\").fetchone()[0]
peak = [l for l in open('/proc/self/status') if l.startswith('VmHWM')][0].split()[1]
print(f"{n} groups in {time.perf_counter() - t:.2f} s, peak {int(peak) / 1024:.0f} MB")
"""
for label, limit in (("unconstrained", ""), ("memory_limit 200MB", "con.execute(\"set memory_limit = '200MB'\")")):
    out = subprocess.run([sys.executable, "-W", "ignore", "-c", CHILD.replace("LIMIT", limit)],
                         capture_output=True, text=True)
    print(f"{label:20} {out.stdout.strip() or out.stderr.strip().splitlines()[-1]}")

unconstrained        51 groups in 3.59 s, peak 217 MB
memory_limit 200MB   51 groups in 3.53 s, peak 217 MB

Measured on the 2.2-million-point US join: unconstrained, 3.59 s at 217 MB; with memory_limit set to 200 MB, 3.53 s at 217 MB — no slower, because the aggregate never needed the memory.

Example 4 — the figures at the top of this page

import time

con = duckdb.connect()
con.execute("set enable_progress_bar = false")
GRID = """select floor(lon / 0.1) * 0.1 as cell_lon, floor(lat / 0.1) * 0.1 as cell_lat,
                 count(*) as n, avg(population) as mean_population
          from {src} where lat is not null group by 1, 2 order by n desc"""
for label, src in (("Parquet", "read_parquet('places.parquet')"),
                   ("1.79 GB text file", "read_csv('allCountries.txt', delim = '\\t', header = false, quote = '', "
                    "columns = {'geonameid': 'BIGINT', 'name': 'VARCHAR', 'ascii': 'VARCHAR', 'alt': 'VARCHAR', "
                    "'lat': 'DOUBLE', 'lon': 'DOUBLE', 'fclass': 'VARCHAR', 'fcode': 'VARCHAR', 'country': 'VARCHAR', "
                    "'cc2': 'VARCHAR', 'a1': 'VARCHAR', 'a2': 'VARCHAR', 'a3': 'VARCHAR', 'a4': 'VARCHAR', "
                    "'population': 'BIGINT', 'elevation': 'VARCHAR', 'dem': 'VARCHAR', 'tz': 'VARCHAR', 'mod': 'VARCHAR'})")):
    t = time.perf_counter()
    top = con.execute(GRID.format(src=src)).fetchone()
    print(f"0.1° grid from {label:18} {time.perf_counter() - t:5.2f} s; busiest cell {top[2]:,} points")

JOIN = """select count(*) filter (where a.name is not null), count(distinct a.adm0_a3), count(*)
          from read_parquet('places.parquet') p
          left join st_read('provinces.shp') a on st_intersects(st_point(p.lon, p.lat), a.geom)"""
child = """
import duckdb, time
con = duckdb.connect(); con.execute("install spatial; load spatial"); con.execute("set enable_progress_bar = false")
t = time.perf_counter()
matched, groups, total = con.execute(JOIN_SQL).fetchone()
peak = int([l for l in open('/proc/self/status') if l.startswith('VmHWM')][0].split()[1]) / 1024
print(f"{total:,} points against 4,596 polygons in {time.perf_counter() - t:.1f} s at {peak:.0f} MB: "
      f"{matched:,} matches in {groups} groups, {total - matched:,} unmatched ({(total - matched) / total:.1%})")
""".replace("JOIN_SQL", repr(JOIN))
print(subprocess.run([sys.executable, "-W", "ignore", "-c", child], capture_output=True, text=True).stdout.strip())
print(f"load average {os.getloadavg()[0]:.1f}")

0.1° grid from Parquet             0.11 s; busiest cell 3,610 points
0.1° grid from 1.79 GB text file   0.62 s; busiest cell 3,610 points
13,472,195 points against 4,596 polygons in 50.9 s at 290 MB: 12,949,690 matches in 249 groups, 522,505 unmatched (3.9%)
load average 16.7

Rounding beats a geometry join

Why rounding beats a geometry join for a grid

A grid cell is defined by arithmetic: floor(x / size). Computing it costs one division and one floor per row, both of which run in the CPU's vector units on a whole batch of values at once.

A spatial join to a grid layer costs geometry construction, a spatial structure, candidate testing and exact predicate evaluation for every point. It produces the same answer for a rectangular grid, and it is two orders of magnitude more work.

Use the join when the cells are not rectangles — real administrative areas, hexagons from a published layer, irregular catchments.

Why streaming makes the memory constant

A GROUP BY needs to hold one entry per group, not one per row. DuckDB reads batches, updates the hash table of groups, and discards the batch.

That is why the 13.5-million-row join and aggregate peaked at 290 MB while GeoPandas needed 23,813 MB: GeoPandas materialised 13.5 million shapely geometries and a joined frame; DuckDB held 249 accumulator rows and a batch.

The consequence is practical. Memory use is bounded by the number of groups, so aggregating a billion rows into a thousand groups is no heavier than aggregating a million.

Why the file format decides the speed

The same aggregation measured 0.62 s from a 1.79 GB TSV and 0.11 s from a 306 MB Parquet file (Example 4). The engine is identical; the difference is that Parquet lets it read two columns instead of nineteen.

Converting once with COPY ... TO ... (FORMAT PARQUET) is a few seconds and makes every subsequent aggregation an order of magnitude faster. For a dataset queried more than twice, it is free.

Why a degree grid is not an equal-area grid

A degree of longitude is 111 km at the equator, 78 km at 45° and 56 km at 60°. A 0.1° grid therefore has cells whose ground area halves between the tropics and northern Europe.

Counts per cell then mean different things at different latitudes, and a density map computed that way shows a latitude gradient that is not in the data. For a single country it rarely matters; for anything continental, project first and grid in metres.

Negative coordinates, nulls and degree grids

• floor() on negative coordinates still bins correctly — it rounds towards negative infinity, which is what a grid wants.
• Null coordinates silently form a group; filter them explicitly.
• The busiest cell is worth looking at. Measured, one 1° cell held 33,850 points — a data-density artefact, not a real hotspot.
• A degree grid is not equal-area. Project for anything continental.
• group by 1, 2 refers to output positions and is idiomatic in DuckDB.
• Build the cell polygons after aggregating, not before.
• median() in DuckDB is exact, which is convenient and costs more than an approximation on huge groups.
• Convert to Parquet once if you will run more than a couple of these.

Internal links

• DuckDB spatial explained: a spatial database that is just a file — the engine
• Columnar or row storage: why DuckDB is fast on wide tables — where the 17× comes from
• How to run a spatial join in DuckDB — aggregating into real areas
• DuckDB's memory model: streaming, spilling and limits — why the memory stays flat
• How to create a fishnet grid in GeoPandas — the geometry version
• How to make a hexbin map of points — hexagonal binning
• The modifiable areal unit problem explained — why the cell size changes the answer
• Point density explained — what the counts mean
• How to use H3 in DuckDB — hexagonal cells instead of a degree grid

FAQ

How do I aggregate millions of points into a grid?

floor(x / size) * size as the group key. Measured on 13,472,195 points, a 0.1° grid took 0.11 s from a Parquet file with no geometry involved.

How much memory does it need?

Bounded by the number of groups, not the number of rows. The 13.5-million-point join and aggregate peaked at 290 MB against GeoPandas' 23,813 MB.

Should I use a spatial join instead?

Only when the cells are not rectangles — real administrative areas, published hexagons, irregular catchments. For a rectangular grid, rounding gives the same answer far more cheaply.

What cell size should I use?

Sweep several and look at the distribution. A median of one point per occupied cell means the grid is finer than the data supports; tens per cell usually reads well.

Does a grid in degrees distort the result?

Yes, away from the equator: a 0.1° cell is 11.1 km wide at the equator and 5.6 km at 60° north. Project and grid in metres for anything continental.

Why is it so much faster from Parquet?

Because a columnar file lets the engine read the two columns the query names. The same aggregation took 0.62 s from a 1.79 GB text file and 0.11 s from a 306 MB Parquet file.