Geospatial: monitoring yes, GIS no, geo-AI partly#

A raster is what this engine was built for: a grid, a rule applied to every cell, and a result that is the same everywhere. That makes it a good fit for monitoring and a poor fit for GIS — and a partial fit for geo-AI. The three are worth separating.

Job	Fit
monitoring: threshold alarms, change against a reference, band arithmetic (NDVI-like), area accounting, propagation models (fire, flood, dispersion), Monte Carlo over scenarios	strong
GIS: reprojection, vector geometry, topology, map algebra on mosaics, tiling huge rasters	not provided
geo-AI: inference of a CNN/U-Net/ViT on tiles, training on imagery	partial — small models are possible, convolution layers are not provided

What a monitoring run looks like#

Two synthetic scenes ship with the repository (a baseline and a current acquisition, with cells without data):

python examples/make_geo_rasters.py                 # writes the two Esri ASCII grids
vision info examples/data/geo_current.asc           # resolution, cell size, origin, missing cells
vision run examples/monitor.vis --receipt r.json --sign key.pem

loaded scene.now  <- examples/data/geo_current.asc (ASCII grid, cell 10.0, no data = -9999.0; 6 cells became 0)
alerts = 398    hot_cells = 87    changed_cells = 398

398 cells × 100 m² = 3.98 ha in alarm, and — this is the part that matters — the run produced a signed receipt, which a third party checks without a GPU:

vision verify r.json --device numpy
  MATCH: 6 outputs (maps and measures) bit-identical to the receipt

An alarm set that a regulator, an insurer or the other side of a dispute can re-derive on any machine is not something a raster stack gives you. That is the reason to use this engine for monitoring: not speed, but that the alarm is reproducible and attributable.

What it reads, and what it does not#

Reads	Does not
`.asc` (the Esri ASCII grid GDAL exports from any GeoTIFF), `.npy`, `.csv`, `.png`	reprojection: the coordinate reference system is not read, so bring the raster in the CRS you want and use `cellsize`/origin as parameters
`NODATA_value`, reported and turned into an error until the program declares what to put there	vector data: no polygons, no points, no topology
the cell size and the origin, so an area is a count times the cell area	tiling: the whole raster is resident, so a mosaic is cut into windows outside and processed one at a time

A 96×96 scene is free; on the card measured here a 4096×4096 raster runs at 21 Gpx/s, and the limit is memory: four bytes per cell per channel, so a 20 000×20 000 scene is 1.6 GB per channel and does not fit beside a model.

Geo-AI: what is realistic#

The tensor side has autodiff, matmul, attention and optimisers, so a small model over a tile — a classifier, a regression, a small transformer — can be built and trained deterministically. What is not provided is the machinery geo-AI usually means: convolution and pooling layers, pretrained weights, data augmentation, tiling and stitching of a mosaic.

The honest architecture for a study that needs both:

the model runs in whatever stack the team already uses (PyTorch, ONNX, a hosted service);
the rules, indices and physics run here, where they are reproducible and carry a receipt;
the two are bound together by the report: the model’s output is an input layer, and the result is the auditable part.

Zonal statistics with geopandas and rasterstats#

The vector libraries and this engine meet at one type: a mask. A polygon cannot enter the language (there is no geometry), but the rasterisation of it can: a channel where 1 means “the cell centre is inside the zone”, which is the same convention rasterstats uses with all_touched=False.

# geopandas / rasterio: the polygon becomes a mask
from rasterio.features import rasterize
mask = rasterize([(polygon, 1)], out_shape=(rows, cols), transform=transform, fill=0, dtype="float32")
np.save("examples/data/geo_zones.npy", mask)

// .vis: the statistics inside the zone, with the receipt
measure mean_zone = mean(t) on scene where zone == 1
measure max_zone  = max(t)  on scene where zone == 1
measure hot_zone  = count(t > 55) on scene where zone == 1

Measured on the synthetic scene that ships with the repository (two zones):

zone	area	mean (this engine)	mean (numpy)	hottest cell	cells above 55
1 — square	8.86 ha	19.10752825	19.107527	30.4	0
2 — circle	10.20 ha	36.69917159	36.699173	81.1	87

The two agree to about one part in a million, and not bit for bit — which is the point worth understanding. numpy and rasterstats accumulate in float64 with their own order; this engine accumulates in its own fixed order, which is the same on every device. Neither is more accurate for a mean of 10⁶ cells; only one of them can be re-run on another machine and give the same bits, and that one can carry a signed receipt:

vision verify receipt.json --device numpy
  MATCH: 8 outputs (maps and measures) bit-identical to the receipt

Spatial Monte Carlo: what works, and three things not to promise#

The scenarios N block runs the same initial state N times, and the built-in scenario is the index that makes each realisation different. Measured on a 6 GB card (256×256, 40 steps, one measure per step):

100 scenarios, 268.7 M pixel updates: 93 ms of GPU work
burnt: 4000 measurements | p5 0.296734 | median 0.299522 | p95 0.302375 | min 0.294998 | max 0.304001

So the distribution is there, with percentiles, from one run. Three corrections before this becomes a sales sentence:

1. A fixed seed makes N identical scenarios. rand(x, y, 1) gives the same field for every scenario — measured: 100 scenarios, one value, identical to the seventh digit. The index must enter the draw: rand(x, y, scenario).

2. “Thousands of scenarios in fractions of a second” is not true today. The scenarios run one after the other, each with a full copy of the state. Measured: 100 scenarios = about 0.75 s of work (plus 0.6 s to compile the program), so 1000 scenarios are about 7 s. Respectable, not a fraction of a second.

3. Percentiles are per measure, not per coordinate. p5, median and p95 summarise the scalar measure over the scenarios. A risk map — the percentile of each cell across the scenarios — is the artefact an insurer or a site-selection study actually buys, and it does not exist yet. Until it does, the honest phrasing is “the distribution of the indicator”, not “the risk surface”.

And one phrase to retire: “mathematically unquestionable”. The receipt proves that the computation is reproducible and attributable; it says nothing about whether the model is right. The physical validity comes from the separate validation work (the published references the kit compares against), and the two claims must be kept apart — a receipt on a wrong model is still a wrong model, only a repeatable one.

The 3-D side#

canvas [x, y, z] exists and the conformance kit exercises it. The cost is cubic in memory: 256³ is 67 MB per channel, 512³ is 537 MB, 1024³ is 4.3 GB, so the card decides the resolution. The grid and the operators are provided; the physics (groundwater flow, plume dispersion, mining topology) is the program written on top of them — nothing in the engine knows what water or gas does.

Image models: the convolution is a stencil#

A convolutional layer is, by definition, a weighted sum of neighbours — the same object the language is built from. So a frozen layer can be executed here, with its weights loaded from a file:

param N = 96
canvas img[N, N] edge clamp { g }
canvas w[3, 3] { kx  ky }

kernel sobel on img {
    var gx = 0.0
    var gy = 0.0
    for i in 0 to 3 {
        for j in 0 to 3 {
            let v = img.g[x + i - 1, y + j - 1]
            gx += w.kx[i, j] * v
            gy += w.ky[i, j] * v
        }
    }
    g = sqrt(gx * gx + gy * gy)
}

python examples/make_conv_weights.py     # the 3x3 kernels as .npy
vision compare examples/conv_edge.vis    # every device against the numpy reference

Measured:

what	numpy	RTX 4050	Intel iGPU	llvmpipe	agreement
Sobel 3×3, 96²	4.2 ms	0.4 ms	1.0 ms	0.7 ms	bit-identical
Sobel 3×3, 1024²	164.5 ms	1.5 ms (692 Mpx/s)	7.9 ms	10.6 ms	bit-identical
chain 1→4→4→1, 96²	21.7 ms	1.0 ms	3.2 ms	1.4 ms	bit-identical

110× faster than numpy on one conv layer, 22× on a three-layer chain, with the same bits on every device — the property a computer-vision stack does not offer: one source, an RTX, an Intel iGPU and a CPU, the same output and a receipt for it.

python -m integrations.torch.conv_to_vis --layers 1,4,4,1 --size 96 --out model.vis   # from weights
python -m integrations.torch.check_against_torch          # the same chain in PyTorch (needs torch)

Against PyTorch on the same chain and the same weights: 5.96e-07 maximum relative difference — the float32 summation order of cuDNN against the generated kernels, not an error. Reproducing that comparison found four defects, and each is now a test:

defect	what it looked like
`for i in 0 to 3`	four taps in this language, not three — the window shifted and the numbers were merely plausible
transposing the kernel	unnecessary: the loader already maps an array’s [row, column] to the canvas [x, y]; the transpose mirrored the filter
dropping the `Conv2d` bias	on by default in PyTorch; its absence moved the output by whole units
a `ReLU` after every layer	right for a conv chain, wrong for a logits head — hence `relu=[..., False]`

Chains that change size: max-pool#

A classifier pools, and pooling halves the canvas — so a chain that pools needs a canvas per stage, not the two alternating ones. emit_staged writes that: conv (with bias and ReLU) and max-pool, each stage its own canvas.

python -m integrations.torch.conv_to_vis        # conv-only chains (alternating canvases, less memory)
# or, from Python, for a chain with pooling:
#   emit_staged([{"kind": "conv", "cout": 4, ...}, {"kind": "pool", "k": 2}, {"kind": "conv", "cout": 2, ...}],
#               size=96, out_path="model.vis")

Measured on examples/cnn_pool.vis (conv 1→4, pool 2, conv 4→2, 96² → 48²):

	numpy	RTX 4050	Intel iGPU	llvmpipe	agreement
conv → pool → conv	10.5 ms	0.8 ms	2.2 ms	1.8 ms	bit-identical

The pool reads a k × k window with x * k + i, so the same inclusive-bound trap applies: two taps for k = 2 means for i in 0 to 1, and a test pins it.

The boundary: inference yes, training not yet#

	Status
inference of a frozen layer (conv, activation, pooling, residual)	works today: the weights enter as canvases from `.npy`, and a 3-D canvas holds the shape a weight tensor needs
multi-layer nets	expressible, one kernel per stage; the channel count is the limit (29 storage buffers on Metal)
training (gradients, optimiser)	no: the language has no autodiff, and the tensor side has autodiff but no convolution layer

So “train an image model here” is false today; “run an image model with bit-identical results anywhere” is true and measured. Training a CNN would need conv + im2col + backward on the tensor side — and there the int8 matrix path (which beats cuBLAS on this card) is what makes it interesting, not the f32 one.

Where the language is underexplored, and it is not insurance#

synthetic training data: a physics simulator that produces images — and produces the same dataset on every machine, which is what stops dataset drift between the machine that trained a model and the one that reproduces it;
domain randomisation with scenarios: thousands of scene variations, each reproducible, for training robustness;
reproducible pre- and post-processing for imaging that ends up in a file — medical, microscopy, satellite, forensics: the chain that produced an image must be re-runnable, not merely documented;
inference attestation: “this mask came from these weights on this image”, checkable without a GPU;
edge deployment: one source; iGPU, discrete GPU or CPU; the same bits.

Library	What it brings	What this engine adds
geopandas (1.2, BSD-3)	vector data, spatial joins, CRS handling, I/O through pyogrio	the arithmetic inside a zone, reproducible and receipted
rasterio (1.5+, BSD-3)	raster I/O, windows, mosaics, reprojection (warp)	the compute on each window, with the same bits everywhere
rasterstats (BSD-3, beta)	zonal statistics and point queries in one call	the same statistics, plus the receipt and cross-device equality
rioxarray / xarray	labelled dimensions and lazy computation	the kernel, not the bookkeeping
shapely / pyproj	geometry and CRS transforms	nothing — this engine does neither

The rule that keeps this honest: a polygon enters as a mask, a raster enters as an array, and the number that leaves carries the receipt. Nothing in the engine depends on any of those libraries; the adapters that ship with it are optional, and one of them (the raster interop) works with geonumpy, rasterio and rioxarray without importing any of them.