⚙️ How-To DuckDB & analytic SQL Performance & scale DuckDB for Spatial Data · step 8 of 20 Updated

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. Aggregating 13,464,017 real points measured:

                        time      peak memory
DuckDB                 61.8 s        305 MB
GeoPandas             117.5 s      4,739 MB

The memory is what bites: 4.7 GB for a job that produced 251 rows of output. 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,464,017 rows: 0.06 s from Parquet, against 1.01 s from the equivalent 1.79 GB text file. The busiest cell held 33,850 points.

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

Flow from coordinate columns through rounding and grouping to output cells.
Use a spatial join only when the cells are not rectangles.

Step-by-step solution

1. Decide whether you need geometry at all

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

Method How Cost
Rectangular grid floor(x / size) * size arithmetic only
Named areas spatial join to polygons a spatial join
Hexagons / H3 index function, or a join to a hex layer between the two

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,464,017 points against 4,596 polygons in 61.8 s at 305 MB, producing 12,942,217 matches across 251 groups. The 521,800 unmatched points — 3.9% — were offshore and small-island features.

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.

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 20,000 cell polygons is instant. Building 13 million point geometries to throw them away is not.

Grid of cell counts, medians and busiest cells at five grid sizes.
The cell count is also the memory: a hash table holds one entry per cell.

Code examples

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,422 occupied cells at 1.0°  (mean 574.8, 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,}")

Measured on 13,464,017 points:

  size      cells  mean/cell   median        max
  0.02  6,414,445        2.1        1        631
  0.10    961,781       14.0        5      3,610
  0.25    224,369       60.0       14      6,722
  0.50     71,921      187.2       33     14,678
  1.00     23,422      574.8       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

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

Explanation

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 aggregate peaked at 305 MB while GeoPandas needed 4,739 MB: GeoPandas materialised 13.5 million shapely geometries and a joined frame; DuckDB held 251 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 1.01 s from a 1.79 GB TSV and 0.06 s from a 446 MB Parquet file. 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.

Bar chart of grid cell ground width shrinking with latitude.
Fine for one country; a bias, not a rounding error, for anything continental.

Edge cases or notes

  • 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.

FAQ

How do I aggregate millions of points into a grid?

floor(x / size) * size as the group key. Measured on 13,464,017 points, that took 0.06 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 aggregate peaked at 305 MB against GeoPandas' 4,739 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 1.01 s from a 1.79 GB text file and 0.06 s from a 446 MB Parquet file.