Compact and connected plans, and zones

The plans of part 4 meet every target at the least cost, and nothing else. Nothing tells the solver that two neighbouring cells are worth more than two cells far apart, so the plans come in 30 pieces, some of them a single square kilometre. In this part you will give plans a shape: compact with a boundary penalty (and calibrate it with a dashboard), free of isolated cells with neighbour constraints, in one piece with contiguity, and connected along a species' habitat with a connectivity penalty. Then you will change the planning units themselves, to H3 hexagons, and let a plan choose between two management zones, strict protection and sustainable use.

The idea#

Why shape matters#

Reserves are not abstract sets of cells. A fragmented plan has more edge, exposed to wind, fire, hunting and invasive species; small pieces hold small populations that wink out; every separate piece needs its own fences, patrols and negotiations; and animals must cross the land between them. Spatial reserve design asks for plans that are compact (little edge for their area) and connected (pieces close together or linked), at the least extra cost (Williams, ReVelle & Levin, 2005). The question is always how much compactness is worth: a plan of one round block is easy to manage and may miss half the features.

Penalties and constraints#

There are two ways to ask for a shape, and they behave differently:

• a penalty adds a price to the objective for what you do not want (boundary, broken links). The solver trades it against the cost: a preference, as strong as its weight. Its value must be calibrated;
• a constraint is a rule the plan must obey (every selected cell has two selected neighbours; the plan is one piece). No value to calibrate, but every rule raises the cost, and some make the problem much harder to solve.

Constraints on shape: neighbours (every selected unit has at least \(k\) selected neighbours) removes isolated cells; contiguity (all selected units form one connected patch, Önal & Briers, 2006) forces a single network; feature contiguity asks it of every feature's units separately.

Connectivity penalties generalise the boundary: a matrix of how strongly each pair of units is linked (by adjacency, by distance, by a dispersal probability, by the habitat between them), and a reward for selecting strongly linked pairs together. Asymmetric versions express links with a direction: a reach of river upstream of another (Beger et al., 2010; Hermoso et al., 2011).

Planning units and zones#

Planning units need not be cells. Properties, watersheds or hexagons are often better: they match how land is managed, and hexagons have equal areas and six equal neighbours. H3 is a global grid of hexagons at sixteen resolutions (resolution 7 cells cover about 5 km², resolution 6 about 36 km²). prioritize.h3_units covers a cost map with them; prioritize.optimize takes any polygons as units and sums every feature over each one exactly.

Management zones let each unit take one of several uses, each with its own cost and its own value for the features (Watts et al., 2009, Marxan with Zones; Hanson et al., 2025). A feature may count fully under strict protection, half under sustainable use; the use zone may cost less, since some income continues. The plan decides, for every unit, which zone (or none).

When to use it#

• Almost every plan on fine units needs a boundary penalty: without it, plans scatter. Calibrate it.
• Neighbour constraints when isolated cells make no sense on the ground but several pieces are acceptable.
• Contiguity for small problems or few features close together, or when one network is a requirement (a corridor system); it is slow on large ones.
• Connectivity penalties when you know how a feature moves: rivers, dispersal distances, least-cost distances (prioritize.decay_connectivity turns any distances into probabilities).
• Polygons or hexagons when decisions are made on properties or management units, or for regional and global analyses (hexagons of equal area).
• Zones for zoning plans: marine spatial plans, land-use plans with protection, sustainable use and restoration.

What you need#

• The features and the cost of the previous parts: serra-verde.zip.
• For hexagons, a cost map in a known coordinate system (Serra Verde's is UTM zone 23 South).
• For zones, a cost map per zone and the zones' contributions to the features.
• The boundary and adjacency of raster cells are computed by prioritize.optimize itself; for other matrices, prioritize.boundary_matrix, prioritize.adjacency_matrix, prioritize.connectivity_matrix, prioritize.proximity_matrix and prioritize.decay_connectivity make them as tables (i, j, x).

Build the model#

A boundary penalty#

format = "terraflux/1"

[[nodes]]
id = "features"
op = "raster.read"
params = { path = "features.tif" }

[[nodes]]
id = "cost"
op = "raster.read"
params = { path = "cost.tif" }

[[nodes]]
id = "plan"
op = "prioritize.optimize"
inputs = { features = "features.raster", cost = "cost.raster" }
params = { objective = "min_set", relative_targets = 0.3, boundary_penalty = 0.01, edge_factor = 0.5, gap = 0.01 }

[[nodes]]
id = "evaluation"
op = "prioritize.evaluate"
inputs = { features = "features.raster", cost = "cost.raster", solution = "plan.solution" }
params = { relative_targets = 0.3 }

[[plots]]
id = "show"
layers = [{ from = "plan.solution" }]

The spatial setting takes several penalties and constraints at once, in order: spatial = [{ kind = "boundary", penalty = 0.01, edge_factor = 0.5 }, { kind = "neighbors", k = 1 }].

Neighbours and contiguity#

Constraints of shape add many variables (contiguity adds a flow through the network) and are best tried on coarser units first. prioritize.coarsen sums the cost and the features over blocks of factor × factor cells, prioritize.refine brings the coarse plan back to the fine grid:

format = "terraflux/1"

[[nodes]]
id = "features"
op = "raster.read"
params = { path = "features.tif" }

[[nodes]]
id = "cost"
op = "raster.read"
params = { path = "cost.tif" }

[[nodes]]
id = "blocks"
op = "prioritize.coarsen"
inputs = { cost = "cost.raster", features = "features.raster" }
params = { factor = 3 }

[[nodes]]
id = "plan"
op = "prioritize.optimize"
inputs = { features = "blocks.features", cost = "blocks.cost" }
params = { objective = "min_set", relative_targets = 0.2, neighbors = 2, gap = 0.01, time_limit = 120 }

[[nodes]]
id = "fine"
op = "prioritize.refine"
inputs = { selection = "plan.solution", like = "cost.raster" }

neighbors = 2 asks for two selected neighbours (of the four sharing an edge) for every selected block; clamp (true) lowers the number for blocks that have fewer neighbours, on the coast. contiguity = true asks for one patch instead (with adjacency_data for polygons).

Connectivity along a species' habitat#

[[nodes]]
id = "tapir"
op = "raster.calc"
inputs = { f = "blocks.features" }
params = { expression = "f.lowland_tapir" }

[[nodes]]
id = "links"
op = "prioritize.connectivity_matrix"
inputs = { units = "blocks.cost", y = "tapir.map" }

[[nodes]]
id = "plan"
op = "prioritize.optimize"
inputs = { features = "blocks.features", cost = "blocks.cost", connectivity_data = "links.table" }
params = { objective = "min_set", relative_targets = 0.2, connectivity_penalty = 0.005, gap = 0.01 }

prioritize.connectivity_matrix gives every pair of adjacent units the length of their shared boundary times the mean of y in the two: here the tapir's habitat suitability, so the links run along its rivers. Its table (i, j, x) goes to the connectivity_data input, and connectivity_penalty weighs it. A positive value rewards selecting linked pairs together, unlike the boundary penalty, which punishes the outline. For one-way links, asym_connectivity_data and asym_connectivity_penalty.

Hexagons#

format = "terraflux/1"

[[nodes]]
id = "features"
op = "raster.read"
params = { path = "features.tif" }

[[nodes]]
id = "cost"
op = "raster.read"
params = { path = "cost.tif" }

[[nodes]]
id = "hexagons"
op = "prioritize.h3_units"
inputs = { cost = "cost.raster" }
params = { resolution = 7 }

[[nodes]]
id = "plan"
op = "prioritize.optimize"
inputs = { polygons = "hexagons.units", features = "features.raster" }
params = { objective = "min_set", relative_targets = 0.3, cost_column = "cost", gap = 0.01 }

[outputs.selection]
from = "plan.selection"

prioritize.h3_units returns the hexagons as polygons (units, with their H3 index and their cost, the costs of the cells under each weighted by the share of each cell inside it) and a table. Connect units to the polygons input of prioritize.optimize and name the cost column (cost_column = "cost"): the features are summed over each hexagon exactly, and the boundaries and adjacency come from the geometries. The plan is the selection table, a row per hexagon in the order of the polygons; write the polygons with vector.write to map it. Any polygons work the same way: properties, watersheds, municipalities read with vector.read.

Zones#

format = "terraflux/1"

[[nodes]]
id = "features"
op = "raster.read"
params = { path = "features.tif" }

[[nodes]]
id = "cost"
op = "raster.read"
params = { path = "cost.tif" }

[[nodes]]
id = "use_cost"
op = "raster.calc"
inputs = { c = "cost.raster" }
params = { expression = "c * 0.45" }

[[nodes]]
id = "zone_costs"
op = "stack.of"
inputs = { protect = "cost.raster", use = "use_cost.map" }

[[nodes]]
id = "plan"
op = "prioritize.optimize"
inputs = { zone_costs = "zone_costs.stack", features = "features.raster" }
params = { objective = "min_set", relative_targets = 0.3, contribution = { use = 0.5 }, gap = 0.01 }

[[plots]]
id = "show"
layers = [{ from = "plan.solution" }]

mandatory_allocation = true puts every unit in some zone (for plans that zone a whole region); boundary_zones says how much an edge between two zones counts in the boundary penalty.

Read the results#

Compactness and its price#

Enlarge

PlayBoundary penalty (thousand USD per metre)

Unlike Marxan's runs in part 3, the cost now only rises with the penalty: the solver finds the cheapest plan for every value. The first steps are almost free: a penalty of 0.01 removes 31 % of the boundary and more than half the patches for 1.2 % more cost. Beyond it, compactness gets dear: halving the boundary costs 7.6 %.

A penalty near the bend, here 0.005 to 0.01, is a defensible default. Part 6 finds it reproducibly with prioritize.calibrate (the knee of the curve, or TOPSIS). Report the value, how it was chosen, and the cost and boundary apart.

Rules on the shape#

Enlarge

Contiguity is the most demanding rule, in money and in time: the network must thread from the mountains to the coast through dear land, and the solver cannot prove its plan within 23 % of the best in three minutes. Use it on small problems, or when the rule is a requirement; elsewhere a boundary penalty gets most of the compactness for a fraction of the cost. Report the gap: a contiguous plan "within 23 % of the optimum" is still a valid plan, but not a proven cheapest one.

Connectivity along the rivers#

The connectivity of a plan, connectivity in the summary of prioritize.evaluate given the same matrix on its connectivity_data input, is the sum of the links between selected units. Like the boundary penalty, the value of the connectivity penalty is an exchange rate between money and links, and needs the same calibration.

Hexagons#

Larger units cost more for the same targets: a unit is bought whole, with the poor land inside it. They also give plans that are compact by construction. Choose the unit by the decision: cells for analysis, hexagons or properties for plans that will be implemented, and report the units with the results.

Zones#

Enlarge

Enlarge

Without the 20 % of strict protection, the plan has 452 strictly protected cells and 735 of sustainable use, for 36.2 million USD: cheaper per cell, sustainable use covers more land. The strict targets cost little here (0.4 %), but they move 81 cells into strict protection. Zoning plans make visible a choice that a plan of protected areas hides: how much of the targets may be met by land that is used.

Settings that change the answer#

• The penalty's value, as the slider shows; it has units and depends on the landscape: calibrate it, never copy it.
• edge_factor: 0 when the edge of the study area is not a real edge (a political border through continuous forest), 1 when it is (a coast that the plan should not hug); the default 0.5 is in between.
• The size of the units, for constraints (blocks of 3 km here) and for the plan itself (hexagons).
• The contributions of zones, which are assumptions about how much a use harms each feature: take them from studies or experts, and test other values.

Common pitfalls#

• Reporting the objective as the cost. With penalties, the objective includes them; read the cost from prioritize.evaluate.
• A penalty copied from another study, or chosen for the map it gives. Calibrate on the curve.
• Contiguity on a large problem without a time limit: it may run for hours. Set time_limit, read the gap, and try coarser units.
• Zone layers named differently in zone_costs, contribution and manual_targets: the names tie them together.
• Comparing boundaries computed with different edge factors: Marxan's maps use 1 by default, Terraflux's 0.5.

Going further#

• Part 6, Evaluating, comparing and reporting plans: calibration with the knee and TOPSIS, trade-off curves, and a checklist.
• Part 7, Multi-action planning: actions against threats, with the boundary of each action.
• Connectivity from least-cost or resistance distances: Distance and restoration and prioritize.decay_connectivity.
• The operators: prioritize.optimize, prioritize.h3_units, prioritize.connectivity_matrix, prioritize.coarsen, prioritize.refine.

References#

• Beger, M., Linke, S., Watts, M., Game, E., Treml, E., Ball, I., & Possingham, H. P. (2010). Incorporating asymmetric connectivity into spatial decision making for conservation. Conservation Letters, 3(5), 359–368. https://doi.org/10.1111/j.1755-263X.2010.00123.x
• Beyer, H. L., Dujardin, Y., Watts, M. E., & Possingham, H. P. (2016). Solving conservation planning problems with integer linear programming. Ecological Modelling, 328, 14–22. https://doi.org/10.1016/j.ecolmodel.2016.02.005
• Hanson, J. O., Schuster, R., Strimas-Mackey, M., Morrell, N., Edwards, B. P. M., Arcese, P., Bennett, J. R., & Possingham, H. P. (2025). Systematic conservation prioritization with the prioritizr R package. Conservation Biology, 39(1), e14376. https://doi.org/10.1111/cobi.14376
• Hermoso, V., Linke, S., Prenda, J., & Possingham, H. P. (2011). Addressing longitudinal connectivity in the systematic conservation planning of fresh waters. Freshwater Biology, 56(1), 57–70. https://doi.org/10.1111/j.1365-2427.2009.02390.x
• Önal, H., & Briers, R. A. (2006). Optimal selection of a connected reserve network. Operations Research, 54(2), 379–388. https://doi.org/10.1287/opre.1060.0272
• Watts, M. E., Ball, I. R., Stewart, R. S., Klein, C. J., Wilson, K., Steinback, C., et al. (2009). Marxan with Zones: software for optimal conservation based land- and sea-use zoning. Environmental Modelling & Software, 24(12), 1513–1521. https://doi.org/10.1016/j.envsoft.2009.06.005
• Williams, J. C., ReVelle, C. S., & Levin, S. A. (2005). Spatial attributes and reserve design models: a review. Environmental Modeling & Assessment, 10(3), 163–181. https://doi.org/10.1007/s10666-005-9007-5