A Rusty Geospatial Adventure

Laurențiu Nicola

Motivation 🌳

I'm not from Bucharest, been here since COVID
This city isn't green enough
Can we measure it, though?
We'll try, using 🛰️!

Copernicus Programme

Joint Earth Observation program between EU, ESA, EUMETSAT, ECMWF, JRC, EEA, EMSA and others
Satellites, ground stations and sensors
The data is available to anyone, for free and under an open license

Sentinel-2 Constellation

Constellation of three satellites
~5 days revisit period (but beware of clouds!)
~110 × 110 km² scenes in a grid
We're in tile 35TMK now

Sentinel-2 Multi-Spectral Instrument

4 bands @ 10 m (RGB, NIR)
6 bands @ 20 m (NIR, SWIR)
3 bands @ 60 m (aerosols, water vapor, cirrus)

Source: Copernicus Sentinel-2A Calibration and Products Validation Status

Recognizing Trees Vegetation from Space

Vegetation is green, but so can be water, roofs, shadows &c.
Why are plants green?
What is green anyway?
Sunlight is white and has all the visible colors
Green things reflect green and absorb the other colors

Chlorophyll and the NDVI

Vegetation is green because of chlorophyll
Chlorophyll absorbs blue and red
Chlorophyll reflects green and infrared
Look for things that absorb red and reflect infrared
Normalized Difference Vegetation Index (NDVI) $NDVI = \frac{NIR - Red}{NIR + Red}$

Interlude: True Color Composite

Sentinel-2 true color (RGB) composite, 35TMK, 2024-06-23

Interlude: False Color Composite

Sentinel-2 false color (NIR-R-G) composite, 35TMK, 2024-06-23

NDVI

Shows photosynthetically active vegetation
From −1 to 1, but negative values are not interesting
Bare soils $\sim 0.2$ , sparse vegetation $\sim 0.3$ , dense vegetation $> 0.5$

Crop Phenology

Source: TimeMatch: Unsupervised Cross-Region Adaptation by Temporal Shift Estimation

NDVI (3 Date Composite)

Red / magenta: high NDVI in spring, early crops
Green: high NDVI in summer, summer crops
Yellow: high NDVI in spring and summer, trees
Blue / cyan: high NDVI in summer and autumn, irrigated crops
Black: not vegetation

Sentinel-2 NDVI composite, 35TMK, 2024-04-14, 2024-07-08, 2024-10-26

Sentinel-2 Product Structure

                        S2B_MSIL2A_20241024T091939_N0511_R093_T35TMK_20241024T105557.SAFE
                        ├── GRANULE
                        │   └── L2A_T35TMK_A039873_20241024T091939
                        │       └── IMG_DATA
                        │           ├── R10m
                        │           │   ├── T35TMK_20241024T091939_B02_10m.jp2
                        │           │   ├── T35TMK_20241024T091939_B03_10m.jp2
                        │           │   ├── T35TMK_20241024T091939_B04_10m.jp2
                        │           │   └── T35TMK_20241024T091939_B08_10m.jp2
                        │           └── R20m
                        │               ├── T35TMK_20241024T091939_B05_20m.jp2
                        │               ├── T35TMK_20241024T091939_B06_20m.jp2
                        │               ├── T35TMK_20241024T091939_B07_20m.jp2
                        │               ├── T35TMK_20241024T091939_B11_20m.jp2
                        │               ├── T35TMK_20241024T091939_B12_20m.jp2
                        │               ├── T35TMK_20241024T091939_B8A_20m.jp2
                        │               └── T35TMK_20241024T091939_SCL_20m.jp2
                        └── MTD_MSIL2A.xml

Lots of files elided, but the reflectance bands are under R10m and R20m
Reflectance is the measure of the proportion of incoming radiation that gets reflected by a material
In Sentinel-2 products, multipled by 10 000 and with a 1000 offset added
Look in the MTD_MSIL2A.xml file to get the -1000 offsets
SCL is a mask (clouds, cloud shadows, snow, water &c.)

Why Rust?

High-level, you can use it like Python, nobody will complain
Easy to use, cf. NameError after running for half a day
Fast, on-par with C, even naïve code might be fast enough
Memory-safe, no SIGSEGVs, no data races
Good libraries, and can call into any C library
Good tooling and documentation

Reading Rasters in Rust

We use GDAL, the geospatial I/O library
New (GDAL 3.10) dataset thread-safety support, coming soon to Rust
Can read regions of the input, to save memory

                                use gdal::{Dataset, GdalOpenFlags::*};

let open_flags = GDAL_OF_RASTER | GDAL_OF_THREAD_SAFE;
let ds = Dataset::open_ex(
    path,
    DatasetOptions {
        open_flags,
        ..Default::default()
    },
)?;

let mut data = ds
    .rasterband(1)?
    .read_as(window, window_size, window_size, None)?
    .to_array()?;
data.mapv_inplace(|x| x + offset);


fn ndvi(b4: i16, b8: i16, scl: u8) -> f32 {
    if matches!(scl, 0 | 1 | 6 | 8 | 9 | 10 | 11) {
        return f32::NAN;
    }
    let (b4, b8) = (b4 as f32, b8 as f32);
    if b8 + b4 == 0.0 {
        f32::NAN
    } else {
        (b8 - b4) / (b8 + b4).clamp(-1.0, 1.0)
    }
}

Raster Layouts

Generally, rasters use one of two pixel data layouts:
- Tiles
- Strips: rows written in order
Tiled layouts are usually more efficient because they allow processing of rectangular regions
This is especially true for visualization
Sentinel-2 uses JPEG 2000 (tiled, small, but very slow to decode)

A Rough Sketch

Query products for desired time interval (e.g. CDSE STAC API)
Take rectangular blocks from each input
- first block from each input
- second block from each input
- …
Read the blocks on a thread pool, send on a channel (queue)
Single receiver thread reads dequeues blocks, computes mean NDVI, writes to output
Don't compute the mean directly, but rather the sum and value count
This means we can process the read regions immediately

Some Results

On my PC, 2 minutes for 215 products over one year
18.5 GB on disk, 58.3 GB uncompressed
6-8 GB RAM, depending on region size and number of threads
Basically all spent in OpenJPEG, but can saturate 32 cores
On cloud VM, April - November, 100 products, 11 minutes per tile stack
3 for processing, 8 minutes to open the files 💀
After cropping, mosaicking, conversion to 8-bit and compression: 3.4 GB

Demo

Potential Improvements

Median NDVI (can't compute online)
Combine with different indices (e.g. EVI)
Multi-date classification
Use a proper dataset for classification (e.g. ESA CCI-LC)

GeoRust

Umbrella organization covering different Rust GIS projects
Geometry: geo, geos^*, robust
Spatial indexing (k-NN query): rstar
Coordinate transforms: proj^*, geodesy^†
Formats: gdal^*, netcdf^*, gpx, geojson, wkt, wkb, kml, osm, transit, shapefile, STAC^‡, PgSTAC^‡, GeoTIFF^‡
Utilities: geozero, geocoding, geohash, polyline
Other tools: rinex