Five changes account for most of the speed-up.
Plotline's decision layer re-checks a registered plot 44× faster than the first version, and a new boundary 5.2× faster. A few changes produced those numbers; they are described here in the order they were made.
- Five changes: keep only what the decision reads; find a registered plot's cells once; decode on the packed bit-planes; scan small polygons in plain Python; drop NumPy for re-checks. A leaner shared reducer helped all of them.
- Two separate scoreboards. Registered plots: 24× at the start, 44× now (83× without NumPy). New boundaries: 2.2× at the start, 5.2× now. All measured hot, against the retained original path (P0). The two are never combined into one multiple.
- No phase changed an answer, a stored byte or a file format. The same proof, 744,732 of 744,732 cells, was re-run after every phase.
Each multiple depends on its baseline and its clock
Every multiple here is measured against the original query path from version 1, retained as the baseline and called P0. It reads three Hansen bands and builds a 101-bin canopy histogram. P0 is re-measured in every benchmark run, so a multiple can drift about half a percent without any code change.
Each multiple also belongs to one of three clocks, and multiples from different clocks are not comparable.
The layer is already loaded and the code has warmed up. One pass over the 17 supplied plots is timed and repeated hundreds of times. This is the per-plot cost once the layer is local.
A new Python process loads the layer and answers the 17 plots. Imports are done first; start-up is not counted. A smaller layer shows up here.
Python is launched from scratch and prints the answer. The time includes interpreter start-up (about 14.6 ms) and importing NumPy (about 25 ms). A do-nothing child in the harness already costs 20.1 ms.
Path codes
Each path is a different way of answering the same question. The suffix d means the decision level (s is the slightly larger screening level).
| Code | What it is | Answers |
|---|---|---|
| P0 | The retained v1 query on the original layers. Baseline for every multiple. | any plot |
| P1 | The original layers, asked the same EUDR question as the decision layer. The like-for-like baseline for the representation step. | any plot |
| P3 | The decision store with the general polygon selection. | any plot |
| P4 | The decision store with exact fast selection: the pure-Python scan and a fused reduction. | new boundaries, one at a time |
| P4B | The decision store with the batch evaluator: many new polygons in one pass. | new boundaries, many per call |
| P5d | The decision store plus a registered-plot index (cells found once). | registered plots only |
| P5Ld | The registered-plot re-check with no NumPy: standard library only. | registered plots only |
| P6d | The registered-plot answer table: stored outputs, looked up by key. | registered plots only |
Phases
The decision-layer work was followed by four short rounds of speed-up.
Keep only what the decision reads. Registered-plot index and answer table. Benchmarked once, as the starting point.
Same code, made faster: decode on packed planes, a quicker reducer, a quicker polygon selection.
Shaped to the workload: scan small polygons in plain Python, a batch evaluator, a NumPy-free re-check.
Two small wins: a reducer twice as lean, and no hidden import in the batch path.
A trusted-integer reducer entry, a smarter rule for when to scan in Python, and a research pass.
Registered plots went 24× to 44×. New boundaries went 2.2× to 5.2×.
Each bar is a multiple against P0, measured hot on the 17 plots. The two groups have separate scales: a registered-plot re-check and a new-boundary check are different jobs and are never quoted as one number.
Decision layers · starting pointA small decision layer and a registered-plot index.
The answer table, about 420× at the start and 410× now, is left off this chart because it is a cache (see breakthrough 2).
All rows by phase
The same data as a table, one column per phase. A green bar and a bold value mark a cell that improved on the previous column by more than about 3%; everything else is flat within noise. Lower is better for rows marked ↓.
Multiples are against P0 in the same run. The answer table reads 410× in Phase 4, not 417×, because the table time (0.206 µs per plot) did not change within noise and P0 re-measured 0.5% lower. It is not a slowdown. New boundaries read 5.2× rather than 5.3× for the same reason.
What each of the five changes contributes
The changes are listed in the order they were made; each answers a question the previous one leaves open. The cards show which clock each one moves. The numbers are step-on-step multiples, so they do not multiply together.
1.6× hot · 7.0× fresh 2 · IndexFind a registered plot's cells once24× → 44× hot
answer table ~410× (a cache) 3 · DecodeDecode on the packed bit-planesstore read 2.5 → 0.67 ms
peak memory 43 → 37 MB 4 · SelectionScan small polygons in plain Pythonnew boundaries 2.2× → 5.2×
batch of 500 ~15.8× 5 · No NumPyDrop NumPy for re-checks; lean reducer44× → 83× hot
end-to-end 50.0 → 24.9 ms
| Step (decision level) | Hot, 17 plots | Fresh process | End-to-end |
|---|
What each step is worth, measured one step at a time on the same machine and questions. The hot and end-to-end columns are rarely both large.
Keep only what the decision reads.
The EUDR screen asks whether a cell was forest on 31 December 2020 and whether that forest has been lost since. Nothing in that question reads the 2000 canopy percentage, which is 72% of Hansen's bytes, or the loss years before the cut-off. The original query still loaded three Hansen bands and built a 101-bin canopy histogram on every call.
The decision layer drops all of that and keeps one small band: a 20-value code per cell. The original layers are 623,237 B; the decision store is 6,146 B, 99.0% smaller.
Bytes the query has to hold
FigureFour layers, plus a canopy histogram that no question reads.
The shapes drawn are illustrative, not the real rasters.
Go deeper: what the 20 codes are technical
Each Hansen cell becomes one of 20 codes: unknown, water, or land crossed with three forest states (not forest, forest, forest unknown in 2020) and six loss states (none, 2021 to 2025). The decision level additionally forgets the loss year on land that was not forest in 2020. One shared reducer turns a 20-bin histogram into the contract answer.
Find a registered plot's cells once, then reuse them.
EUDR plots are registered once and re-checked every time a layer updates. Which cells sit inside a plot depends on the plot's boundary and the grid, not on the layer's values. So the expensive part, finding the cells, can be done once and saved as a short list of cell numbers. A re-check then counts the codes at those cells and turns the count into the answer.
That took the re-check from 85.4 µs per plot to 3.55 µs at the start (24×), and to 1.90 µs now (44×). The index gave the initial gain; a leaner reducer (breakthrough 5) improved it further.
Time spent on 5 re-checks of one plot per-plot hot times, to scale
FigureScan the polygon once.
Registering a plot costs one selection, once; that is not timed separately, so it is not in the bars. The grid and polygon are invented.
Go deeper: how the index and table work, and what hashing costs technical
The registry stores each plot's in-chip flat cell indices (row × width + column, int32, sorted). Evaluating a plot is a bincount of the codes at those indices, passed to the shared reducer. The answer table returns the stored output only on an exact plot, geometry, layer and level match, otherwise nothing.
The realistic workflow stores each plot's geometry hash at registration. If a service instead re-hashes the polygon on every call, cost rises: P5d 1.90 to 5.97 µs, P6d 0.21 to 3.65 µs per plot. The 44× and 410× figures use the precomputed key.
Decode on the packed bit-planes.
The decision store keeps the map as a few yes/no planes (is it unknown? water? forest?) packed eight cells to a byte, plus a short list of loss-year events. The first reader unpacked every plane into 744,732 separate cells before it checked or used any of them.
Phase 1 does its checks on the packed bytes (93 KB), places the sparse loss years only where an event is, and finds ranked event cells by counting set bits. The file bytes and format did not change. Reading the decision store went from 2.5 ms to 0.67 ms (3.8×).
Reading the decision store
FigureBefore: expand every plane to one cell per byte.
The picture is schematic: the real chip is 1,233 × 604 cells.
Go deeper: the three Phase 1 changes technical
The store decoder runs its plane checks on strip-order packed bytes, and only the dense forest plane is expanded. A side effect: the validator's crafted-blob memory case now peaks near 6 MB instead of about 50 MB. The other two changes: the histogram reduction (1.6 to 2.2× on its own) and polygon selection, which now builds the polygon straight from one coordinate array and paints ambiguous-cell rectangles without np.add.at. It calls a few low-level Shapely functions directly, so it depends on the pinned shapely 2.1.2.
For a small polygon, a plain Python loop beats NumPy.
For a new boundary there is no saved list of cells, so they must be found. The rule is the centre rule: a cell counts if its centre is inside the plot. A scanline does it: for each map row, find where the boundary crosses the row; the centres between crossings are inside.
A plot is only about 120 cells. NumPy pays a fixed cost on every call, and at that size the cost exceeds the work. Phase 2 therefore scans small windows in plain Python and builds the histogram straight from the runs of selected centres. Ambiguous centres, such as an edge exactly through a centre, still go to exact GEOS geometry tests, so the result stays exact.
New boundaries went 3.0× to 5.0× in Phase 2 (2.2× at the start, 5.2× now).
Cost per new polygon 500 synthetic polygons, hot
FigureSweep one row at a time.
Synthetic shapes of 0.5 to 20 ha; the batch figure is the batch time divided by 500.
np.unique imports a module on first use (about 2 ms); a small sort-and-dedupe helper fixed that.Go deeper: when the Python scan is used, and the Phase 4 tuning technical
Windows up to 500 rows × edges (and 4,096 cells) always use the Python scan. Phase 4 also lets it take 500 to 2,000, unless the boundary changes rows more than 150 times in total (a "sawtooth"), where the NumPy scan is faster. The reason: a 224-shape sweep found the Python scan costs about 0.17 µs per row an edge crosses, and rows × edges predicts its time poorly (R² of 0.66, against 0.99 for a fit on edges and row crossings). The two scans were already proven identical. On the 17 plots every polygon was already under the first threshold, so the Phase 4 rule change did not touch them.
Drop NumPy for re-checks, and make the shared reducer lean.
This is the only change that moves the whole-process time. A fresh Python process spends about 25 ms just importing NumPy. A registered re-check only needs the codes at a few hundred cells, so it never needed NumPy.
The NumPy-free re-check reads and validates the store using only the standard library, then decodes codes only at the registered cells, straight from the packed planes. It never builds the 744,732-cell raster. End-to-end time for a registered re-check halves, 50.0 to 24.9 ms, and peak memory drops from 37 to 23 MB.
FigureRegistered index: half the bar is the harness floor, a third is imports.
The third segment is computed by subtraction, not measured separately.
A reducer twice as lean
Every live answer ends in one function, the reducer, that turns the 20-bin histogram into the contract answer. In Phase 3 it unpacks the 20 numbers by name instead of slicing, zipping and calling a nested helper: about 1.6 to 0.8 µs per call. That alone lifted the registered index from 33× to 44× and the NumPy-free re-check from 44× to 65×. In Phase 4 a trusted-integer entry skips a list conversion that callers who already hold ints do not need: the NumPy-free re-check goes from 65× to 83× (1.31 to 1.02 µs per plot).
The reducer alone per call
# Phase 3: unpack the 20 codes by name
_, water, n0, n1, n2, n3, n4, n5, f0, ... = counts
Per plot, by phase µs, lower is better
The NumPy-free re-check begins in Phase 2.
No single stage dominates a new-boundary check
A new boundary now costs about 16 µs. An independent audit of the 17 plots measured each stage. The largest, scanning the rows, is 30% of the total. Every other stage is smaller, so no single remaining change would double the speed.
The audit total (16.35 µs) is a separate measurement from the benchmark's 16.15 µs; they differ by about 1%. It comes from the 17 plots only. The validity check grows with vertex count: about 0.84 µs at 8 vertices, 1.25 at 16, 4.1 at 60 and 12 at 150.
The same paths give different multiples on the three clocks
The bars show the same seven paths on the same 17 plots under each clock. Hot is where the index and the table dominate. Fresh-process time rewards the small layer. End-to-end barely moves, because 20 ms of it is the harness and a fresh Python and NumPy import is most of the rest.
Registered-only paths (P5d, P5Ld, P6d) are shown in their own group so they are not read as a speed-up for a new plot.
Hot
Nothing here is disk-cold: the OS file cache was warm.
No better option was found for small exact plots on a CPU
Phase 4 included a review of papers, forums and engineering options. The centre rule is the standard zonal-statistics default (GDAL, rasterstats). The closest research, Raster Join, does what Plotline does in its exact form (a fast interior pass and exact checks at the border); its bounded form gives up exactness, which this use case cannot accept. Other approaches do not fit plots of about 120 cells.
| Option | Finding | Verdict |
|---|---|---|
| GPU for queries | A plot is about 120 cells; GPU launch and transfer is around 10 µs, about the cost of the whole query. Apple GPUs have no float64, so the exact centre maths would change. GEOS runs on the CPU. | skipped |
| GPU for building the layer | Building is 92% compression, which tries 10 settings per stream and mostly picks LZMA (no GPU version). Free code derivation would make the build only about 1.06× faster. | skipped |
| Python 3.14 JIT | Not in the Homebrew build; where available, 2 to 3% slower. | skipped |
| Tail-calling interpreter | About 2 to 5%, and needs a different interpreter. | skipped |
| Cython on the module | About 6%, and needs per-platform builds. | skipped |
| mypyc · numba · PyPy | mypyc failed to import; numba takes 0.2 to 0.7 s to compile on the first call; PyPy has no shapely or numpy wheels. | skipped |
| GEOS 3.14 / shapely 2.2 | About 3%. | skipped |
| Hand-written validity check | About 8× slower than GEOS in Python, with risky edge cases. | skipped |
| Change-only re-checks | When a layer updates, re-check only the plots whose cells changed. Needs an inverted cell-to-plot index (incremental view maintenance). | roadmap, not built |
| Direct GEOS binding via ctypes | A prototype cut a fresh new-plot process from about 48 to about 20 ms and saved about 2.7 µs per plot, with 0 mismatches on 60,000 degenerate polygons. 3 to 5 hours of careful work. | roadmap, not built |
FigureCompression is 92% of building a layer. A GPU would speed up the 8% that is not compression, so the build would be about 1.06× faster.
A scratch measurement on 6 of the 160 strips of one coastal tile; the estimate of about 7 CPU-minutes per tile and roughly 30 to 60 CPU-hours for about 280 land tiles is rough, since inland tiles encode more slowly.
Answers, stored bytes and file format are identical in every phase
No phase changed any output, any stored byte or any file format. After each one the benchmark was re-run with every timed answer bound to the same expected digests, and the proof was re-run. The table lists what each phase was compared against in addition.
| Phase | Cells identical | Fuzzed polygons | What the new code was compared with |
|---|---|---|---|
| Phase 1 | ✓ 744,732 / 744,732 | ✓ 9,829 / 9,829 | About 16,000 corrupted store blobs (identical decodes and rejections); 180,000 reducer comparisons; about 9,000 polygons against the old fast path and select_cells. 300 engine and 311 app tests; all 17 validators exit 0. |
| Phase 2 | ✓ 744,732 / 744,732 | ✓ 9,829 / 9,829 | 58,896 polygon and grid cases against the old scan and select_cells; about 185,000 batch cases against the single path; about 118,000 corrupted blobs for the NumPy-free path, all rejected. |
| Phase 3 | ✓ 744,732 / 744,732 | ✓ 9,829 / 9,829 | 300,000 random histograms against the old reducer (outputs, key order, errors); a permanent test on 40,000 seeded cases. |
| Phase 4 | ✓ 744,732 / 744,732 | ✓ 9,829 / 9,829 | 14,500 polygons against select_cells; 20,000 cases for the new reducer entry. 361 engine tests; all 17 validators exit 0. |
"Identical" means every output field under the documented centre rule, not accuracy against the real world.
Re-checking 10 million new plots costs about three-tenths of a cent.
The figures are Lambda compute for re-checking 10 million plots over the tropical belt, at assumed list prices and assumed tile counts.
| New-plot line, by phase | Start | Phase 1 | Phase 2 | Phase 3 | Phase 4 |
|---|---|---|---|---|---|
| Lambda compute, 10M new plots | $0.0069 | $0.0053 | $0.0032 | $0.0031 | $0.0031 |
Compute time: original 1,054 s, registry 42 s, new polygons 185 s. The new-plot line prices P4 on the 17 plots, which did not change in Phase 4. Storage for the tropical belt (measured basis) is $0.018 a month for the decision layer against $0.90 for the original layers. Not priced: cold object-store latency, network time and the 15-minute Lambda limit. The per-strip container this assumes is not implemented.
Design limits
- Edge cells are all-or-nothing. Under the centre rule a cell counts for a plot when its centre lies inside the boundary; cells are about 0.076 ha. Partial edge coverage and the 0.5 ha minimum area are not modelled.
- JRC is sampled once per Hansen cell, at the centre of its 3 × 3 block. The alternatives (majority, any forest) give the same decision for all 17 plots, but only one plot has mapped forest.
- The layer is rule-specific. A changed rule or an added layer means re-deriving from the lossless audit tier: seconds per region, about 7 CPU-minutes per tile (estimated).
- The largest speed-ups need registered plots. The 44× to 410× figures apply to plots registered in advance; a new boundary gets 5.2×.
- Storage was measured on a whole tile and query timing on one region.
The speed story in one paragraph
A short account of how the speed was reached.