How We Delivered “Fields of The World” with RasterFlow: A Planetary-Scale GeoAI Pipeline

Inspect the published feature mosaic

The produced feature mosaic is published as a Zarr store on Source Cooperative. Therefore, we can inspect structure, bands, and a spatial subset directly.

import xarray as xr
import rasterix

mosaic_ds = xr.open_zarr(
    "https://data.source.coop/ftw/global-data/features/zarr/alpha/global.zarr"
).pipe(rasterix.assign_index)
mosaic_ds

<xarray.Dataset> Size: 502TB
Dimensions:      (time: 2, band: 10, y: 1566049, x: 4007517)
Coordinates:
  * time         (time) datetime64[ns] 16B 2024-01-01 2025-01-01
  * band         (band) object 80B 's2med_harvest:B02' ... 's2med_planting:N_...
  * y            (y) float64 13MB 83.75 83.75 83.75 ... -56.93 -56.93 -56.93
  * x            (x) float64 32MB -180.0 -180.0 -180.0 ... 180.0 180.0 180.0
    spatial_ref  int64 8B ...
Data variables:
    variables    (time, band, y, x) float32 502TB dask.array<chunksize=(1, 1, 4096, 4096), meta=np.ndarray>
Indexes:
  ┌ x        RasterIndex (crs=None)
  └ y

Persisting an intermediate global Zarr mosaic lets us query any region without reprojecting or resampling.

Training data for the PRUE model is also aligned to EPSG:4326, so feature and prediction stores remain grid-compatible for side-by-side analysis. For data layout details, see Appendix: Zarr internals.

For the examples below, we use the same Japan AOI shown later in the PMTiles preview so the raster and vector artifacts are easy to compare. For more information on xarray, see Appendix: Xarray notes.

import matplotlib.pyplot as plt
import numpy as np

bbox = {
    "x": slice(140.24, 140.34),
    "y": slice(37.36, 37.26),
}
subset = mosaic_ds.sel(x=bbox["x"], y=bbox["y"])
geo_aspect = 1 / np.cos(np.deg2rad(float(subset.y.mean())))

fig, axes = plt.subplots(1, 2, figsize=(9, 9), constrained_layout=True)
for ax, band_prefix, title in [
    (axes[0], "s2med_planting", "Planting RGB (PMTiles AOI, 2025)"),
    (axes[1], "s2med_harvest", "Harvest RGB (PMTiles AOI, 2025)"),
]:
    subset["variables"].sel(
        time="2025-01-01",
        band=[f"{band_prefix}:B04", f"{band_prefix}:B03", f"{band_prefix}:B02"],
    ).plot.imshow(ax=ax, robust=True)
    ax.set_aspect(geo_aspect)
    ax.set_title(title)

plt.show()

See additional information in Appendix: Zarr internals and Appendix: Xarray notes.

Inspect the published prediction mosaic

The prediction Zarr is also published on Source Cooperative. Let’s verify shape, bands, and value ranges.

The shape, bounds, and time remain aligned with the feature mosaic; however, the bands now represent predicted classes.

pred_ds = xr.open_zarr(
    "https://data.source.coop/ftw/global-data/predictions/zarr/alpha/global.zarr"
).pipe(rasterix.assign_index)
pred_ds

<xarray.Dataset> Size: 151TB
Dimensions:      (time: 2, band: 3, y: 1566049, x: 4007517)
Coordinates:
  * time         (time) datetime64[ns] 16B 2024-01-01 2025-01-01
  * band         (band) <U20 240B 'non_field_background' ... 'field_boundaries'
  * y            (y) float64 13MB 83.75 83.75 83.75 ... -56.93 -56.93 -56.93
  * x            (x) float64 32MB -180.0 -180.0 -180.0 ... 180.0 180.0 180.0
    spatial_ref  int64 8B ...
Data variables:
    variables    (time, band, y, x) float32 151TB dask.array<chunksize=(1, 1, 2048, 2048), meta=np.ndarray>
Indexes:
  ┌ x        RasterIndex (crs=None)
  └ y

Let’s plot the three labels for the same PMTiles-area AOI so the raster predictions line up with the vector preview later in the post:

import numpy as np

bbox = {
    "x": slice(140.24, 140.34),
    "y": slice(37.36, 37.26),
}
pred_subset = pred_ds.sel(x=bbox["x"], y=bbox["y"])
geo_aspect = 1 / np.cos(np.deg2rad(float(pred_subset.y.mean())))

grid = pred_subset["variables"].plot.imshow(col="band", row="time", figsize=(12, 8))
for ax in grid.axes.flat:
    ax.set_aspect(geo_aspect)

plt.show()

Zarr metadata for those curious is available in Appendix: Zarr internals.

Building Feature COGs

Internally, build_mosaic orchestrates four key steps:

1. Seasonal DOY filtering — candidate scenes are filtered by day-of-year windows derived from latitude-aware seasonal heuristics.
2. Quality masking — cloud/quality masks suppress invalid observations.
3. Greedy scene selection — we iteratively pick scenes that maximize uncovered valid area so we reach broad coverage with a bounded scene count.
4. Temporal compositing — selected scenes are collapsed into median composites for downstream model features.

The current DOY heuristic and greedy objective are intentionally conservative and production-stable, but there is room for improvement (for example, region-specific phenology priors, dynamic cloud climatology, and cost-aware objective tuning).

We also produce a COG index, which is used with GDAL’s GTI driver to build the underlying global mosaic.

import geopandas as gpd
import fsspec

with fsspec.open("https://data.source.coop/ftw/global-data/features/cogs/alpha/index.parquet") as f:
    index = gpd.read_parquet(f)
index.head(3)

Mosaicking notes (GTI + Ray)

After feature COGs are written, we build a global virtual mosaic index and then materialize aligned arrays on the target grid. This relies on GDAL’s GTI driver as the index layer and mosaicking.

Grid configuration for this run:

• CRS: EPSG:4326
• Resolution: 8.983119e-5 degrees (~10 m at equator)
• Resampling: cubic

The GTI-backed approach keeps COG indexing separate from downstream chunk/block execution and maps cleanly to distributed scheduling.

Zarr internals

Zarr is the storage layer that makes the global feature and prediction mosaics usable as products rather than just intermediate files.

Two layout choices matter here:

• Logical chunks are the units we want readers and compute stages to address.
• Shards are the larger storage objects we actually write to object storage.

We split them on purpose. Small logical chunks keep regional reads and retry boundaries precise, while larger shards prevent the store from exploding into billions of tiny objects. Larger shards also lower overall S3 PUT and LIST costs and reduce request-rate pressure.

For the feature store, the logical chunk is (1, 1, 4096, 4096) and shards are written at (1, 1, 12288, 12288), so each object packs a 3 x 3 tile of logical chunks. For the prediction store, the logical chunk is (1, 1, 2048, 2048) and shards are written at (1, 3, 8192, 8192), so each object packs all three class bands across a 4 x 4 spatial tile of logical chunks.

The tradeoff is deliberate: larger shards improve object-store efficiency and reduce listing overhead, but smaller logical chunks still give us better partial reads, safer retries, and more control over Ray execution geometry.

The feature Zarr mosaic metadata is as follows:

import zarr

group = zarr.open("https://data.source.coop/ftw/global-data/features/zarr/alpha/global.zarr")
group["variables"].metadata.__dict__

{'shape': (2, 10, 1566049, 4007517),
 'data_type': Float32(endianness='little'),
 'chunk_grid': RegularChunkGrid(chunk_shape=(1, 1, 12288, 12288)),
 'chunk_key_encoding': DefaultChunkKeyEncoding(separator='/'),
 'codecs': (ShardingCodec(chunk_shape=(1, 1, 4096, 4096), codecs=(BytesCodec(endian=<Endian.little: 'little'>), ZstdCodec(level=0, checksum=False)), index_codecs=(BytesCodec(endian=<Endian.little: 'little'>), Crc32cCodec()), index_location=<ShardingCodecIndexLocation.end: 'end'>),),
 'dimension_names': ('time', 'band', 'y', 'x'),
 'fill_value': np.float32(nan),
 'attributes': {'coordinates': 'spatial_ref', '_FillValue': 'AAAAAAAA+H8='},
 'storage_transformers': (),
 'extra_fields': {}}

The prediction zarr mosaic metadata is as follows:

group = zarr.open("https://data.source.coop/ftw/global-data/predictions/zarr/alpha/global.zarr")
group["variables"].metadata.__dict__

{'shape': (2, 3, 1566049, 4007517),
 'data_type': Float32(endianness='little'),
 'chunk_grid': RegularChunkGrid(chunk_shape=(1, 3, 8192, 8192)),
 'chunk_key_encoding': DefaultChunkKeyEncoding(separator='/'),
 'codecs': (ShardingCodec(chunk_shape=(1, 1, 2048, 2048), codecs=(BytesCodec(endian=<Endian.little: 'little'>), ZstdCodec(level=0, checksum=False)), index_codecs=(BytesCodec(endian=<Endian.little: 'little'>), Crc32cCodec()), index_location=<ShardingCodecIndexLocation.end: 'end'>),),
 'dimension_names': ('time', 'band', 'y', 'x'),
 'fill_value': np.float32(nan),
 'attributes': {'coordinates': 'spatial_ref', '_FillValue': 'AAAAAAAA+H8='},
 'storage_transformers': (),
 'extra_fields': {}}

Xarray notes

We use rasterix to attach analytical coordinates lazily so opening the global feature store does not require eagerly reading large coordinate vectors every time.

Why this matters:

• x length: 1,566,049
• y length: 4,007,517
• At float64 precision, coordinates alone are roughly (1,566,049 + 4,007,517) * 8 ~= 45 MB before reading any model features.

Operationally, our xarray pattern is:

1. Open lazily from object storage.
2. Slice first (.sel) to bound I/O.
3. Materialize only when plotting/validating.

This keeps exploratory analysis responsive while preserving direct comparability between feature and prediction stores on the same grid.

The Inference Config

The InferenceConfig is the contract between model packaging and distributed execution. It tells RasterFlow how to read features, how large each inference window should be, what runtime to launch, and how overlapping predictions should be merged back into the output mosaic.

In practice, the fields fall into a few groups:

• model_path and actor select the model artifact and the inference runtime implementation.
• patch_size and clip_size define the read window and how much border area is discarded to suppress edge artifacts.
• features and labels define the exact tensor interface between the Zarr store and the model.
• device and max_batch_size control throughput versus memory pressure on the target worker.
• merge_mode defines how overlapping patch predictions are blended into the final raster.

from rasterflow_remote.data_models import MergeModeEnum, InferenceActorEnum

{
    # path to the pt2 archive on hugging face
    'model_path': 'https://huggingface.co/wherobots/prue-pt2/resolve/main/prue-efnetb7.pt2',
    # an enum used to specify how to handle the model for inference
    'actor': InferenceActorEnum.SEMANTIC_SEGMENTATION_PYTORCH,
    # patch size determines the size of the inputs provided to the model
    'patch_size': 256,
    # clip size determines the size we clip off the borders to reduce edge effects during inference
    'clip_size': 32,
    # device to run inference on, e.g. 'cuda' or 'cpu'
    'device': 'cuda',
    # list of features in order as required by the model. These names map to the input Zarr band labels.
    'features': ['s2med_harvest:B04',
    's2med_harvest:B03',
    's2med_harvest:B02',
    's2med_harvest:B08',
    's2med_planting:B04',
    's2med_planting:B03',
    's2med_planting:B02',
    's2med_planting:B08'],
    # the named labels for the predicted classes. These will be the output band labels in the output Zarr.
    'labels': ['non_field_background', 'field', 'field_boundaries'],
    # the maximum batch size assuming our default A10 GPU with 24GB of memory.
    'max_batch_size': 128,
    # the method to use when merging overlapping predictions to reduce edge effects at the patch level
    'merge_mode': MergeModeEnum.WEIGHTED_AVERAGE
}

How Ray is used for inference

We chose Ray as the data-plane engine because the workload is embarrassingly parallel at block level, but still needs data-aware scheduling and actor pools for GPU inference.

Ray is used as a bounded block-processing engine over sharded Zarr storage.

Key mapping model:

• Logical chunks define storage math and indexability.
• Shards reduce object count by packing multiple chunks per object.
• Ray blocks define execution granularity and may span multiple chunks.

This decoupling lets us tune compute batch shape without rewriting storage layout. In the data pipeline, Ray Data uses Arrow-backed transport so tabular metadata/state can move across stages with minimal copying.

Practical takeaway: throughput scales with actor parallelism while per-task memory remains tied to block geometry, not global mosaic size.

Distributed inference (wherobots-rasterflow)

The InferenceConfig above feeds directly into block planning, actor setup, and overlap-aware merging. Internally, predict_mosaic orchestrates three phases:

1. Block planning — derive execution blocks from chunk metadata and resource limits.
2. Read/Infer/Write pipeline — execute overlap-aware reads, actor inference, and deterministic writes.
3. Retry-safe commit — persist block outputs with stable target regions.

# conceptual : chunk/shard metadata -> Ray blocks

chunk_meta = read_zarr_chunk_metadata(store)
blocks = plan_ray_blocks(chunk_meta, target_rows_per_block)

ray_ds = ray.data.from_items(blocks)
result = (
    ray_ds
    .map(read_zarr_block_with_overlap)
    .map(run_inference_actor_pool)
    .map(clip_and_write_block)
)

Our distributed approach to vectorizing

vectorize_mosaic operates at block level so geometry generation scales linearly with partition count rather than requiring one global polygonization pass.

High-level flow:

1. Read prediction blocks from Zarr.
2. Apply thresholding per block.
3. Polygonize per block.
4. Emit GeoParquet partitions and merge metadata.

This keeps memory bounded and allows retries/reprocessing for individual partitions.