Create a Hexagon Grid and Count Points in PyQGIS

Forty thousand crime reports, traffic collisions or tree records plotted as dots tell you almost nothing: the dense areas become a solid smear and every symbol overlaps its neighbours. Aggregating them into a grid of equal cells turns the smear into a readable density surface, and because every cell has the same area, the counts are directly comparable — something administrative boundaries of wildly different sizes can never offer. Hexagons are the usual choice because every neighbour is the same distance away and the grid has no strong horizontal or vertical lines to draw the eye.

This recipe belongs to Vector Data Manipulation in PyQGIS. It creates a hexagon grid over a study area, counts points per cell with and without weights, handles empty and partial cells, and styles the counts as a map.

From overlapping dots to comparable cellsLeft panel: a dense cloud of small orange dots where the centre is a solid mass and individual density differences are invisible. Right panel: the same area covered by a hexagon grid, each cell shaded by the number of points it contains, from pale for few points to dark green for many, revealing a hotspot and a secondary cluster.Same 40,000 points, now comparabledots: density is guessworkhexagons: counts per equal area

Prerequisites

  • QGIS 3.40 LTR or newer, or the QGIS 4 series.
  • A point layer and a study-area polygon, both in a projected CRS with metre units. A grid built in degrees has cells whose ground area shrinks towards the poles, which defeats the purpose.
  • A rough idea of the scale readers care about — a street, a neighbourhood, a district — because that sets the cell size.

Choose a cell size

There is no correct cell size, only one that suits the question and the data. Too small and most cells hold zero or one point, so the map shows noise; too large and genuine patterns are averaged away. A useful starting point is a size at which the median non-empty cell holds somewhere between five and fifty points.

Cell size changes the storyThree panels. With 100 metre hexagons, 82 percent of cells are empty and the rest hold one or two points, so the pattern is speckle. With 400 metre hexagons the median occupied cell holds 14 points and two clusters are visible. With 1,600 metre hexagons the whole town falls into six cells and the clusters disappear.Too fine is noise, too coarse is a blur100 m82%cells emptymedian occupied: 1speckle400 m14median occupied cell31% emptytwo clusters visible1,600 m6cells cover the townmedian occupied: 5,900pattern averaged away

import math
import statistics
import processing
from qgis.core import QgsProject

incidents = QgsProject.instance().mapLayersByName("incidents_2025")[0]
study = QgsProject.instance().mapLayersByName("borough_boundary")[0]
assert not incidents.crs().isGeographic()

def trial(size):
    grid = processing.run("native:creategrid", {
        "TYPE": 4, "EXTENT": study.extent(), "HSPACING": size, "VSPACING": size,
        "HOVERLAY": 0, "VOVERLAY": 0, "CRS": incidents.crs(), "OUTPUT": "TEMPORARY_OUTPUT",
    })["OUTPUT"]
    counted = processing.run("native:countpointsinpolygon", {
        "POLYGONS": grid, "POINTS": incidents, "FIELD": "n", "OUTPUT": "TEMPORARY_OUTPUT",
    })["OUTPUT"]
    counts = [f["n"] for f in counted.getFeatures()]
    occupied = [c for c in counts if c > 0]
    return (len(counts), 1 - len(occupied) / len(counts),
            statistics.median(occupied) if occupied else 0)

for size in (100, 200, 400, 800, 1600):
    cells, empty, median = trial(size)
    print(f"{size:>5} m: {cells:>6} cells, {empty:5.0%} empty, median occupied {median}")

Breakdown: TYPE 4 is a hexagon grid (0 point, 1 line, 2 rectangle, 3 diamond). For hexagons, HSPACING and VSPACING are the distances between cell centres; setting them equal gives regular hexagons. Running a handful of sizes on temporary outputs takes seconds on most datasets and turns the choice into an informed one rather than a guess. Write down the size you settle on and why — it is a methodological decision a reader is entitled to question.

Build the grid over the study area

A grid built on the study area's bounding box covers far more than the area itself: corners of sea, neighbouring districts, empty land that will show as zero and make the map look emptier than it is. Keep only cells that intersect the study area, then count.

SIZE = 400

grid = processing.run("native:creategrid", {
    "TYPE": 4, "EXTENT": study.extent(), "HSPACING": SIZE, "VSPACING": SIZE,
    "HOVERLAY": 0, "VOVERLAY": 0, "CRS": study.crs(), "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

in_area = processing.run("native:extractbylocation", {
    "INPUT": grid, "PREDICATE": [0], "INTERSECT": study, "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

hexes = processing.run("native:countpointsinpolygon", {
    "POLYGONS": in_area, "POINTS": incidents, "WEIGHT": "", "CLASSFIELD": "",
    "FIELD": "n", "OUTPUT": "/data/work/incidents_hex400.gpkg",
})["OUTPUT"]

Breakdown: Extracting cells that intersect the boundary — rather than clipping them — keeps every hexagon whole, so all cells have equal area and counts stay comparable; the boundary is drawn on top to show where the study area ends. Clipping instead would leave partial cells along the edge with less area and artificially low counts. The extraction is the select-by-location pattern with predicate 0, intersect. Grid cells carry left, top, right, bottom and id fields, which are handy for joining results from different years onto the same cells.

Keep whole cells at the edgeLeft: hexagons along a coastline kept whole, extending slightly past the boundary, each with the same area so counts compare fairly. Right: the same hexagons clipped to the coastline; edge cells shrink to a fraction of their area and show misleadingly low counts, which would need dividing by area to correct.Clipped cells undercount; whole cells compare fairlywhole cells + outlineevery cell the same areaclipped to coastlineedge cells look quieter than they are

Weighted counts and category counts

A plain count treats every point as one. Often a point stands for more — a collision with several casualties, a household with several residents — or you want counts split by type.

casualties = processing.run("native:countpointsinpolygon", {
    "POLYGONS": in_area, "POINTS": incidents, "WEIGHT": "casualties",
    "FIELD": "casualties", "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

types = processing.run("native:countpointsinpolygon", {
    "POLYGONS": in_area, "POINTS": incidents, "CLASSFIELD": "category",
    "FIELD": "n_categories", "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

by_type = processing.run("native:joinbylocationsummary", {
    "INPUT": in_area, "JOIN": incidents, "PREDICATE": [0],
    "JOIN_FIELDS": ["category"], "SUMMARIES": [0, 3],
    "DISCARD_NONMATCHING": False, "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]

Breakdown: WEIGHT sums a numeric field instead of counting points, so the output field holds total casualties per cell. CLASSFIELD does something people often do not expect: it counts the number of distinct values of that field in each cell, not the number of points per value — a measure of variety, useful and easily misread. For counts per category, native:joinbylocationsummary with count and unique-values summaries, or a filtered count per category, gives the breakdown. Note that a point exactly on the shared edge of two hexagons is counted in both; with real-world coordinates this is rare, but with points snapped to a grid that coincides with the hexagon grid it is not.

Normalise, class and style the result

Raw counts mostly map where people are. Dividing by a denominator — population, road length, number of properties — turns a density map into a rate map, and a sensible classification stops a single extreme cell from flattening everything else.

rates = processing.run("native:joinbylocationsummary", {
    "INPUT": hexes, "JOIN": QgsProject.instance().mapLayersByName("address_points")[0],
    "PREDICATE": [0], "SUMMARIES": [0], "DISCARD_NONMATCHING": False,
    "OUTPUT": "TEMPORARY_OUTPUT",
})["OUTPUT"]
rates = processing.run("native:fieldcalculator", {
    "INPUT": rates, "FIELD_NAME": "per_1000_addr", "FIELD_TYPE": 0, "FIELD_PRECISION": 2,
    "FORMULA": 'CASE WHEN "fid_count" >= 20 THEN 1000.0 * "n" / "fid_count" END',
    "OUTPUT": "/data/work/incidents_hex400_rates.gpkg",
})["OUTPUT"]

Breakdown: Counting address points per cell gives a denominator on exactly the same geometry as the incident counts, and the rate expression divides the two. Cells with fewer than twenty addresses get a null rate rather than a number, because a rate built on three households swings wildly on a single event and would dominate any colour scale; the threshold is a judgement worth stating on the map. The summary field name — fid_count here — depends on which field the count summary is taken from, so print the output's field names once before relying on it.

from qgis.core import (
    QgsVectorLayer, QgsGraduatedSymbolRenderer, QgsClassificationQuantile,
    QgsStyle, QgsRendererRange, QgsFillSymbol,
)

layer = QgsVectorLayer(hexes, "incidents per 400 m hexagon", "ogr")
layer.setSubsetString('"n" > 0')

renderer = QgsGraduatedSymbolRenderer("n")
renderer.setClassificationMethod(QgsClassificationQuantile())
renderer.updateClasses(layer, 6)
renderer.updateColorRamp(QgsStyle.defaultStyle().colorRamp("Greens"))
for i, r in enumerate(renderer.ranges()):
    symbol = r.symbol().clone()
    symbol.symbolLayer(0).setStrokeColor(symbol.color().darker(115))
    renderer.updateRangeSymbol(i, symbol)
layer.setRenderer(renderer)
layer.setOpacity(0.85)
QgsProject.instance().addMapLayer(layer)

Breakdown: Hiding zero cells with a subset string keeps the map focused on where incidents happened; show them in a neutral grey instead if the absence of events is itself the finding. Quantile classes put an equal number of cells in each colour, which suits skewed count data far better than equal intervals, where one busy town centre cell would push every other cell into the lightest class. Slightly darker outlines separate adjacent cells of the same class. For choosing breaks more deliberately — natural breaks, or fixed thresholds agreed with the people who use the map — see creating a choropleth map and classifying with natural breaks.

QGIS version compatibility

native:creategrid has been available since QGIS 3.0 and native:countpointsinpolygon since 3.10; earlier releases use qgis:countpointsinpolygon with the same parameters. native:joinbylocationsummary replaced its qgis: predecessor in 3.20. QgsClassificationQuantile and setClassificationMethod date from 3.10. Everything shown runs unchanged on 3.40, 3.44 and the QGIS 4 series.

Troubleshooting

  • Hexagons look squashed. HSPACING and VSPACING differ, or the layer is in a geographic CRS.
  • Every count is zero. Points and grid are in different CRSs; the grid CRS parameter did not match the points.
  • The grid covers the sea. It was built on the extent; extract cells intersecting the study area.
  • Totals exceed the number of points. Points on shared cell edges are counted twice, or a weight field was used.
  • The map is one colour with a single dark cell. Equal-interval classes on skewed data; use quantiles or natural breaks.

Conclusion

Work in a projected CRS, test a few cell sizes and pick the one where occupied cells hold a meaningful number of points, keep whole cells that intersect the study area, and count with weights or categories as the question needs. Then normalise where a denominator exists, classify with quantiles or natural breaks, and style the result so readers compare equal areas rather than administrative accidents.

Frequently Asked Questions

Why hexagons instead of squares? Every hexagon neighbour is equidistant, cells approximate circles more closely, and the grid lacks long straight lines that draw the eye. Squares are fine when the output must align with a raster.

Can I produce the same thing as a raster? Yes — a heatmap or a point-density raster, but those smooth across cell boundaries. Hexagons keep discrete, countable units that can be joined and tabulated.

How do I compare two years on the same grid? Build the grid once, save it, and count each year's points into copies of it; join the results on the grid cell id.

Is there a hexagon index like H3 in QGIS? Not built in. Plugins provide H3; for most maps a Processing grid in a local projected CRS is simpler and has equal-area cells by construction.