Add balanced service area partitioning

[brendancol]

Author of Proposal: @brendan

Reason or problem

allocation assigns each cell to the nearest source by cost distance. That works when proximity is what you care about, but not when workload matters. Five fire stations on a cost surface: pure nearest-source allocation might give one station a tiny cheap patch and leave another covering a huge high-friction zone.

Right now the only way to get balanced territories is to write your own rebalancing loop on top of cost_distance, and it's easy to get wrong.

Proposal

Add a balanced_allocation function that runs cost_distance per source, then iteratively adjusts additive biases until territories have roughly equal cost-weighted area.

Design:

1. Run cost_distance once per source to get N individual cost surfaces.
2. Assign each cell to the source with the lowest biased cost: argmin(cost[i] + bias[i]).
3. Sum the friction values within each territory to get cost-weighted area.
4. Bump biases up for territories that are too large, down for too-small ones (proportional update).
5. Repeat 2-4 until the worst-case deviation from the mean is below tolerance, or max_iterations is hit.
6. Return the allocation as a DataArray of source IDs.

Step 1 is the expensive part, and it's handled by the existing cost_distance function, so all four backends (NumPy, CuPy, Dask+NumPy, Dask+CuPy) work out of the box. The balancing loop is just argmin + masked sums, no custom kernels needed.

Usage:

from xrspatial import balanced_allocation

result = balanced_allocation(
    raster,       # source locations (non-zero = source ID)
    friction,     # cost surface (positive values)
    tolerance=0.05,       # stop when within 5% of the mean
    max_iterations=100,
)

Why bother:

Balanced territory partitioning comes up a lot in logistics, emergency response planning, and delivery routing. Wrapping it in one function call saves people from reimplementing the same bias-adjustment loop.

Stakeholders and impacts

New public function in a new module (xrspatial/balanced_allocation.py). No changes to existing functions. Only depends on cost_distance.

Drawbacks

• Running cost_distance N times (once per source) gets expensive as N grows. With hundreds of sources you'd want a tight max_cost bound to keep it tractable.
• The balancing loop materialises territory weights (scalar sums) each iteration, so Dask can't stay fully lazy through the whole computation.

Alternatives

• Add a weights parameter to the existing allocation function. The iterative rebalancing changes the semantics enough that a separate function makes more sense.
• Weighted Voronoi generators. Doesn't account for friction surfaces at all.

Unresolved questions

• Should the function return just the allocation, or also the per-territory weights and convergence info?
• Is a capacity parameter (per-source max workload) worth adding, or is equal-split good enough for a first pass?

Additional context

This is sometimes called "balanced Voronoi partitioning" or "capacitated territory design" in OR literature. The proportional bias update typically converges in under 20 iterations for well-separated sources.