Spatial Indexing

When you query a traditional (relational) database column, a B-tree index usually handles the search efficiently.
Spatial data, however, is not just a single dimension; it’s at least two (latitude, longitude) and often more (3D coordinates, time, attributes). This makes spatial indexing a multidimensional problem. Without spatial indexing, every query would require scanning the entire dataset something that quickly becomes unmanageable with millions of points.

It’s not something most of us will ever encounter, but with all things involving data it’s good to have an understanding of what possible solutions exist and how you could potentially leverage them.

Why Spatial Indexing is Different

Spatial indexes must handle:

• Proximity: “Find all points within 10 km of X”
• Containment: “Find all objects inside polygon Y”
• Intersection: “Find all roads crossing river Z”

Unlike a normal index on a text or number field, spatial indexes must consider geometry and spatial relationships between multiple coordinates.

R-Tree – The Workhorse of Spatial Indexing

Invented by: Antonin Guttman (1984)
Used in: PostGIS, Oracle Spatial, SQLite/SpatiaLite, MySQL GIS

Concept:

• Stores Minimum Bounding Rectangles (MBRs) that enclose each geometry.
• Index structure is hierarchical — MBRs are nested.
• Queries quickly eliminate branches of the tree that don’t intersect the search area.

SQL Example: PostGIS R-Tree with GiST

-- Create an R-Tree index on the geometry column 
CREATE INDEX idx_parks_geom ON parks USING GIST (geom);

-- Create an R-Tree index on the geometry column 
CREATE INDEX idx_parks_geom ON parks USING GIST (geom);

Python Example — Query using R-Tree in Shapely + Rtree

from rtree import index from shapely.geometry 
import Point 

# Create an R-tree index 
idx = index.Index() 
points = [Point(0,0), Point(1,1), Point(5,5)] 
for i, p in enumerate(points): 
  idx.insert(i, (p.x, p.y, p.x, p.y))
   
# Query for points within a bounding box r
esults = list(idx.intersection((0,0,2,2))) 
print(results) # [0, 1]

from rtree import index from shapely.geometry 
import Point 

# Create an R-tree index 
idx = index.Index() 
points = [Point(0,0), Point(1,1), Point(5,5)] 
for i, p in enumerate(points): 
  idx.insert(i, (p.x, p.y, p.x, p.y))
   
# Query for points within a bounding box r
esults = list(idx.intersection((0,0,2,2))) 
print(results) # [0, 1]

QuadTree – Divide and Conquer

Used in: GIS tile rendering, game engines, image processing, Google Maps tiling

Concept:

• Recursively divides space into four quadrants (NE, NW, SE, SW).
• Efficient for datasets where points are unevenly distributed.
• Works well for point data, less so for complex polygons.

Python Example — QuadTree Point Search:

class QuadTree: 
  def __init__(self, bounds, capacity): 
    self.bounds = bounds # (xmin, ymin, xmax, ymax) 
    self.capacity = capacity 
    self.points = [] 
    self.divided = False 
  
  def insert(self, point): 
    if not self._in_bounds(point): 
      return False 
    if len(self.points) < self.capacity: self.points.append(point) 
      return True 
    else: 
      if not self.divided: 
        self._subdivide() 
    return 
      (
        self.nw.insert(point) or 
        self.ne.insert(point) or 
        self.sw.insert(point) or 
        self.se.insert(point)
      ) 
        
    def _in_bounds(self, point): 
      x, y = point 
      xmin, ymin, xmax, ymax = self.bounds 
      return xmin <= x <= xmax and ymin <= y <= ymax 
      
    def _subdivide(self): 
      xmin, ymin, xmax, ymax = self.bounds 
      xmid = (xmin + xmax) / 2 
      ymid = (ymin + ymax) / 2 
      self.nw = QuadTree((xmin, ymid, xmid, ymax), self.capacity) 
      self.ne = QuadTree((xmid, ymid, xmax, ymax), self.capacity) 
      self.sw = QuadTree((xmin, ymin, xmid, ymid), self.capacity) 
      self.se = QuadTree((xmid, ymin, xmax, ymax), self.capacity) 
      self.divided = True

class QuadTree: 
  def __init__(self, bounds, capacity): 
    self.bounds = bounds # (xmin, ymin, xmax, ymax) 
    self.capacity = capacity 
    self.points = [] 
    self.divided = False 
  
  def insert(self, point): 
    if not self._in_bounds(point): 
      return False 
    if len(self.points) < self.capacity: self.points.append(point) 
      return True 
    else: 
      if not self.divided: 
        self._subdivide() 
    return 
      (
        self.nw.insert(point) or 
        self.ne.insert(point) or 
        self.sw.insert(point) or 
        self.se.insert(point)
      ) 
        
    def _in_bounds(self, point): 
      x, y = point 
      xmin, ymin, xmax, ymax = self.bounds 
      return xmin <= x <= xmax and ymin <= y <= ymax 
      
    def _subdivide(self): 
      xmin, ymin, xmax, ymax = self.bounds 
      xmid = (xmin + xmax) / 2 
      ymid = (ymin + ymax) / 2 
      self.nw = QuadTree((xmin, ymid, xmid, ymax), self.capacity) 
      self.ne = QuadTree((xmid, ymid, xmax, ymax), self.capacity) 
      self.sw = QuadTree((xmin, ymin, xmid, ymid), self.capacity) 
      self.se = QuadTree((xmid, ymin, xmax, ymax), self.capacity) 
      self.divided = True

Geohash – Encoding Location as a String

Used in: Elasticsearch, MongoDB, DynamoDB, Redis

Concept:

• Encodes lat/long into a single base32 string.
• Nearby locations share common prefixes.
• Ideal for sharding/distribution and proximity search.

Python Example — Geohash Encoding

import geohash 
lat, lon = 38.8977, -77.0365 

# White House 
code = geohash.encode(lat, lon, precision=7) 
print(code) # e.g., "dqcjqcp"

import geohash 
lat, lon = 38.8977, -77.0365 

# White House 
code = geohash.encode(lat, lon, precision=7) 
print(code) # e.g., "dqcjqcp"

Space-Filling Curves – Hilbert & Z-order

Concept:

• Map 2D coordinates into 1D values while preserving locality.
• Used in Google S2 Geometry library, Apache HBase, and Cassandra.

Why useful:

• Works well for distributed storage systems.
• Reduces index fragmentation for range queries.

Spatial indexing is the reason modern GIS queries return results in milliseconds instead of hours.
Choosing the right index depends on data shape, query patterns, and storage architecture.

If your dataset is growing into billions of points or polygons, you’ll almost always want:

• R-Tree for general-purpose geometry queries.
• Geohash or Hilbert curves for distributed/cloud-native systems.

References

1. Guttman, A. (1984). R-Trees: A Dynamic Index Structure for Spatial Searching. ACM SIGMOD.
2. Samet, H. (1990). The Design and Analysis of Spatial Data Structures. Addison-Wesley.
3. PostGIS Docs — Indexing: https://postgis.net/workshops/postgis-intro/indexing.html
4. Google S2 Geometry Library: https://s2geometry.io/