End-to-End Geospatial Processing with EORST

// 2026-05-07 · Updated 2026-05-11 · 4 min read

This post walks through a complete geospatial processing pipeline using rss_core and eorst. We’ll query satellite imagery from DEA spanning two UTM zones, and show how eorst handles the automatic reprojection.

New to remote sensing? Start with Getting Started with Geospatial Rust — covers satellite platforms, spectral bands, indices, and why cloud detection matters.

The Setup

Add the dependencies to your Cargo.toml:

[dependencies]
rss_core = "0.4"
eorst = "0.3"
anyhow = "1.0"
chrono = "0.4"
log = "0.4"
env_logger = "0.11"

Step 1: Query Satellite Imagery

We’ll query DEA for Sentinel-2 data over a region in Queensland that spans two UTM zones (tiles 54KYB and 55KBS). Using a bounding box instead of individual tiles:

use rss_core::{
    query::ImageQueryBuilder,
    qvf::Collection,
    DEA,
    utils::{Cmp, Intersects},
};
use chrono::NaiveDate;
use anyhow::Result;
use std::path::PathBuf;

fn main() -> Result<()> {
    // tiles 54KYB and 55KBS are adjacent in Queensland
    // combined bbox spans UTM zones 54 and 55:
    // - 54KYB: x 142.92 → 143.90 (UTM zone 54)
    // - 55KBS: x 144.48 → 145.08 (UTM zone 55)
    let bbox = Intersects::Bbox(vec![
        142.92, -21.69,  // xmin, ymin
        145.08, -20.78,  // xmax, ymax
    ]);

    let query = ImageQueryBuilder::new(
        DEA,
        Collection::Sentinel2,
        bbox,
    )
    .start_date(NaiveDate::parse_from_str("2024-01-01", "%Y-%m-%d").unwrap())
    .end_date(NaiveDate::parse_from_str("2024-06-01", "%Y-%m-%d").unwrap())
    .cloudcover((Cmp::Less, 10))
    .canonical_bands(["red", "nir"])  // Semantic band names!
    .build();

    let tmp_dir = PathBuf::from("/tmp/satellite_data");
    std::fs::create_dir_all(&tmp_dir)?;

    // Download the imagery via STAC
    let stac = query.get(&tmp_dir, None, None)?;
    println!("Downloaded {} items", stac.items.len());

    Ok(())
}

Notice .canonical_bands(["red", "nir"]) — this is the semantic band resolution. rss_core translates "red" to provider-specific band names:

DEA: nbart_red
Microsoft Planetary Computer: band4
Element84: B04

This means your code isn’t tied to one provider. Switch catalogs without changing your query.

Why a Bounding Box?

Sentinel-2 uses a tiling system based on MGRS (Military Grid Reference System). Our region of interest sits across two adjacent tiles that span different UTM zones:

Tile	UTM Zone	Longitude Range
54KYB	54	142.92°E → 143.90°E
55KBS	55	144.48°E → 145.08°E

A bounding box naturally captures all overlapping tiles — no need to list them explicitly. And here’s the magic: eorst handles the automatic reprojection. When data comes from different UTM zones, it aligns everything to a common CRS automatically.

Step 2: Build a RasterDataset

use eorst::{RasterDatasetBuilder, RasterDataset, types::BlockSize};

fn build_dataset(stac: &rss_core::stac::FeatureCollection) -> Result<RasterDataset<u16>> {
    let rds: RasterDataset<u16> = RasterDatasetBuilder::from_stac_query(stac)
        .block_size(BlockSize { rows: 2048, cols: 2048 })
        .build();

    println!("{}", rds);
    Ok(rds)
}

The RasterDataset is a virtual view over your data — it handles the multi-temporal stacking, reprojection, and block-based access automatically.

Step 3: Process with the Apply API

Now the key part — processing with .apply() and the Select trait:

use eorst::{Select, RasterDataset, RasterDatasetBuilder};
use ndarray::{Array4, Zip};
use std::path::PathBuf;

fn process(rds: &RasterDataset<u16>) -> Result<()> {
    let output = PathBuf::from("/tmp/ndvi_output.tif");
    
    rds.apply(
        |block| {
            // Select layers by semantic name - no more hardcoded indices!
            let red = block.select_layers(&["red"])?;
            let nir = block.select_layers(&["nir"])?;

            // Compute NDVI: (NIR - Red) / (NIR + Red)
            // Range: -1.0 (water) to +1.0 (healthy vegetation)
            // Scaled to i16: -10000 to +10000
            // See: https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index
            let mut ndvi = Array4::<i16>::zeros(red.data.raw_dim());
            Zip::from(&mut ndvi)
                .and(&red.data)
                .and(&nir.data)
                .for_each(|out, &r, &n| {
                    let val = (n as f32 - r as f32) / (n as f32 + r as f32 + 1e-10);
                    *out = (val * 10000.0) as i16;
                });

            Ok(ndvi)
        },
        4,  // number of threads
        &output,
    )?;
    
    println!("NDVI computed: {:?}", output);
    Ok(())
}

The Select trait lets you pick bands by name. This is cleaner than .slice(s![0, .., ..]) — your code says what it means.

Other Indices You Can Compute

The same pattern works for other spectral indices:

// NDWI (Normalized Difference Water Index) - https://www.opendatacube.org/ndwi
// Highlights open water: positive = water, built-up = near 0 or negative
let green = block.select_layers(&["green"])?;
let ndwi = (&green.data - &nir.data) / (&green.data + nir.data + 1e-10);

// NBR (Normalized Burn Ratio) - https://landscapetoolbox.org/remote-sensing-methods/normalized-burn-ratio/
// Detects burn scars: high positive = healthy vegetation, low/negative = burned areas
let swir = block.select_layers(&["swir_2"])?;
let nbr = (&nir.data - &swir.data) / (&nir.data + &swir.data + 1e-10);

Step 4: Export

The .apply() method handles export automatically. Results are written to a GeoTIFF with proper georeferencing — load it in QGIS, ArcGIS, or Python.

What Just Happened?

Query: We asked DEA for Sentinel-2 data over a bounding box spanning tiles 54KYB and 55KBS (two UTM zones!)
Download: rss_core fetched the data via STAC and cached it locally
Build: eorst created a virtual RasterDataset — automatically handling the reprojection from two CRSs
Apply: We computed NDVI using semantic band selection
Export: Results written to a single aligned GeoTIFF

Why This Matters

The canonical band system means your code isn’t tied to one provider. Query DEA today, switch to Planetary Computer tomorrow — the same .canonical_bands(["red", "nir"]) call works everywhere.

The Select trait makes your code readable. Instead of guessing band indices, you name what you want — and the compiler helps you catch typos.

How fast is eorst compared to Python? See Benchmarking Rust vs Python+Dask for NDVI: 23× Faster, and Why That Matters for Carbon — a head-to-head comparison of eorst/Rust against Python+rioxarray on real satellite imagery. Simple NDVI is 2.9× faster; add FMask cloud masking and it’s 23× faster.