0
# Spatial Indexing
1

2
R-tree spatial indexing for efficient spatial queries on large collections of geometries. The STRtree (Sort-Tile-Recursive tree) provides fast spatial queries by organizing geometries into a hierarchical bounding box structure.
3

4
## Capabilities
5

6
### STRtree Index
7

8
High-performance spatial index for geometry collections.
9

10
```python { .api }
11
class STRtree:
12
    def __init__(self, geometries, node_capacity=10):
13
        """
14
        Create spatial index from geometry collection.
15
        
16
        Parameters:
17
        - geometries: array-like of geometries to index
18
        - node_capacity: maximum number of geometries per tree node
19
        """
20
    
21
    @property
22
    def geometries(self):
23
        """Array of indexed geometries."""
24
    
25
    def query(self, geometry, predicate=None):
26
        """
27
        Query index for geometries that satisfy spatial predicate.
28
        
29
        Parameters:
30
        - geometry: query geometry
31
        - predicate: spatial predicate ('intersects', 'within', 'contains', etc.)
32
                    If None, returns geometries with overlapping bounding boxes
33
        
34
        Returns:
35
        ndarray: indices of matching geometries
36
        """
37
    
38
    def query_items(self, geometry, predicate=None):
39
        """
40
        Query index and return actual geometries instead of indices.
41
        
42
        Parameters:
43
        - geometry: query geometry
44
        - predicate: spatial predicate
45
        
46
        Returns:
47
        ndarray: matching geometries
48
        """
49
    
50
    def nearest(self, geometry, max_distance=None, return_distance=False, exclusive=False):
51
        """
52
        Find nearest geometries to query geometry.
53
        
54
        Parameters:
55
        - geometry: query geometry
56
        - max_distance: maximum search distance (None for unlimited)
57
        - return_distance: if True, return distances along with geometries
58
        - exclusive: if True, exclude query geometry from results
59
        
60
        Returns:
61
        ndarray or tuple: nearest geometry indices, optionally with distances
62
        """
63
```
64

65
**Usage Example:**
66

67
```python
68
import shapely
69
import numpy as np
70

71
# Create a collection of geometries to index
72
np.random.seed(42)
73
points = shapely.points(np.random.rand(1000, 2) * 100)
74
polygons = [shapely.Point(x, y).buffer(2) for x, y in np.random.rand(100, 2) * 100]
75

76
# Create spatial index
77
tree = shapely.STRtree(polygons)
78
print(f"Indexed {len(tree.geometries)} polygons")
79

80
# Query geometries intersecting a region
81
query_region = shapely.box(20, 20, 30, 30)
82
intersecting_indices = tree.query(query_region, predicate='intersects')
83
print(f"Found {len(intersecting_indices)} intersecting polygons")
84

85
# Get actual geometries instead of indices
86
intersecting_geoms = tree.query_items(query_region, predicate='intersects')
87
print(f"Retrieved {len(intersecting_geoms)} geometry objects")
88

89
# Find nearest polygons to a point
90
query_point = shapely.Point(50, 50)
91
nearest_indices = tree.nearest(query_point)
92
print(f"Nearest polygon index: {nearest_indices[0] if len(nearest_indices) > 0 else 'None'}")
93

94
# Get nearest with distance
95
nearest_with_dist = tree.nearest(query_point, return_distance=True)
96
if len(nearest_with_dist[0]) > 0:
97
    index, distance = nearest_with_dist[0][0], nearest_with_dist[1][0]
98
    print(f"Nearest polygon at index {index}, distance {distance:.2f}")
99
```
100

101
### Spatial Predicates for Queries
102

103
Available predicates for spatial index queries.
104

105
```python { .api }
106
class BinaryPredicate:
107
    intersects = 1      # Geometries that intersect query
108
    within = 2          # Geometries within query geometry
109
    contains = 3        # Geometries that contain query geometry
110
    overlaps = 4        # Geometries that overlap query
111
    crosses = 5         # Geometries that cross query
112
    touches = 6         # Geometries that touch query
113
    covers = 7          # Geometries that cover query
114
    covered_by = 8      # Geometries covered by query
115
    contains_properly = 9  # Geometries that properly contain query
116
```
117

118
**Usage Example:**
119

120
```python
121
import shapely
122

123
# Create test geometries
124
large_polygon = shapely.box(0, 0, 20, 20)
125
small_polygons = [
126
    shapely.box(5, 5, 8, 8),      # inside
127
    shapely.box(15, 15, 25, 25),  # overlapping
128
    shapely.box(25, 25, 30, 30),  # outside
129
    shapely.box(-5, 10, 5, 15),   # touching edge
130
]
131

132
# Index the small polygons
133
tree = shapely.STRtree(small_polygons)
134

135
# Different query types
136
within_results = tree.query(large_polygon, predicate='within')
137
overlaps_results = tree.query(large_polygon, predicate='overlaps')  
138
touches_results = tree.query(large_polygon, predicate='touches')
139

140
print(f"Polygons within large polygon: {len(within_results)}")
141
print(f"Polygons overlapping large polygon: {len(overlaps_results)}")
142
print(f"Polygons touching large polygon: {len(touches_results)}")
143
```
144

145
### Performance Optimization
146

147
Spatial indexing provides significant performance improvements for large datasets.
148

149
**Usage Example:**
150

151
```python
152
import shapely
153
import numpy as np
154
import time
155

156
# Create large dataset
157
np.random.seed(42)
158
n_geoms = 10000
159
centers = np.random.rand(n_geoms, 2) * 1000
160
radii = np.random.rand(n_geoms) * 5 + 1
161
polygons = [shapely.Point(x, y).buffer(r) for (x, y), r in zip(centers, radii)]
162

163
print(f"Created {len(polygons)} test polygons")
164

165
# Query region
166
query_box = shapely.box(400, 400, 600, 600)
167

168
# Method 1: Brute force (no index)
169
start_time = time.time()
170
brute_force_results = []
171
for i, poly in enumerate(polygons):
172
    if query_box.intersects(poly):
173
        brute_force_results.append(i)
174
brute_force_time = time.time() - start_time
175

176
# Method 2: Using spatial index
177
start_time = time.time()
178
tree = shapely.STRtree(polygons)
179
index_build_time = time.time() - start_time
180

181
start_time = time.time()
182
indexed_results = tree.query(query_box, predicate='intersects')
183
index_query_time = time.time() - start_time
184

185
print(f"Brute force: {len(brute_force_results)} results in {brute_force_time:.3f}s")
186
print(f"Index build time: {index_build_time:.3f}s")
187
print(f"Index query: {len(indexed_results)} results in {index_query_time:.3f}s")
188
print(f"Speedup: {brute_force_time / index_query_time:.1f}x (excluding build time)")
189
print(f"Results match: {set(brute_force_results) == set(indexed_results)}")
190
```
191

192
### Advanced Usage Patterns
193

194
Complex spatial queries and analysis patterns.
195

196
**Usage Example:**
197

198
```python
199
import shapely
200
import numpy as np
201

202
# Create urban planning scenario: buildings and points of interest
203
np.random.seed(42)
204

205
# Buildings (rectangles)
206
building_coords = np.random.rand(500, 2) * 1000
207
buildings = [
208
    shapely.box(x-5, y-5, x+5, y+5) 
209
    for x, y in building_coords
210
]
211

212
# Points of interest
213
poi_coords = np.random.rand(100, 2) * 1000
214
pois = shapely.points(poi_coords)
215

216
# Create indexes
217
building_tree = shapely.STRtree(buildings)
218
poi_tree = shapely.STRtree(pois)
219

220
# Analysis: Find buildings near each POI
221
analysis_results = []
222
for i, poi in enumerate(pois):
223
    # 50-unit buffer around POI
224
    poi_buffer = shapely.buffer(poi, 50)
225
    
226
    # Find nearby buildings
227
    nearby_building_indices = building_tree.query(poi_buffer, predicate='intersects')
228
    
229
    analysis_results.append({
230
        'poi_index': i,
231
        'poi_location': (poi.x, poi.y),
232
        'nearby_buildings': len(nearby_building_indices),
233
        'building_indices': nearby_building_indices
234
    })
235

236
# Summary statistics
237
nearby_counts = [r['nearby_buildings'] for r in analysis_results]
238
print(f"POI analysis complete:")
239
print(f"Average buildings within 50 units: {np.mean(nearby_counts):.1f}")
240
print(f"Max buildings near one POI: {np.max(nearby_counts)}")
241
print(f"POIs with no nearby buildings: {sum(1 for c in nearby_counts if c == 0)}")
242

243
# Find POI clusters: POIs within 100 units of each other
244
cluster_threshold = 100
245
poi_clusters = []
246

247
for i, poi in enumerate(pois):
248
    cluster_buffer = shapely.buffer(poi, cluster_threshold)
249
    nearby_poi_indices = poi_tree.query(cluster_buffer, predicate='intersects')
250
    
251
    # Remove self from results
252
    nearby_poi_indices = nearby_poi_indices[nearby_poi_indices != i]
253
    
254
    if len(nearby_poi_indices) > 0:
255
        poi_clusters.append({
256
            'center_poi': i,
257
            'cluster_size': len(nearby_poi_indices) + 1,  # Include self
258
            'cluster_members': nearby_poi_indices
259
        })
260

261
print(f"Found {len(poi_clusters)} POI clusters")
262
if poi_clusters:
263
    largest_cluster = max(poi_clusters, key=lambda x: x['cluster_size'])
264
    print(f"Largest cluster has {largest_cluster['cluster_size']} POIs")
265
```