cuRaster

cuRaster is a high-performance Python library for GPU-accelerated raster processing. It reads GeoTIFF files (locally or directly from S3), executes band-math algebra on the GPU, and optionally reprojects, clips, and streams results — all through a clean, lazy pipeline API.

Table of Contents

• Installation
• Requirements (building from source)
• Quick Start
• API Reference
  • curaster.open()
  • Chain.algebra()
  • Chain.clip()
  • Chain.reproject()
  • Chain.get_info()
  • Chain.save_local()
  • Chain.save_s3()
  • Chain.to_memory()
  • Chain.iter_begin()
  • Chain.focal()
  • Chain.terrain()
  • Chain.texture()
  • Chain.zonal_stats()
  • curaster.open_stack() / StackChain
  • RasterResult
  • ChunkQueue
• Examples
• Building from Source
• License

Installation

pip install curaster

Wheels are pre-built for Linux and Windows, Python 3.9–3.13. A compatible NVIDIA GPU and driver must be present at runtime (the CUDA runtime is bundled in the wheel).

Requirements (building from source)

Quick Start

import curaster

# Compute NDVI and save to a local GeoTIFF
curaster.open("landsat.tif") \
    .algebra("(B5 - B4) / (B5 + B4)") \
    .save_local("ndvi.tif")

API Reference

curaster.open(path)

Open a GeoTIFF and return a lazy Chain. No GPU work happens here.

chain = curaster.open("input.tif")
chain = curaster.open("s3://my-bucket/data/scene.tif")   # S3 direct-read
chain = curaster.open("/vsis3/my-bucket/data/scene.tif") # GDAL vsis3 URI

Returns Chain

Chain.algebra(expression)

Append a band-math operation. Bands are referenced as B1, B2, … (1-indexed).

chain.algebra("(B5 - B4) / (B5 + B4)")      # NDVI
chain.algebra("B1 * 0.0001")                  # Scale factor
chain.algebra("(B3 + B2 + B1) / 3")          # Visible mean
chain.algebra("B4 > 0.3")                     # Boolean mask (1.0 or 0.0)
chain.algebra("(B4 > 0.2) * B4")             # Apply mask conditionally

Supported operators: + - * / > < >= <= == !=

Returns a new Chain (original is unmodified)

Chain.clip(geojson)

Clip the output to a polygon. Pixels outside the polygon are set to zero.

import json

aoi = json.dumps({
    "type": "Polygon",
    "coordinates": [[[10.0, 52.0], [11.0, 52.0], [11.0, 53.0], [10.0, 53.0], [10.0, 52.0]]]
})

chain.algebra("(B5 - B4) / (B5 + B4)").clip(aoi)

Returns a new Chain

Chain.reproject(target_crs, ...)

Reproject the output to a different coordinate reference system.

chain.reproject("EPSG:4326")                             # Auto pixel size
chain.reproject("EPSG:3857", res_x=10.0, res_y=10.0)   # Fixed 10 m resolution
chain.reproject("EPSG:4326", resampling="nearest")      # Nearest-neighbour

# Fixed output extent (in target CRS units)
chain.reproject(
    "EPSG:4326",
    res_x=0.0001, res_y=0.0001,
    te_xmin=9.5, te_ymin=51.5,
    te_xmax=10.5, te_ymax=52.5
)

Returns a new Chain

Chain.get_info()

Return metadata for the output raster without executing the pipeline.

info = curaster.open("scene.tif").reproject("EPSG:4326").get_info()
print(info)
# {'width': 4096, 'height': 3072, 'geotransform': [...], 'crs': 'GEOGCS[...]'}

Returns dict with keys width, height, geotransform (list of 6 floats), crs (WKT string)

Chain.save_local(path, verbose=False)

Execute the pipeline and write a Float32 tiled GeoTIFF to disk.

curaster.open("scene.tif") \
    .algebra("(B5 - B4) / (B5 + B4)") \
    .save_local("ndvi.tif", verbose=True)

Chain.save_s3(s3_path, verbose=False)

Execute the pipeline and upload the result directly to S3.
AWS credentials must be set via environment variables (AWS_ACCESS_KEY_ID, AWS_SECRET_ACCESS_KEY, AWS_DEFAULT_REGION).

curaster.open("scene.tif") \
    .algebra("(B5 - B4) / (B5 + B4)") \
    .save_s3("/vsis3/my-bucket/output/ndvi.tif")

Chain.to_memory(verbose=False)

Execute and return all pixels as a RasterResult object. Raises RuntimeError if the result would exceed 75 % of available RAM — use iter_begin() for large rasters.

result = curaster.open("scene.tif") \
    .algebra("B1 * 0.0001") \
    .to_memory()

import numpy as np
arr = result.data()          # numpy array, shape (height, width), dtype float32
print(arr.mean(), arr.std())
print(result.width, result.height, result.proj)

Returns RasterResult

Chain.iter_begin(buf_chunks=4)

Start background execution and return a ChunkQueue for memory-efficient streaming. Each chunk covers a horizontal strip of the output.

queue = curaster.open("huge_scene.tif") \
    .algebra("(B5 - B4) / (B5 + B4)") \
    .iter_begin(buf_chunks=8)

while True:
    chunk = queue.next()
    if chunk is None:
        break
    # chunk = {'y_offset': int, 'width': int, 'height': int, 'data': np.ndarray}
    process(chunk["data"], chunk["y_offset"])

Returns ChunkQueue

RasterResult

Returned by to_memory().

ChunkQueue

Returned by iter_begin(). Processing runs on a background thread.

Each chunk dict has keys: y_offset (int), width (int), height (int), data (np.ndarray float32).

Examples

NDVI — local file, save to disk

import curaster

curaster.open("landsat8_sr.tif") \
    .algebra("(B5 - B4) / (B5 + B4)") \
    .save_local("ndvi.tif", verbose=True)

S3 direct-read → S3 write

import curaster, os

# Credentials are read from the environment automatically
curaster.open("s3://my-bucket/scenes/LC08_2024_scene.tif") \
    .algebra("(B5 - B4) / (B5 + B4)") \
    .save_s3("/vsis3/my-bucket/output/ndvi.tif")

Clip to area of interest

import curaster, json

aoi = json.dumps({
    "type": "Polygon",
    "coordinates": [[[13.3, 52.4], [13.5, 52.4], [13.5, 52.6], [13.3, 52.6], [13.3, 52.4]]]
})

curaster.open("sentinel2.tif") \
    .algebra("(B8 - B4) / (B8 + B4)") \
    .clip(aoi) \
    .save_local("ndvi_berlin.tif")

Reproject to WGS84 with fixed resolution

import curaster

curaster.open("utm_scene.tif") \
    .algebra("(B4 - B3) / (B4 + B3)") \
    .reproject("EPSG:4326", res_x=0.0001, res_y=0.0001) \
    .save_local("ndvi_wgs84.tif")

Full pipeline: S3 → algebra → clip → reproject → S3

import curaster, json

aoi = json.dumps({
    "type": "Polygon",
    "coordinates": [[[10.0, 52.0], [11.0, 52.0], [11.0, 53.0], [10.0, 53.0], [10.0, 52.0]]]
})

curaster.open("s3://my-bucket/raw/sentinel2.tif") \
    .algebra("(B8 - B4) / (B8 + B4)") \
    .clip(aoi) \
    .reproject("EPSG:4326", res_x=0.0001, res_y=0.0001) \
    .save_s3("/vsis3/my-bucket/processed/ndvi_reprojected.tif")

Inspect output metadata before running

import curaster

info = curaster.open("scene.tif") \
    .reproject("EPSG:4326", res_x=0.0001) \
    .get_info()

print(f"Output will be {info['width']} × {info['height']} pixels")
print(f"CRS: {info['crs'][:60]}...")

Load into numpy / xarray

import curaster
import numpy as np
import xarray as xr

result = curaster.open("scene.tif") \
    .algebra("(B5 - B4) / (B5 + B4)") \
    .to_memory()

arr = result.data()   # shape (H, W), dtype float32
arr[arr == -9999.0] = np.nan

gt = curaster.open("scene.tif").get_info()["geotransform"]
xcoords = gt[0] + np.arange(result.width)  * gt[1]
ycoords = gt[3] + np.arange(result.height) * gt[5]

da = xr.DataArray(arr, dims=["y", "x"], coords={"x": xcoords, "y": ycoords})
print(da)

Streaming large rasters chunk-by-chunk

import curaster
import numpy as np

output = np.zeros((10000, 10000), dtype=np.float32)

queue = curaster.open("massive_scene.tif") \
    .algebra("(B5 - B4) / (B5 + B4)") \
    .iter_begin(buf_chunks=6)

while True:
    chunk = queue.next()
    if chunk is None:
        break
    y0 = chunk["y_offset"]
    h  = chunk["height"]
    output[y0 : y0 + h, :] = chunk["data"]

print("Done. Mean NDVI:", output.mean())

Boolean / conditional expression

import curaster

# Mask pixels where NIR reflectance > 0.3, zero elsewhere
curaster.open("scene.tif") \
    .algebra("(B5 > 0.3) * B5") \
    .save_local("nir_high_mask.tif")

# Multi-band composite score
curaster.open("scene.tif") \
    .algebra("(B5 - B4) / (B5 + B4) + (B3 - B2) / (B3 + B2)") \
    .save_local("composite_score.tif")

Chain.focal(stat, radius=3, shape="square", clamp_border=True)

Apply a moving-window focal statistic.

curaster.open("dem.tif") \
    .focal("mean", radius=5) \
    .save_local("dem_smoothed.tif")

curaster.open("dem.tif") \
    .focal("median", radius=3, shape="circle") \
    .save_local("dem_median.tif")

Returns a new Chain

Chain.terrain(metrics=["slope"], unit="degrees", sun_azimuth=315.0, sun_altitude=45.0, method="horn")

Compute terrain derivatives from a DEM.

curaster.open("dem.tif") \
    .terrain(["slope", "aspect", "hillshade"]) \
    .save_local("terrain.tif")

# All supported metrics:
curaster.open("dem.tif") \
    .terrain(["slope", "aspect", "hillshade", "tri", "tpi",
              "roughness", "prof_curv", "plan_curv", "total_curv"],
             unit="degrees") \
    .save_local("terrain_all.tif")

Output is a multi-band GeoTIFF with one band per metric in the order given.

Returns a new Chain

Chain.texture(features=[], window=11, levels=32, direction_mode="average", log_scale=False, val_min=0.0, val_max=0.0)

Compute GLCM (Grey-Level Co-occurrence Matrix) Haralick texture features.

curaster.open("sar.tif") \
    .texture(["contrast", "homogeneity", "entropy"], window=15, levels=64, log_scale=True) \
    .save_local("texture.tif")

# All 18 features, 4-direction average
curaster.open("image.tif") \
    .texture(window=11, levels=32) \
    .save_local("texture_full.tif")

Returns a new Chain

Chain.zonal_stats(geojson, stats=["mean", "std", "min", "max", "count", "sum"], band=1, verbose=False)

Compute per-polygon zonal statistics over any raster. Terminal operation — returns results immediately.

import json, curaster

aoi = json.dumps({
    "type": "MultiPolygon",
    "coordinates": [...]
})

results = curaster.open("ndvi.tif").zonal_stats(aoi, stats=["mean", "std", "min", "max"])
for r in results:
    print(r.zone_id, r.mean, r.std_dev, r.min, r.max)

Returns list[ZoneResult] where each has .zone_id, .count, .mean, .std_dev, .min_val, .max_val, .sum

curaster.open_stack(files) / StackChain

Open a multi-temporal stack of aligned GeoTIFF files and reduce them to a single raster.

import curaster

stack = curaster.open_stack(["s2_20230601.tif", "s2_20230701.tif", "s2_20230801.tif"])

# Temporal difference (last - first)
stack.temporal("diff").save_local("diff.tif")

# Linear trend slope (change per scene)
stack.temporal("trend", time_values=[0.0, 30.0, 60.0]).save_local("trend.tif")

# Mean of all scenes
stack.temporal("mean").save_local("mean.tif")

# You can chain further operations after temporal reduction
stack.temporal("diff") \
    .clip(aoi_geojson) \
    .save_local("diff_clipped.tif")

StackChain.temporal() Returns a Chain (can chain algebra, clip, save_local, etc.)

All scenes in the stack must have the same width, height, and CRS. Use .reproject() on each Chain before stacking if misaligned.

Performance & Benchmarks

The following benchmarks demonstrate cuRaster's performance for both Local storage and direct S3 reads.

Hardware Specifications:

• Instance: AWS g4dn.xlarge
• Compute: 4 vCPUs, 16 GiB RAM (Intel Xeon 2.5 GHz)
• GPU: 1x NVIDIA T4 Tensor Core (16 GiB VRAM, Compute Capability 7.5)
• On-Demand Cost: €0.563 / hour

We test across varying raster sizes:

• S: 2048 × 2048
• M: 4096 × 4096
• L: 8192 × 8192
• XL: 16384 × 16384
• XXL: 32768 × 8192

Each cell displays the Processing Time alongside the Estimated Compute Cost for that single operation.

Local GeoTIFF Operations

These tests read files directly from the local NVMe SSD.

S3 Direct-Read Operations

These tests read data dynamically over the network from an AWS S3 bucket using GDAL's virtual file system and libcurl with HTTP Range requests.

# 1. Install Python build dependencies
pip install build pybind11 scikit-build-core setuptools_scm

# 2. Install system libraries (Ubuntu/Debian)
sudo apt install libgdal-dev libzstd-dev libssl-dev libcurl4-openssl-dev libomp-dev

# 3. Build the wheel
python -m build --wheel

# 4. Install the built wheel
pip install wheelhouse/*.whl

Or build directly with CMake for development:

mkdir build && cd build
cmake .. -Dpybind11_DIR=$(python -c "import pybind11; print(pybind11.get_cmake_dir())")
make -j$(nproc)

# Copy the .so into your working directory
cp curaster*.so ..

License

See LICENSE for details.