Estimating a retrieval radius without retrieving

A seven-step geometric binary search search sizes a spatial radius query in 1 to 3 seconds, against 8 to 25 seconds for running the pipeline to find out.

Concentric rings expanding from a single point over dark terrain contours, one ring highlighted in amber

A radius query needs a radius. The obvious way to pick one is to try a radius, run the pipeline, see whether it produced enough supply, and adjust. That works and it takes 8 to 25 seconds per attempt. Doing it inside a conversation where the user is waiting is not viable.

The problem

The retrieval step pulls every forest cluster within some distance of a facility and scores it. The radius has to be large enough that the clusters inside it can meet annual demand, and no larger, because every extra kilometre adds clusters that will be scored and then discarded.

Get it too small and the pipeline runs, finds insufficient supply, and has to run again. Get it too large and you have retrieved and scored tens of thousands of clusters you did not need. Both cost a round trip the user can feel.

The circularity is the awkward part. You cannot know the right radius until you know what is inside it, and finding out what is inside it is the expensive operation you are trying to size.

What the estimator does

Break the circularity by answering a cheaper question first. The full pipeline asks what is the delivered cost of every cluster within r. The estimator asks how much biomass is within r, which is a single aggregate over an indexed column and does not require scoring anything.

Define the target radius as the smallest radius whose total biomass meets demand with a buffer:

r* = min { r : sum of biomass within r >= 1.5 * annual_demand }

The 1.5 buffer exists because not every cluster inside the radius will be usable. Some are infeasible for every harvesting system, some are too expensive to make the cut. Sizing to exactly demand guarantees a shortfall.

Then search geometrically rather than linearly:

def estimate_radius(facility, demand, buffer=1.5,
                    r0=10, gamma=1.5, r_max=150):
    target = buffer * demand
    r = r0
    while r <= r_max:
        available = sum_biomass_within(facility, r)
        if available >= target:
            return r
        r = min(r * gamma, r_max)
    return r_max  # flag insufficient supply

Starting at 10 km with a growth factor of 1.5 reaches the 150 km ceiling in seven steps. Seven queries against a pre-aggregated table, 1 to 3 seconds total, against 8 to 25 seconds for one attempt at the full pipeline.

Why the obvious alternatives are worse

A fixed radius. Pick 50 km and skip the estimation. This was measured as an ablation and it fails in both directions. At dense locations like Quincy and Mount Shasta, where the true radius is around 24 km, a fixed 50 km over-retrieves by two to four times, so every downstream computation runs on a cluster set several times larger than needed. At Fresno, where the true radius is 56 km, it under-retrieves, the supply curve comes up short, and the agent has to re-call with a larger buffer, adding a full round trip. Overall completion drops from 96.4% to 94.5% and supply failures triple from one in 55 to three in 55.

A true binary search. Bisection needs a known upper bound that is guaranteed to satisfy the constraint, and here the only such bound is the 150 km ceiling, which is far from the answer in most cases. Bisecting from 10 to 150 wastes its early steps in territory the geometric search skips. Growth by a constant factor spends its steps near the answer instead.

Estimating from density. Divide demand by regional biomass per square kilometre and solve for the radius. Fast, and wrong wherever biomass is unevenly distributed, which is everywhere with terrain in it. The estimator queries actual sums rather than assuming uniformity.

The index trap

One detail matters more than the algorithm.

PostGIS will happily accept this:

SELECT SUM(biomass_bdt) FROM mv_cluster_supply
WHERE ST_DWithin(geom::geography, :point::geography, :radius_metres);

It returns the correct answer. It also runs slowly, because the cast to geography produces an expression the GiST index on geom cannot serve, so the planner falls back to a sequential scan over 11.9 million rows. Inside a seven-iteration loop, that is seven sequential scans.

The fix is to stay in the indexed type and convert the units instead:

SELECT SUM(biomass_bdt) FROM mv_cluster_supply
WHERE ST_DWithin(geom, ST_MakePoint(:lon, :lat), :radius_km / 111.32);

One degree of latitude is about 111.32 km. Dividing gives a radius in degrees, which ST_DWithin can answer from the index directly.

Two query plans side by side, one showing a sequential scan over the full table and the other showing an index scan on the GiST index, with execution times
Same result, same predicate. The cast is the only difference between an index scan and a full table scan.

The approximation is not free. Degrees of longitude shrink toward the poles, so a degree-radius circle is really an ellipse, and it is slightly wrong in the east-west direction. At California's latitudes the error is under 0.3%, which on a 40 km radius is about 120 metres. The estimator is sizing a search boundary that already carries a 50% buffer, so 120 metres does not change any decision. Applying the same trick to a final distance calculation would be a mistake.

This is the general shape of the problem, not a PostGIS quirk. Any predicate that wraps an indexed column in a function or a cast is likely to be unindexable. Check the query plan rather than trusting that a spatial function uses the spatial index.

What it costs

The estimator is an approximation and it is wrong sometimes. Its buffer of 1.5 is a guess about what fraction of clusters will survive scoring, and in terrain where a large share of clusters are infeasible for ground-based equipment, that guess is too low. When the supply curve comes up short, the agent re-calls with a buffer of 2.0 and retrieval runs again.

That fallback fired on 8% of benchmark queries, concentrated in High Sierra locations where ground-based infeasibility exceeds 40%. Those queries pay the extra round trip, so the estimator costs them time rather than saving it. The other 92% get their radius in 1 to 3 seconds.

An obvious improvement is to make the buffer terrain-dependent rather than constant, since the ground-accessible fraction is already stored per hex. That has not been built.

More on how procurement radius is reported at /docs/understanding-supply.

References

  1. PostGIS Project Steering Committee. ST_DWithin. PostGIS documentation.

  2. PostgreSQL Global Development Group. GiST Indexes. PostgreSQL documentation.