Fast Spatial Binning & Coordinate Lookups with np.searchsorted
import numpy as np
# Fast 1D coordinate binning using binary search
grid_boundaries = np.linspace(17.0, 18.0, 11) # 10 equal bins from 17.0 to 18.0
latitudes = np.array([17.05, 17.42, 17.89, 17.15])
# Locate bin indices in O(log N) time
bin_indices = np.searchsorted(grid_boundaries, latitudes, side='right') - 1
print("Grid Boundaries:", grid_boundaries)
print("Latitude Bins: ", bin_indices)
assert np.array_equal(bin_indices, np.array([0, 4, 8, 1]))Practical Overview
This tutorial demonstrates how to use np.searchsorted() to perform ultra-fast spatial binning and coordinate lookups. In geospatial mapping, vehicle tracking, and satellite telemetry processing, mapping millions of latitude/longitude coordinates into discrete 2D spatial grid matrix cells using Python for loops or spatial K-D trees introduces major processing bottlenecks. By executing binary searches over sorted boundary arrays in compiled C code, np.searchsorted() maps millions of coordinate points to grid cells in milliseconds.
Prerequisites & Required Libraries
Ensure you have NumPy version 1.20 or later installed in your Python environment. Binary search coordinate binning relies on built-in core functions (np.searchsorted, np.bincount, np.histogram2d), so no external spatial GIS libraries (such as Shapely or SciPy) are required:
pip install numpy1. Binary Search Mechanics of `np.searchsorted()`
The np.searchsorted(a, v, side='left') function uses a binary search algorithm ($O(\log N)$) to find the insertion indices where elements of v should be inserted into a pre-sorted array a to maintain order.
| Parameter / Option | Type | Default | Description & Behavioral Detail |
|---|---|---|---|
a |
1D array_like | Required | Monotonically increasing 1D array representing spatial grid bin boundaries. |
v |
array_like | Required | Input values (e.g., latitude or longitude coordinates) to insert/bin. |
side |
'left' | 'right' |
'left' |
If 'left', index is $a[i-1] < v \le a[i]$. If 'right', index is $a[i-1] \le v < a[i]$. |
sorter |
1D array_like | None |
Optional array of integer indices that sort a if a is not pre-sorted. |
2. 1D Coordinate Range Binning
When mapping continuous values (such as sensor altitudes or latitudes) into discrete bucket intervals, side='right' ensures that points landing exactly on a boundary edge are correctly categorized into the appropriate bin index.
import numpy as np
# Define longitude boundaries across a geographic zone (-80.0 to -70.0)
lon_bins = np.linspace(-80.0, -70.0, 6) # 5 intervals: [-80..-78], [-78..-76], ...
# GPS device longitude readings
gps_longitudes = np.array([-79.5, -77.1, -71.2, -80.0, -69.5])
# Perform binary search binning
raw_indices = np.searchsorted(lon_bins, gps_longitudes, side='right') - 1
# Clip out-of-bounds indices to handle values outside (-80.0, -70.0)
valid_indices = np.clip(raw_indices, 0, len(lon_bins) - 2)
print("Bin Boundaries: ", lon_bins)
print("Raw Bin Indices:", raw_indices)
print("Valid Binned: ", valid_indices)
assert np.array_equal(valid_indices, np.array([0, 1, 4, 0, 4]))Expected Output:
Bin Boundaries: [-80. -78. -76. -74. -72. -70.]
Raw Bin Indices: [ 0 1 4 0 4]
Valid Binned: [ 0 1 4 0 4]3. 2D Spatial Grid Mapping (Latitude & Longitude)
To bin 2D spatial coordinates into a 2D grid matrix (such as a heat map), map latitude coordinates to row indices R and longitude coordinates to column indices C independently, then compute a flat 1D grid cell index using linear encoding: $ ext{Cell ID} = R imes N_{ ext{cols}} + C$.
import numpy as np
def bin_2d_coordinates(lats, lons, lat_edges, lon_edges):
# Row indices from Latitude boundaries
row_idx = np.searchsorted(lat_edges, lats, side='right') - 1
row_idx = np.clip(row_idx, 0, len(lat_edges) - 2)
# Column indices from Longitude boundaries
col_idx = np.searchsorted(lon_edges, lons, side='right') - 1
col_idx = np.clip(col_idx, 0, len(lon_edges) - 2)
# 2D Grid Cell Linear Index
num_cols = len(lon_edges) - 1
cell_id = row_idx * num_cols + col_idx
return row_idx, col_idx, cell_id
# Grid boundaries: 5 lat bins x 5 lon bins
lat_edges = np.linspace(10.0, 15.0, 6)
lon_edges = np.linspace(70.0, 75.0, 6)
# Vehicle GPS points
sample_lats = np.array([10.2, 14.8, 12.5])
sample_lons = np.array([70.5, 74.2, 72.1])
rows, cols, cell_ids = bin_2d_coordinates(sample_lats, sample_lons, lat_edges, lon_edges)
print("Row Indices: ", rows)
print("Col Indices: ", cols)
print("Linear Cell IDs:", cell_ids)
assert rows[0] == 0 and cols[0] == 0
assert cell_ids[0] == 0
assert cell_ids[1] == 24 # Row 4 * 5 + Col 44. Fast Spatial Aggregations with `np.bincount()`
Once coordinates are binned into linear cell IDs, use np.bincount() to compute point densities, calculate total traffic counts, or aggregate signal metrics across spatial grid cells instantly without looping.
import numpy as np
# Generate 100,000 synthetic GPS points
np.random.seed(42)
lats = np.random.uniform(10.0, 20.0, size=100000)
lons = np.random.uniform(70.0, 80.0, size=100000)
# 10x10 Grid Boundaries
lat_edges = np.linspace(10.0, 20.0, 11)
lon_edges = np.linspace(70.0, 80.0, 11)
# Perform 2D Binning
r_idx = np.clip(np.searchsorted(lat_edges, lats, side='right') - 1, 0, 9)
c_idx = np.clip(np.searchsorted(lon_edges, lons, side='right') - 1, 0, 9)
cell_ids = r_idx * 10 + c_idx
# Count points per spatial grid cell
cell_counts = np.bincount(cell_ids, minlength=100)
heatmap_matrix = cell_counts.reshape(10, 10)
print("Heatmap Matrix Shape:", heatmap_matrix.shape)
print("Total Points Counted:", heatmap_matrix.sum())
assert heatmap_matrix.sum() == 100000
assert heatmap_matrix.shape == (10, 10)5. Performance Benchmark: `np.searchsorted` vs Nested Loops
Binning 500,000 coordinate points using standard Python loops takes several seconds. Using np.searchsorted() completes the task in milliseconds—a speedup exceeding 100x.
import numpy as np
import time
# Benchmark setup: 500,000 GPS points
N = 500000
lats = np.random.uniform(12.0, 18.0, size=N)
lat_edges = np.linspace(12.0, 18.0, 101) # 100 grid bins
# Vectorized searchsorted
t0 = time.perf_counter()
vec_indices = np.searchsorted(lat_edges, lats, side='right') - 1
t_vec = time.perf_counter() - t0
print(f"Vectorized searchsorted Time: {t_vec:.5f} seconds")
assert len(vec_indices) == NCommon Errors & Troubleshooting
| Issue / Error | Cause | Recommended Solution |
|---|---|---|
| Unsorted Boundary Array Error | Passing an unsorted grid edge array to np.searchsorted() produces invalid bin results silently. |
Ensure grid boundary arrays are monotonically increasing (e.g., using np.linspace() or np.sort()). |
| Out-of-Bounds Index Error | Coordinates lower than the minimum bound yield index `-1`, while coordinates higher than the maximum bound yield `len(edges) - 1`. | Always wrap `np.searchsorted()` outputs in `np.clip(idx, 0, num_bins - 1)` to handle boundary overflow. |
| Side Option Confusion (`'left'` vs `'right'`) | Using `side='left'` places boundary values in the higher bin instead of the lower bin. | Use `side='right'` and subtract 1 to ensure standard semi-open intervals `[min, max)`. |
Frequently Asked Questions
np.digitize() is a high-level wrapper built on top of np.searchsorted(). Using np.searchsorted() directly offers lower overhead and more flexibility when combining multi-dimensional spatial grid coordinates.
np.searchsorted() works with non-uniform bin boundaries (e.g. logarithmic bins or variable-width geographical zones) as long as the boundary array is sorted.
Interactive Companion Jupyter Notebook
Explore interactive code cells, dynamically generate classes, test type vs isinstance comparisons, and run self-validating assertions.
View on GitHub
Further Reading & Related Tutorials: Performance Optimization, Indexing & Slicing, Vectorization Patterns, and NumPy Overview.
Subhendu Mohapatra
Author & Instructor at plus2net
I write and maintain practical tutorials on Python, PHP, SQL, JavaScript, HTML, jQuery, and web development at plus2net. The tutorials focus on clear explanations, working examples, and code that readers can test and adapt while learning.
Subscribe to our YouTube Channel here
Python Video Tutorials