An R-tree indexes things that have size: parcels, roads, building footprints, bounding boxes of detected objects, time ranges paired with price ranges. Point structures like the k-d tree split space at coordinates, which works for points but breaks down for an object that straddles the split line. The R-tree takes the opposite approach. It groups nearby objects, stores each group's minimum bounding rectangle (MBR), and lets those rectangles overlap. It is a balanced, page-oriented tree like a B-tree, which is why it became the standard spatial index in databases from PostGIS to SQLite.
This article builds an R-tree from first principles: the invariants, search, Guttman's insertion and split heuristics, the R* improvements, best-first nearest-neighbour search, deletion, and bulk loading with STR and Hilbert ordering. The Python code was run on a worked example of eight shops, and the outputs shown are the real ones. It closes with how the production implementations expose it, the failure modes that make an R-tree slow, and a checklist. If you index points only, compare it with the KD-tree article first; for one-dimensional intervals, an interval tree is simpler.
Structure and invariants
An R-tree of order (m, M) satisfies these invariants, from Antonin Guttman's 1984 SIGMOD paper:
- Every node holds between m and M entries, except the root, which holds at least two unless it is a leaf. M is chosen so a node fills a disk page or a few cache lines, and m is at most M/2.
- A leaf entry is (MBR, object id). An internal entry is (MBR, child pointer), where the MBR tightly bounds every rectangle in that child.
- All leaves are at the same depth, so a tree of N objects has height at most about log base m of N.
The key difference from a B-tree is that sibling MBRs may overlap, and nothing in the invariants limits by how much. Search correctness does not depend on overlap; search cost does, because a query point inside two overlapping MBRs must descend into both. Everything clever about R-tree variants is about keeping overlap and dead space (area inside an MBR that contains no object) small.
Structure and invariants
An R-tree of order (m, M) satisfies these invariants, from Antonin Guttman's 1984 SIGMOD paper:
- Every node holds between m and M entries, except the root, which holds at least two unless it is a leaf. M is chosen so a node fills a disk page or a few cache lines, and m is at most M/2.
- A leaf entry is (MBR, object id). An internal entry is (MBR, child pointer), where the MBR tightly bounds every rectangle in that child.
- All leaves are at the same depth, so a tree of N objects has height at most about log base m of N.
The key difference from a B-tree is that sibling MBRs may overlap, and nothing in the invariants limits by how much. Search correctness does not depend on overlap; search cost does, because a query point inside two overlapping MBRs must descend into both. Everything clever about R-tree variants is about keeping overlap and dead space (area inside an MBR that contains no object) small.
Window search: filter, then refine
A window query descends into every child whose MBR intersects the query rectangle, and reports leaf entries that intersect it. Rectangles are (xmin, ymin, xmax, ymax) tuples throughout.
def area(r): return (r[2]-r[0]) * (r[3]-r[1])
def union(r, s): return (min(r[0],s[0]), min(r[1],s[1]), max(r[2],s[2]), max(r[3],s[3]))
def intersects(r, s): return r[0] <= s[2] and s[0] <= r[2] and r[1] <= s[3] and s[1] <= r[3]
def enlargement(r, s): return area(union(r, s)) - area(r)
class Node:
def __init__(self, leaf):
self.leaf, self.entries = leaf, [] # entries: (rect, child_or_id)
def search(node, q):
out = []
for r, ch in node.entries:
if intersects(r, q):
out += [ch] if node.leaf else search(ch, q)
return outThe result is a candidate set, not a final answer: an MBR intersecting the query does not mean the polygon inside it does. Databases call this the filter step, followed by a refine step that runs the exact geometry test (orientation and point-in-polygon tests) on each candidate. In PostGIS, the && operator is the pure filter, and ST_Intersects uses the index for the filter step and then does the exact test.
Insertion and ChooseSubtree
Insertion descends from the root, choosing one child per level, appends the entry to a leaf, and splits nodes that overflow, propagating splits upward. If the root splits, a new root is created, which is the only way the tree grows taller, so all leaves stay at one depth.
M, m = 4, 2
class RTree:
def __init__(self):
self.root = Node(leaf=True)
def insert(self, rect, oid):
split = self._insert(self.root, rect, oid)
if split: # root overflowed: grow upward
a, b = split
self.root = Node(leaf=False)
self.root.entries = [(mbr(a.entries), a), (mbr(b.entries), b)]
def _insert(self, node, rect, oid):
if node.leaf:
node.entries.append((rect, oid))
else:
# ChooseSubtree: least enlargement, ties broken by smaller area
k = min(range(len(node.entries)),
key=lambda i: (enlargement(node.entries[i][0], rect), area(node.entries[i][0])))
child = node.entries[k][1]
split = self._insert(child, rect, oid)
if split:
a, b = split
node.entries[k] = (mbr(a.entries), a)
node.entries.append((mbr(b.entries), b))
else:
node.entries[k] = (union(node.entries[k][0], rect), child)
return quadratic_split(node) if len(node.entries) > M else NoneHere mbr(entries) is the union of the entries' rectangles. The ChooseSubtree rule (least area enlargement) is Guttman's. It is greedy, and it is the main reason a dynamically built R-tree degrades: early decisions are never revisited.
Splitting nodes: Guttman and R-star
When a node holds M+1 entries, they must be divided into two groups, each with at least m entries. Guttman proposed three strategies. The exhaustive split tries every partition and is exponential. The linear split picks two seeds that are far apart along some axis, then assigns the rest in arbitrary order. The quadratic split, below, picks as seeds the pair that would waste the most area if grouped together, then repeatedly assigns the entry with the strongest preference for one group.
def quadratic_split(node):
E = node.entries
_, i, j = max((area(union(E[i][0], E[j][0])) - area(E[i][0]) - area(E[j][0]), i, j)
for i in range(len(E)) for j in range(i+1, len(E)))
g1, g2 = [E[i]], [E[j]]
rest = [e for k, e in enumerate(E) if k not in (i, j)]
while rest:
if len(g1) + len(rest) == m: g1 += rest; break # guarantee minimum fill
if len(g2) + len(rest) == m: g2 += rest; break
r1, r2 = mbr(g1), mbr(g2)
e = max(rest, key=lambda e: abs(enlargement(r1, e[0]) - enlargement(r2, e[0])))
rest.remove(e)
d1, d2 = enlargement(r1, e[0]), enlargement(r2, e[0])
(g1 if (d1, area(r1), len(g1)) <= (d2, area(r2), len(g2)) else g2).append(e)
a, b = Node(node.leaf), Node(node.leaf)
a.entries, b.entries = g1, g2
return a, bThe R*-tree (Beckmann, Kriegel, Schneider and Seeger, 1990) changed three things and is what most serious implementations use. At the leaf level, ChooseSubtree minimises overlap enlargement rather than area enlargement. The split sorts entries along each axis, picks the axis with the smallest total perimeter over candidate distributions (perimeter favours square-ish rectangles), then the distribution with the least overlap. And on the first overflow at a given level during an insertion, instead of splitting, it removes the entries farthest from the node's centre (the paper recommends 30 percent) and reinserts them from the top. Forced reinsertion lets the tree undo early bad decisions, and the paper reports a minimum fill around 40 percent of M working well. The cost is slower inserts; the benefit is fewer node reads on every query after.
Worked example: eight shops
Insert eight shop footprints in the order A to H, with M = 4 and m = 2: A (1,1)-(2,2), B (2,6)-(3,7), C (6,1)-(7,2), D (7,7)-(8,8), E (1.5,2.5)-(2.5,3.5), F (6.5,6)-(7.5,6.8), G (3,3)-(4,4) and H (5,5)-(5.5,5.5). The first four fill the root leaf; inserting E overflows it. The quadratic split chooses A and D as seeds (the pair with the most wasted area), and the resulting tree, after all eight inserts, is:
root
(1, 1, 4, 7) leaf 1: A (1,1,2,2) E (1.5,2.5,2.5,3.5) B (2,6,3,7) G (3,3,4,4)
(5, 1, 8, 8) leaf 2: D (7,7,8,8) C (6,1,7,2) F (6.5,6,7.5,6.8) H (5,5,5.5,5.5)
search((0, 0, 3, 3)) -> ['A', 'E', 'G'] nodes visited: 2
nearest((6, 6), k=2) -> [('F', 0.5), ('H', 0.707)]The window query visits the root and leaf 1 only; leaf 2's MBR starts at x = 5, so it is pruned without being read. G is reported because its corner (3,3) touches the query's corner, which shows that intersects is closed: decide whether touching counts in your domain. The diagram shows the weakness too: B joined leaf 1, stretching its MBR to the top of the map, so any query in the empty upper-left region reads leaf 1 for nothing. That dead space is what R* reinsertion and bulk loading reduce.
Best-first nearest-neighbour search
Nearest-neighbour search uses a priority queue ordered by MINDIST, the distance from the query point to the nearest point of a rectangle, which is zero when the point is inside. This best-first strategy (Hjaltason and Samet, 1999) pops the closest item; if it is a node, it pushes the node's children; if it is an object, it is the next nearest result. Because MINDIST never overestimates the distance to anything inside a rectangle, objects come out in true distance order and the search reads the fewest possible nodes.
import heapq, itertools, math
def nearest(tree, p, k=1):
def mindist(r):
dx = max(r[0]-p[0], 0, p[0]-r[2])
dy = max(r[1]-p[1], 0, p[1]-r[3])
return dx*dx + dy*dy
tie = itertools.count() # avoids comparing Node objects
heap, out = [(0, next(tie), False, tree.root)], []
while heap and len(out) < k:
d, _, is_obj, item = heapq.heappop(heap)
if is_obj:
out.append((item, math.sqrt(d)))
continue
for r, ch in item.entries:
heapq.heappush(heap, (mindist(r), next(tie), item.leaf, ch))
return outIn the example, the point (6,6) lies inside leaf 2's MBR (distance 0) and 2 units from leaf 1's, so leaf 2 is expanded first; F (0.5 away) and H (0.707) come out before leaf 1 would be opened. As with window queries, the distance here is to the MBR; for real polygons, store the exact distance as the priority of object entries and the order stays correct.
Deletion and CondenseTree
Deletion finds the leaf holding the entry (possibly searching several paths, since MBRs overlap), removes it, and runs CondenseTree upward: nodes below m entries are detached and their entries set aside, path MBRs are tightened, and the orphans are reinserted at their original level. A root left with one child is replaced by it. Reinserting rather than merging gives orphans a better placement. An update is a delete plus an insert, so moving-object workloads often rebuild periodically instead.
Bulk loading: STR and Hilbert packing
When the data is known up front, building by repeated insertion wastes both time and quality. Bulk loading packs leaves full and nearly free of overlap:
- Sort-Tile-Recursive (STR) (Leutenegger, Lopez and Edgington, 1997): with N rectangles and capacity M, you need P = ceil(N/M) leaves. Sort by x-centre, cut into ceil(sqrt(P)) vertical slices, sort each slice by y-centre, and pack runs of M into leaves. Repeat on the leaf MBRs to build each level above. Shapely's
STRtreeis built this way. - Hilbert packing (Kamel and Faloutsos): sort objects by the Hilbert-curve value of their centres and pack consecutive runs. The Hilbert curve keeps nearby points close in the sort order better than a row-by-row or Z-order sort.
A packed tree is smaller (leaves are full rather than half to two-thirds full) and faster to query. The trade-off is that the structure is static: Shapely's STRtree cannot be modified after creation, so mixed workloads typically bulk-load a base and either rebuild periodically or keep a small dynamic tree for recent changes, the same base-plus-delta idea as an LSM tree.
R-trees in production systems
You will rarely ship the code above; you will configure one of these:
| System | How you use it | Notes |
|---|---|---|
| PostGIS | CREATE INDEX ON parcels USING GIST (geom); | R-tree-like index on GiST; && filters, <-> orders by distance for kNN |
| SQLite | CREATE VIRTUAL TABLE idx USING rtree(id, minX, maxX, minY, maxY); | R*Tree module; stores boxes only, join to the real table for refine |
| Shapely 2.x | STRtree(geoms).query(g, predicate="intersects") | Static STR tree; query returns integer indices into the input, not geometries |
| Python rtree | index.Index() | Bindings to libspatialindex; dynamic, supports R* and disk storage |
Index the geometry in the coordinate system you query in. An index on a geometry column does not help a predicate that transforms the column on the fly, since the planner cannot use it for the transformed expression. Check with EXPLAIN that the index is actually used.
Failure modes and tuning
| Symptom | Cause | Remedy |
|---|---|---|
| Queries read most of the tree | High overlap from incremental inserts | Bulk-load or rebuild; use R* splits |
| One huge object slows everything | A long road or country boundary inflates every ancestor MBR | Split large geometries into pieces before indexing |
| Index used, results still slow | Expensive refine step on complex polygons | Simplify for the filter, or pre-subdivide polygons |
| Index ignored by the planner | Predicate on a transformed column or a non-indexable operator | Store the projected geometry; check EXPLAIN |
| Wrong kNN order on polygons | Ranking by MBR distance only | Rank objects by exact distance in the refine step |
| Poor performance in many dimensions | MBRs overlap heavily beyond about ten dimensions | Use an approximate vector index instead |
Tuning order: measure node reads per query first (most libraries expose them, or count them as the worked example does). If reads are high, rebuild with bulk loading before touching node size; larger nodes help disk throughput but make each node's scan slower.
What to do next
To put an R-tree to work:
- Reproduce the worked example: insert the eight shops, print the tree, and confirm the query returns A, E and G in two node reads. Then insert them in a different order and compare.
- Add a node-read counter to search and measure it for your real query mix before and after STR bulk loading.
- In your database, create the spatial index, run
EXPLAINon your three most common queries, and fix any that do not use it. - Find your largest geometries by bounding-box area and test whether subdividing them reduces reads.
- Separate filter from refine in your code, and decide explicitly whether touching counts as intersecting.
- Compare against a point structure from the k-d tree deep dive if your objects are really points.