eo-processor¶

High-performance Rust (PyO3) UDFs for Earth Observation (EO) processing with Python bindings. Fast spectral indices, temporal statistics, masking utilities, and spatial distance functions.

Documentation¶

Full documentation is hosted online: https://bnjam.dev/eo-processor/

The site contains:

Quick start guides and worked examples (NDVI, EVI, compositing, masking)
Complete API reference with function signatures, expected input shapes/dtypes, and return types
CLI usage and examples for batch processing and PNG preview generation
Tutorials for integrating with XArray / Dask and for building reproducible benchmarks
Developer & contribution notes (how to add Rust UDFs, register Python bindings, test, and create type stubs)
Guidance on building from source, running the benchmark harness, and regenerating coverage badges
Release notes / changelog and citation information

Overview¶

eo-processor accelerates common remote sensing computations using safe Rust (no unsafe) exposed via PyO3. All public functions interoperate with NumPy and can be embedded in XArray / Dask pipelines. Rust kernels release Python’s GIL; multi-core parallelism (via Rayon) is leveraged for selected operations (larger temporal aggregations, pairwise distances).

Focus areas:

Spectral & change-detection indices
Temporal statistics & median compositing (1D–4D stacks)
Masking & data quality filtering (value / range / SCL / invalid sentinels)
Pairwise spatial distances (utility layer)
Benchmark harness for reproducible performance measurements

Key Features¶

Rust-accelerated numerical kernels (float64 internal, stable results)
Automatic dimensional dispatch (1D / 2D for spectral indices, 1D–4D for temporal/masking)
Change detection support (ΔNDVI, ΔNBR)
Flexible masking utilities (exact values, ranges, SCL codes)
Median, mean, sample standard deviation over time axis
Pairwise distance functions (Euclidean, Manhattan, Chebyshev, Minkowski)
Type stubs (__init__.pyi) for IDE / mypy
Benchmark script with optional NumPy baseline comparison
Pure CPU, no external network or storage side-effects in core path

Installation¶

PyPI (standard)¶

pip install eo-processor

Optional extras for array ecosystem:

pip install eo-processor[dask]

Using `uv`¶

uv venv
source .venv/bin/activate
uv pip install eo-processor

From Source¶

Requirements:

Python 3.9+
Rust toolchain (rustup recommended)
maturin for building the extension module

git clone https://github.com/BnJam/eo-processor.git
cd eo-processor

pip install maturin
maturin develop --release        # build & install in-place
# or wheel:
maturin build --release
pip install target/wheels/*.whl

Quick Start¶

import numpy as np
from eo_processor import ndvi, ndwi, evi, normalized_difference

nir   = np.array([0.8, 0.7, 0.6])
red   = np.array([0.2, 0.1, 0.3])
blue  = np.array([0.1, 0.05, 0.08])
green = np.array([0.35, 0.42, 0.55])

print(ndvi(nir, red))               # NDVI
print(ndwi(green, nir))             # NDWI
print(evi(nir, red, blue))          # EVI
print(normalized_difference(nir, red))

All inputs may be any numeric NumPy dtype (int/uint/float); internal coercion to float64.

API Summary¶

Function	Purpose
`normalized_difference(a, b)`	Generic normalized difference `(a - b) / (a + b)` with near-zero denominator safeguard
`ndvi(nir, red)`	Normalized Difference Vegetation Index
`ndwi(green, nir)`	Normalized Difference Water Index
`ndsi(green, swir1)`	Normalized Difference Snow Index `(Green - SWIR1)/(Green + SWIR1)`
`evi(nir, red, blue)` / `enhanced_vegetation_index(...)`	Enhanced Vegetation Index (G*(NIR - Red)/(NIR + C1Red - C2Blue + L))
`evi2(nir, red)`	Enhanced Vegetation Index 2-band variant `2.5(NIR - Red)/(NIR + 2.4Red + 1)`
`savi(nir, red, L=0.5)`	Soil Adjusted Vegetation Index `(NIR - Red)/(NIR + Red + L) * (1 + L)`
`nbr(nir, swir2)`	Normalized Burn Ratio `(NIR - SWIR2)/(NIR + SWIR2)`
`ndmi(nir, swir1)`	Normalized Difference Moisture Index `(NIR - SWIR1)/(NIR + SWIR1)`
`nbr2(swir1, swir2)`	Normalized Burn Ratio 2 `(SWIR1 - SWIR2)/(SWIR1 + SWIR2)`
`gci(nir, green)`	Green Chlorophyll Index `(NIR / Green) - 1` (division guard)
`delta_ndvi(pre_nir, pre_red, post_nir, post_red)`	Change in NDVI `(NDVI_pre - NDVI_post)`
`delta_nbr(pre_nir, pre_swir2, post_nir, post_swir2)`	Change in NBR `(NBR_pre - NBR_post)`
`median(arr, skip_na=True)`	Temporal median (time axis) with NaN skipping
`composite(arr, method="median")`	Compositing convenience (currently median only)
`temporal_mean(arr, skip_na=True)`	Mean across time axis
`temporal_std(arr, skip_na=True)`	Sample standard deviation (n-1) across time
`temporal_sum(arr, skip_na=True)`	Sum across time axis
`temporal_composite(arr, weights, skip_na=True)`	Weighted temporal composite
`moving_average_temporal(arr, window, skip_na=True, mode="same")`	Sliding window mean (same/valid edge modes, NaN skip/propagate)
`moving_average_temporal_stride(arr, window, stride, skip_na=True, mode="same")`	Strided moving average (downsampled temporal smoothing)
`pixelwise_transform(arr, scale=1.0, offset=0.0, clamp_min=None, clamp_max=None)`	Per-pixel linear transform with optional clamping
`linear_regression(y)`	Simple linear regression (slope, intercept, residuals) on 1D array
`bfast_monitor(stack, dates, ...)`	BFAST Monitor change detection on time series
`random_forest_train(features, labels, ...)`	Train a Random Forest model
`random_forest_predict(model_json, features)`	Predict using a trained Random Forest model
`complex_classification(blue, green, ...)`	Multi-band classification workflow
`haralick_features(data, window_size, ...)`	Calculate Haralick texture features (GLCM)
`zonal_stats(values, zones)`	Calculate statistics per zone
`binary_dilation(input, kernel_size)`	Binary morphological dilation
`binary_erosion(input, kernel_size)`	Binary morphological erosion
`binary_opening(input, kernel_size)`	Binary morphological opening
`binary_closing(input, kernel_size)`	Binary morphological closing
`euclidean_distance(points_a, points_b)`	Pairwise Euclidean distances
`manhattan_distance(points_a, points_b)`	Pairwise L1 distances
`chebyshev_distance(points_a, points_b)`	Pairwise L∞ distances
`minkowski_distance(points_a, points_b, p)`	Pairwise L^p distances (p ≥ 1)
`mask_vals(arr, values=None, fill_value=None, nan_to=None)`	Mask exact codes, optional fill & NaN normalization
`replace_nans(arr, value)`	Replace all NaNs with `value`
`mask_out_range(arr, min_val=None, max_val=None, fill_value=None)`	Mask values outside `[min, max]`
`mask_in_range(arr, min_val=None, max_val=None, fill_value=None)`	Mask values inside `[min, max]`
`mask_invalid(arr, invalid_values, fill_value=None)`	Mask list of sentinel values (e.g., `0, -9999`)
`mask_scl(scl, keep_codes=None, fill_value=None)`	Mask Sentinel‑2 SCL codes, keeping selected classes
`mask_with_scl(data, scl, mask_codes=None, fill_value=None)`	Apply SCL-based mask to data array

Temporal dimension expectations:

1D: (time,)
2D: (time, band)
3D: (time, y, x)
4D: (time, band, y, x)

Distance functions: input shape (N, D) and (M, D) → output (N, M) (O(N*M) memory/time).

Spectral & Change Detection Indices¶

All indices auto-dispatch 1D vs 2D arrays (matching shapes required).

NDVI¶

(NIR - Red) / (NIR + Red) Interpretation (approximate):

< 0: water / snow
0.0–0.2: bare soil / built surfaces
0.2–0.5: sparse to moderate vegetation
0.5: healthy dense vegetation

NDWI¶

(Green - NIR) / (Green + NIR)

0.3: open water (often 0.4–0.6)
0.0–0.3: moist vegetation / wetlands
< 0.0: dry vegetation / soil

NDSI¶

(Green - SWIR1) / (Green + SWIR1)

Often > 0.3: likely snow/ice
Near 0: mixed pixels / uncertain surfaces
< 0.0: common for non-snow surfaces

EVI¶

G * (NIR - Red) / (NIR + C1*Red - C2*Blue + L) (MODIS constants: G=2.5, C1=6.0, C2=7.5, L=1.0) Improves sensitivity over high biomass & reduces soil/atmospheric noise vs NDVI.

EVI2¶

2.5 * (NIR - Red) / (NIR + 2.4*Red + 1) 2-band EVI variant often used when the blue band is unavailable.

SAVI¶

(NIR - Red) / (NIR + Red + L) * (1 + L) Typical L=0.5. Larger L for sparse vegetation (bright soil), smaller for dense vegetation.

NBR¶

(NIR - SWIR2) / (NIR + SWIR2) Used for burn severity. Compare pre/post via ΔNBR.

NDMI¶

(NIR - SWIR1) / (NIR + SWIR1) Moisture / canopy water content indicator.

NBR2¶

(SWIR1 - SWIR2) / (SWIR1 + SWIR2) Highlights moisture & thermal differences; complementary to NBR/NDMI.

GCI¶

(NIR / Green) - 1 Chlorophyll proxy; division by near-zero guarded to avoid instability.

Change Detection¶

ΔNDVI = NDVI_pre - NDVI_post ΔNBR = NBR_pre - NBR_post Positive ΔNDVI: vegetation loss. Positive ΔNBR: burn severity increase.

Masking Utilities¶

Rust-accelerated preprocessing helpers for quality filtering.

Function	Notes
`mask_vals`	Exact equality masking (codes → `fill_value` or NaN) + optional NaN normalization
`replace_nans`	Force all NaNs to a scalar
`mask_out_range`	Mask outside interval
`mask_in_range`	Mask inside interval
`mask_invalid`	Shorthand for common invalid sentinels
`mask_scl`	Keep only selected Sentinel‑2 SCL classes
`mask_with_scl`	Apply SCL-based mask directly to data array

Example:

import numpy as np
from eo_processor import mask_vals, replace_nans, mask_out_range, mask_scl, mask_with_scl

scl = np.array([4,5,6,8,9])  # vegetation, vegetation, water, cloud (med), cloud (high)
clear = mask_scl(scl, keep_codes=[4,5,6])   # -> [4., 5., 6., nan, nan]
# mask data where SCL is cloud/high cloud (8, 9)
masked_data = mask_with_scl(np.ones(5), scl, mask_codes=[8, 9]) # -> [1., 1., 1., nan, nan]

ndvi = np.array([-0.3, 0.1, 0.8, 1.2])
valid = mask_out_range(ndvi, min_val=-0.2, max_val=1.0)  # -> [nan,0.1,0.8,nan]

arr = np.array([0, 100, -9999, 50])
clean = mask_vals(arr, values=[0, -9999])  # -> [nan,100.,nan,50.]
filled = replace_nans(clean, -9999.0)      # -> [-9999.,100.,-9999.,50.]

Morphological Operations¶

Binary morphological operations for 2D arrays (e.g. masks).

Function	Purpose
`binary_dilation(input, kernel_size)`	Dilate features (expand white regions)
`binary_erosion(input, kernel_size)`	Erode features (shrink white regions)
`binary_opening(input, kernel_size)`	Erosion followed by dilation (remove noise)
`binary_closing(input, kernel_size)`	Dilation followed by erosion (fill holes)

All operations assume the input is a 2D array where values > 0 are treated as True/foreground. The structuring element is a square kernel of size kernel_size.

Temporal Statistics & Compositing¶

Median, mean, sum, and standard deviation across time axis (skip NaNs optional):

import numpy as np
from eo_processor import temporal_mean, temporal_std, temporal_sum, median, temporal_composite

cube = np.random.rand(12, 256, 256)  # (time, y, x)
mean_img  = temporal_mean(cube)      # (256, 256)
std_img   = temporal_std(cube)       # (256, 256)
sum_img   = temporal_sum(cube)       # (256, 256)
median_img = median(cube)

# Weighted composite for 4D arrays (time, bands, y, x)
cube_4d = np.random.rand(5, 4, 256, 256)
weights = np.array([0.1, 0.2, 0.4, 0.2, 0.1])
comp_img = temporal_composite(cube_4d, weights) # (4, 256, 256)

composite(cube, method="median") currently routes to median.

Trend Analysis & Regression¶

eo-processor provides tools for trend analysis, regression, and change detection on time series data.

Function	Purpose
`trend_analysis(y, threshold)`	Detects breaks in a time series by recursively fitting linear models.
`linear_regression(y)`	Simple linear regression on a 1D array (returns slope, intercept, residuals).
`bfast_monitor(stack, dates, ...)`	BFAST Monitor for change detection in time series stacks.

Example:

import numpy as np
from eo_processor._core import trend_analysis
from eo_processor import linear_regression

# Simple linear regression
y_reg = np.array([1.0, 2.1, 2.9, 4.2])
slope, intercept, resid = linear_regression(y_reg)

# Trend Analysis on time series
# Generate some sample data with a break
y = np.concatenate([
    np.linspace(0, 10, 50),
    np.linspace(10, 0, 50)
]) + np.random.normal(0, 0.5, 100)

# Run the trend analysis
segments = trend_analysis(y.tolist(), threshold=5.0)

# Print the results
print("Trend Analysis Results:")
for segment in segments:
    print(
        f"  Start: {segment.start_index}, "
        f"End: {segment.end_index}, "
        f"Slope: {segment.slope:.4f}, "
        f"Intercept: {segment.intercept:.4f}"
    )

# BFAST Monitor (Change Detection)
from eo_processor import bfast_monitor
from datetime import datetime

# stack: (time, y, x), dates: list of python datetime objects matching time axis
# result = bfast_monitor(
#     stack, dates,
#     history_start_date=datetime(2019, 1, 1),
#     monitor_start_date=datetime(2020, 1, 1),
#     h=0.25, alpha=0.05
# )

Classification & Feature Extraction¶

Tools for classification and texture analysis.

Function	Purpose
`random_forest_train(features, labels, ...)`	Train a Random Forest classifier
`random_forest_predict(model_json, features)`	Predict using a trained Random Forest model
`complex_classification(...)`	Multi-band classification workflow
`haralick_features(data, window_size, ...)`	Compute GLCM texture features (Contrast, Homogeneity, etc.)
`zonal_stats(values, zones)`	Compute statistics for defined zones

from eo_processor import haralick_features, zonal_stats
import xarray as xr

# Haralick features on xarray DataArray (Dask-aware)
# data = xr.DataArray(...)
# features = haralick_features(data, window_size=5)

# Zonal statistics
values = np.random.rand(100, 100)
zones = np.random.randint(0, 5, (100, 100))
stats = zonal_stats(values, zones)
# Access results: stats[zone_id].mean, stats[zone_id].sum, etc.

# Random Forest
from eo_processor import random_forest_train, random_forest_predict

# Train (returns JSON model string)
# model_json = random_forest_train(features, labels, n_estimators=100)

# Predict
# predictions = random_forest_predict(model_json, features)

Advanced Temporal & Pixelwise Processing¶

High-performance smoothing and per-pixel transforms for deep temporal stacks and large spatial tiles.

Formulas:

Moving average: MA_t = mean(x_{start..end}) where [start, end] is the window centered (same) or fixed (valid) around t.
Strided moving average: sample MA_{k*stride} for integer k to downsample temporal resolution.
Pixelwise transform: y = clamp(scale * x + offset) (clamping optional).

Example (moving average with edge handling and NaN skipping):

from eo_processor import moving_average_temporal
import numpy as np

series = np.array([1.0, 2.0, 3.0, 4.0])
ma_same  = moving_average_temporal(series, window=3, mode="same")   # length preserved
ma_valid = moving_average_temporal(series, window=3, mode="valid")  # only full windows

3D temporal cube smoothing (deep stack):

cube = np.random.rand(48, 1024, 1024)
smoothed = moving_average_temporal(cube, window=5, mode="same", skip_na=True)

Strided downsampling (reduce temporal resolution):

from eo_processor import moving_average_temporal_stride
downsampled = moving_average_temporal_stride(cube, window=5, stride=4, mode="same")
print(downsampled.shape)  # (ceil(48/4), 1024, 1024)

Pixelwise transform (scale + offset + clamping):

from eo_processor import pixelwise_transform
arr = np.random.rand(2048, 2048)
stretched = pixelwise_transform(arr, scale=1.2, offset=-0.1, clamp_min=0.0, clamp_max=1.0)

Chaining operations (temporal smoothing then per-pixel adjustment):

ma = moving_average_temporal(cube, window=7)
enhanced = pixelwise_transform(ma, scale=1.1, offset=0.05, clamp_min=0.0, clamp_max=1.0)

Performance Notes:

Prefix-sum approach makes moving average O(T) per pixel independent of window size.
Parallelization occurs over spatial/band pixels for 3D/4D arrays.
Strided variant reduces output size for downstream tasks (e.g., model inference, feature extraction).
Pixelwise transforms are single-pass and can be fused with other operations in custom workflows.

Use Cases:

Smoothing noisy temporal reflectance or index stacks prior to trend analysis.
Reducing temporal dimension before ML model training (stride-based smoothing).
Intensity scaling & clamping for visualization or input normalization.

Spatial Distances¶

Pairwise distance matrices:

import numpy as np
from eo_processor import euclidean_distance, manhattan_distance

A = np.random.rand(100, 8)  # (N, D)
B = np.random.rand(250, 8)  # (M, D)

dist_e = euclidean_distance(A, B)    # (100, 250)
dist_l1 = manhattan_distance(A, B)

For large N*M consider spatial indexing or chunking (not implemented).

XArray / Dask Integration¶

import dask.array as da
import xarray as xr
from eo_processor import ndvi

nir_dask  = da.random.random((5000, 5000), chunks=(500, 500))
red_dask  = da.random.random((5000, 5000), chunks=(500, 500))

nir_xr = xr.DataArray(nir_dask, dims=["y", "x"])
red_xr = xr.DataArray(red_dask, dims=["y", "x"])

ndvi_xr = xr.apply_ufunc(
    ndvi,
    nir_xr,
    red_xr,
    dask="parallelized",
    output_dtypes=[float],
)

result = ndvi_xr.compute()

CLI Usage¶

Console script exposed as eo-processor (installed via PyPI):

# Single index
eo-processor --index ndvi --nir nir.npy --red red.npy --out ndvi.npy

# Multiple indices (provide necessary bands)
eo-processor --index ndvi savi ndmi nbr --nir nir.npy --red red.npy --swir1 swir1.npy --swir2 swir2.npy --out-dir outputs/

# Change detection (ΔNBR)
eo-processor --index delta_nbr \
  --pre-nir pre/nir.npy --pre-swir2 pre/swir2.npy \
  --post-nir post/nir.npy --post-swir2 post/swir2.npy \
  --out outputs/delta_nbr.npy

# List supported indices
eo-processor --list

# Apply cloud mask (0=cloud, 1=clear)
eo-processor --index ndvi --nir nir.npy --red red.npy --mask cloudmask.npy --out ndvi_masked.npy

# PNG preview (requires optional Pillow)
eo-processor --index ndvi --nir nir.npy --red red.npy --out ndvi.npy --png-preview ndvi.png

Selected flags:

--savi-l soil brightness factor for SAVI.
--clamp MIN MAX output range clamping.
--allow-missing skip indices lacking required bands instead of error.

Performance¶

Example benchmark (NDVI on a large array):

import numpy as np, time
from eo_processor import ndvi

nir = np.random.rand(5000, 5000)
red = np.random.rand(5000, 5000)

t0 = time.time()
rust_out = ndvi(nir, red)
t_rust = time.time() - t0

t0 = time.time()
numpy_out = (nir - red) / (nir + red)
t_numpy = time.time() - t0

print(f"Rust: {t_rust:.3f}s  NumPy: {t_numpy:.3f}s  Speedup: {t_numpy/t_rust:.2f}x")

Speedups depend on array shape, memory bandwidth, and CPU cores. Use the benchmark harness for systematic comparison.

Benchmark Harness¶

scripts/benchmark.py provides grouped tests:

# Spectral functions (e.g., NDVI, NDWI, EVI, SAVI, NBR, NDMI, NBR2, GCI)
python scripts/benchmark.py --group spectral --height 2048 --width 2048

# Temporal (compare Rust vs NumPy)
python scripts/benchmark.py --group temporal --time 24 --height 1024 --width 1024 --compare-numpy

# Distances
python scripts/benchmark.py --group distances --points-a 2000 --points-b 2000 --point-dim 8

# All groups; write reports
python scripts/benchmark.py --group all --compare-numpy --json-out bench.json --md-out bench.md

Key options:

--functions <list> override group selection.
--compare-numpy baseline timings (speedup > 1.0 ⇒ Rust faster).
--minkowski-p <p> set order (p ≥ 1).
--loops, --warmups repetition control.
--json-out, --md-out artifact outputs.

Test Coverage¶

Regenerate badge after modifying logic/tests:

tox -e coverage
python scripts/generate_coverage_badge.py coverage.xml coverage-badge.svg

Ensure the badge is committed if coverage changes materially.

Contributing¶

Follow repository guidelines (AGENTS.md, copilot instructions). Checklist before proposing a PR:

Implement Rust function(s) (no unsafe)
Register via wrap_pyfunction! in src/lib.rs
Export in python/eo_processor/__init__.py
Add type stubs in python/eo_processor/__init__.pyi
Add tests (tests/test_<feature>.py) including edge cases & NaN handling
Update README (API Summary, examples, formulas)
Run:
- cargo fmt
- cargo clippy -- -D warnings
- cargo test (if Rust tests)
- pytest
- tox -e coverage
- ruff and mypy (if configured)
Update version if public API added (minor bump)
Regenerate coverage badge if changed
Confirm no secrets / large binaries staged

Commit message pattern:

<type>(scope): concise summary

Optional rationale, benchmarks, references

Types: feat, fix, perf, docs, test, chore, build, ci

Example:

feat(indices): add Green Chlorophyll Index (GCI)

Implements 1D/2D dispatch, tests, docs, benchmark entry.

Semantic Versioning¶

Patch: Internal fixes, refactors, docs only
Minor: New functions (backward-compatible)
Major: Breaking changes (signature changes, removals)

Roadmap (Indicative)¶

Additional spectral indices (future: NBR derivatives, custom moisture composites)
Sliding window / neighborhood statistics (mean, variance)
Optional multithread strategies for very large temporal cubes
Expanded masking (boolean predicate composition)
Extended change metrics (ΔNDMI, fractional vegetation cover)

(Items requiring strategic design will request human review before implementation.)

Scientific Citation¶

@software{eo_processor,
  title   = {eo-processor: High-performance Rust UDFs for Earth Observation},
  author  = {Ben Smith},
  year    = {2025},
  url     = {https://github.com/BnJam/eo-processor}
}

License¶

MIT License – see LICENSE.

Disclaimer¶

Core library focuses on computational primitives. It does NOT perform:

Sensor-specific radiometric calibration
Atmospheric correction
CRS reprojection / spatial indexing
Cloud/shadow detection algorithms beyond simple masking
Data acquisition / I/O orchestration

Integrate with domain tools (rasterio, xarray, dask, geopandas) for full pipelines.

Support¶

Open issues for bugs or enhancements. Provide:

Reproducible snippet
Input shapes / dtypes
Expected vs actual output
Benchmark data (if performance-related)

Acknowledgements¶

Built with PyO3, NumPy, ndarray, and Rayon. Thanks to the scientific EO community for standardized index formulations.

Enjoy fast, reproducible Earth Observation processing!

eo-processor¶

Documentation¶

Overview¶

Key Features¶

Installation¶

PyPI (standard)¶

Using uv¶

From Source¶

Quick Start¶

API Summary¶

Spectral & Change Detection Indices¶

NDVI¶

NDWI¶

NDSI¶

EVI¶

EVI2¶

SAVI¶

NBR¶

NDMI¶

NBR2¶

GCI¶

Change Detection¶

Masking Utilities¶

Morphological Operations¶

Temporal Statistics & Compositing¶

Trend Analysis & Regression¶

Classification & Feature Extraction¶

Advanced Temporal & Pixelwise Processing¶

Spatial Distances¶

XArray / Dask Integration¶

CLI Usage¶

Performance¶

Benchmark Harness¶

Test Coverage¶

Contributing¶

Semantic Versioning¶

Roadmap (Indicative)¶

Scientific Citation¶

License¶

Disclaimer¶

Support¶

Acknowledgements¶

Using `uv`¶