R-tree
Also known as: R*-tree · STR · Sort-Tile-Recursive · spatial index · Guttman · best-first search
A spatial index for objects that have size (segments, buildings, plots): each node stores a rectangle that contains its whole subtree, and a query descends only where the rectangle touches the window. On the world’s coastline: a map window, three ways to build the tree, and the distance from each of 45 thousand towns to the sea over the Earth’s surface.
A map on the screen shows a piece of the world, and only what falls into it needs drawing. The world’s coastline at a scale of 1:50 million is 58,987 segments. You can check each one for each position of the window, but every time the map moves that is another 58,987 checks. An R-tree puts neighbouring segments into groups of 16, draws a rectangle around each group, puts the groups into groups again, and so on up to a single root. A 10° × 5° window over the straits between two seas touches 499 segments, and to find them the tree checks 847 rectangles: the rest is discarded in whole branches.
Thin rectangles are the nodes of the chosen level, light blue ones the nodes the query opened, red ones the segments found. Switch the way of building. The packed tree and the trees assembled by one-by-one insertion find the same things in the window but look completely different: with packing, leaves appear in the empty ocean as strips tens of degrees of latitude tall; with insertion along the coast, the leaves hug the line.
The idea
A segment, a building or a plot has size, and a single point does not describe it. But it can always be covered by the smallest rectangle with sides along the axes, from the minima and maxima of its coordinates. Checking whether two such rectangles intersect takes four comparisons.
An R-tree resembles a B-tree: all leaves are at the same depth, and a node has from m to M children. But instead of keys a node stores a rectangle that contains its whole subtree. If the window does not touch a node’s rectangle, it touches no object of the subtree, and there is no need to descend into the node.
That is all a window search needs: a walk from the root that enters only nodes touching the window. Unlike a k-d tree, an R-tree does not divide the plane: the rectangles of neighbouring nodes may overlap. An object then lies in exactly one leaf (it never has to be cut or duplicated), but a query sometimes descends into several branches at once, and the tree has no worst-case guarantees. The quality of the tree, how small and how little overlapping the node rectangles are, decides everything.
How to build it
One-by-one insertion (Antonin Guttman, 1984). A new object descends from the root into the child whose rectangle it would enlarge least, and on a tie into the smaller one. In a leaf the object is added; if there are now M + 1 children, the node splits in two, and the rectangles are recomputed upwards. The split is the main subtlety. Guttman’s quadratic split chooses two “seeds”, the pair that together would leave the most empty area, and then distributes the rest one at a time: first the one for which the choice of group matters most, into the group it enlarges less. This is how a tree is built when objects arrive one at a time and the tree must work all along.
Packing, if all the objects are known in advance. The STR method (Sort-Tile-Recursive): sort the objects by the x of their centre, cut them into √(n/M) strips with equal numbers of objects, sort each strip by y, and take M in a row per leaf. Then the same with the leaves, and so on up to the root. That is two sorts per level, and the nodes are full.
The coastline shows that the choice is not obvious. Packing cuts the world into strips by the number of segments, not by width. Where coasts are dense, a strip is narrow and the leaves small. In the ocean the strip between 80° and 93° E takes in both the coast at 65° S and the islands near the equator, and a leaf stretches over 72° of latitude. All the leaves of the packed tree together cover 31,406 square degrees; the tree built by inserting in the order of walking the coast covers 3767, since neighbouring vertices of the polyline lie close anyway. Insertion in random order is in between, at 20,680.
But leaf area is not yet speed. On a thousand random 2° × 1° windows on the coasts from 10° W to 40° E and from 50° to 72° N, the packed tree opens the fewest nodes, 3.7, since every node in it is full, but checks more rectangles: 59 against 43 for insertion along the coast and 40 for insertion in random order. The insertion trees’ nodes are 60–70% full; there are more of them, but they are tighter.
Code
The nearest coast
The second question for a spatial index is what lies nearest to a point. For each of the 45,227 towns with ten thousand people or more, let us find the nearest point of the coastline, over the Earth’s surface rather than on a flat map. There are 58,999 coast points and 45,227 towns, and brute force is two and a half billion distances.
Best-first search keeps a heap of nodes and points, each with its own bound, no greater than the true distance to anything inside. The closest by bound always comes out of the heap first. A node is opened, and its children go into the heap with their bounds. The first point taken out of the heap is the answer. The town farthest from the sea on this map is Karamay: 2422 kilometres to the coast, and for it the search measured 287 distances to nodes and points, and on average over all towns 132 instead of 58,999.
Let a node’s bound be no greater than the distance from the query to any point of its subtree. Then the point that best-first search takes out of the heap first is the nearest to the query.
Let the point p at distance d be taken out. The heap gives up the smallest bound, so everything left in it has a bound of at least d. Any point r not yet taken out lies in the subtree of some element of the heap: initially the heap holds the root, and an opened node is replaced by its children, so the subtrees of the heap’s elements cover everything not yet opened. The distance to r is no less than that element’s bound, and the bound is no less than d. So no point is nearer than p.
What remains is the bound. A node’s rectangle is given in longitudes and latitudes, while distance is measured along a great circle, and here you cannot take the nearest corner: the point of the rectangle nearest to a town may be in the middle of a side.
Let the town be at the point (λ, φ), and the rectangle span longitudes [λ1, λ2] and latitudes [φ1, φ2]. If λ is inside [λ1, λ2], the nearest point is on the same meridian, and the distance is |φ − φ*|, where φ* is φ clamped to [φ1, φ2]. Otherwise the nearest point lies on the side with the smaller difference of longitudes Δλ, and its latitude is either α, where tan α = sin φ / (cos φ · cos Δλ), if α falls in [φ1, φ2], or one of the ends.
By the spherical law of cosines the angle d between the town and the point (λ′, φ′) satisfies cos d = sin φ sin φ′ + cos φ cos φ′ cos(λ′ − λ). Smaller d means larger cosine. For fixed φ′ the factor cos φ cos φ′ is non-negative, so the cosine is larger the smaller |λ′ − λ| is (around the circle, from 0 to π). If the town’s longitude is inside the rectangle, we can take λ′ = λ, and then cos d = cos(φ′ − φ): nearest is φ′ clamped to [φ1, φ2]. If it is outside, then on the longitude range [λ1, λ2], which does not contain λ, the difference around the circle is smallest at one of the ends, and at every latitude it is the same, nearer end.
On this side Δλ is constant, and cos d = sin φ · sin φ′ + cos φ cos Δλ · cos φ′ = A cos(φ′ − α), where A cos α = cos φ cos Δλ and A sin α = sin φ. The latitudes of the side occupy a range no longer than π inside [−π/2, π/2], and on such a range A cos(φ′ − α) increases towards α and decreases away from it. So the maximum is at α if it is inside, and otherwise at one of the ends. This is not just a bound but the exact distance to the rectangle, and all the more it is no greater than the distance to any point inside.
The listing takes the 6283 towns with a hundred thousand people or more, to finish within a second, and checks the first answer against brute force over all coast points. The scene computes every town; the five farthest from the sea are towns of one area between 44° and 47° N, from 2422 to 2394 kilometres. Their nearest point is sometimes on a gulf at the mouth of a northern river, sometimes on the delta of a southern one: the distances to the two coasts are almost equal. Distance is measured to the vertices of the coastline, not to the segments between them. Half the segments of this map are shorter than 7.6 kilometres, so the line itself can be a few kilometres nearer than the nearest vertex.
What it costs
A window query costs the number of opened nodes times M, plus the number of objects found. If the window is small and the node rectangles hardly overlap, a node or two per level are opened: O(log n). The worst case is worse: nodes that all touch the window and contain nothing in it force the whole tree open. A guarantee, on the order of √n + k block reads per query, comes from the priority R-tree (2004), but in practice packing or an R*-tree is usually enough.
Scanning all segments grows linearly, the tree barely at all: from 922 to 58,987 segments the time for the windows grows fivefold in Go and even less in Python. On the whole map the tree is almost three hundred times faster than brute force in Python and over forty times in Go: a scan of simple comparisons is very cheap in Go. Building by packing is sorts, O(n log n); one-by-one insertion is more expensive: a descent through the tree and a quadratic split every few insertions.
Where it is easy to go wrong
- Strict inequalities. Rectangles that touch do intersect; a segment along a meridian has a rectangle of zero width, a point one of zero area. With strict comparisons they will never be found.
- Map distance instead of distance on the sphere. A degree of longitude is shorter near the pole than at the equator, and across the 180th meridian a coast at the other edge of the map is actually close. The bound to a rectangle must be computed in the same metric as the answer, or best-first search will not return the nearest.
- The nearest corner of a rectangle as the bound. It is not a lower bound: the nearest point may be in the middle of a side, and the node with the answer comes out of the heap too late.
- Insertion without regard to order. A tree built by insertion depends on the order the objects came in. If the objects can be ordered so that neighbours come in a row (along a line, along a Hilbert curve), the tree comes out noticeably better.
- Deletion. A node left with fewer than m children after a deletion is dissolved entirely in Guttman’s scheme, and its objects are reinserted. Simply leaving a half-empty node is possible, but the tree degrades.
Where it works
Databases. The spatial index in PostGIS is an R-tree built on the generalised search tree GiST; SQLite has the R*Tree module, and the spatial indexes in MySQL and Oracle are R-trees too. A query for “objects in this rectangle” first selects candidates by the rectangles in the index and only then checks the exact geometry.
Geometry libraries. GEOS and JTS, and through them Shapely and many GIS programs, have STRtree, an R-tree packed by the STR method. It is built once over all the geometries and asked: what intersects this, what is nearest to this.
Maps and games. Drawing what is visible, picking an object with a click, finding collisions among moving objects: everywhere the answer needed is “what is near this rectangle”. For moving objects the tree has to be rebuilt often, and there a grid or a tree that is easy to adjust is more common.
History
Antonin Guttman proposed the R-tree in 1984 as an index for computer-aided design systems and geographic data, and right away with three ways of splitting a node: exhaustive, quadratic and linear. In 1990 Norbert Beckmann, Hans-Peter Kriegel, Ralf Schneider and Bernhard Seeger described the R*-tree: when choosing a child and when splitting, it takes into account not only area but also overlap and perimeter, and an overflowing node first tries to hand some objects over for reinsertion. Ibrahim Kamel and Christos Faloutsos ordered objects along the Hilbert curve in 1993–1994, and Scott Leutenegger, Mario Lopez and Jeffrey Edgington proposed STR packing in 1997. Best-first nearest-neighbour search over such trees was described by Gísli Hjaltason and Hanan Samet in 1995.
Problems
Many windows over many objects: one query for each position of the map.
What falls into the windowData: Coastlines: Natural Earth, public domain Places: GeoNames, CC BY 4.0