def run_raster_recipe(
recipe_id: str,
inputs: Mapping[str, Path],
output_path: Path,
*,
input_units: Mapping[str, str] | None = None,
parameters: Mapping[str, Any] | None = None,
plan: str | None = None,
profile: str | None = None,
scratch_dir: Path | None = None,
block_size: int = 512,
hash_assets: bool = False,
overwrite: bool = False,
) -> RasterRecipeResult:
"""Execute an allowlisted recipe over aligned numeric rasters.
Inputs must already share a CRS, transform, dimensions, and pixel grid. This
function intentionally does not reproject or interpolate hydraulic results.
Source surfaces should first be generated by RASMapper/RasProcess.
"""
import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.shutil import copy as copy_raster
from rasterio.windows import Window
recipe = get_raster_recipe(recipe_id)
supplied_roles = set(inputs)
required_roles = set(recipe.input_roles)
if supplied_roles != required_roles:
raise ValueError(
f"Recipe {recipe_id} requires inputs {sorted(required_roles)}; "
f"received {sorted(supplied_roles)}"
)
if block_size < 64 or block_size > 4096:
raise ValueError("block_size must be between 64 and 4096 pixels")
if recipe.requires_synchronized_profile:
normalized_profile = str(profile or "").strip().lower()
if not normalized_profile or normalized_profile in {"max", "maximum", "min", "minimum"}:
raise ValueError(
f"Recipe {recipe_id} requires one synchronized timestep/profile; "
"independent Max/Min rasters are not valid inputs"
)
output_path = Path(output_path)
if output_path.suffix.lower() not in {".tif", ".tiff"}:
raise ValueError("Raster recipe output must be a .tif or .tiff")
if output_path.exists() and not overwrite:
raise FileExistsError(f"Raster recipe output already exists: {output_path}")
output_path.parent.mkdir(parents=True, exist_ok=True)
input_paths = {role: Path(inputs[role]).resolve() for role in recipe.input_roles}
missing = [str(path) for path in input_paths.values() if not path.is_file()]
if missing:
raise FileNotFoundError(f"Raster recipe inputs do not exist: {', '.join(missing)}")
sources = {role: rasterio.open(path) for role, path in input_paths.items()}
try:
reference = sources[recipe.input_roles[0]]
_validate_grids(sources, reference)
units = _resolve_units(recipe, sources, input_units or {})
effective_parameters = dict(recipe.parameter_defaults)
effective_parameters.update(parameters or {})
_validate_parameters(recipe, effective_parameters)
output_units = _output_units(recipe, units)
scratch_parent = Path(scratch_dir).resolve() if scratch_dir else output_path.parent
scratch_parent.mkdir(parents=True, exist_ok=True)
with tempfile.TemporaryDirectory(prefix="ras2cng-recipe-", dir=scratch_parent) as temp_name:
temporary = Path(temp_name) / "calculated.tif"
profile_options = reference.profile.copy()
profile_options.update(
driver="GTiff",
count=1,
dtype=recipe.output_dtype,
nodata=recipe.output_nodata,
tiled=True,
blockxsize=_valid_tiff_block(block_size),
blockysize=_valid_tiff_block(block_size),
compress="ZSTD",
predictor=2 if recipe.categorical else 3,
BIGTIFF="IF_SAFER",
)
aggregate = _StatisticsAccumulator()
with rasterio.open(temporary, "w", **profile_options) as destination:
destination.update_tags(
recipe_id=recipe.recipe_id,
recipe_version=recipe.version,
units=output_units,
interpolation_authority="RASMapper/RasProcess source rasters",
)
destination.update_tags(1, units=output_units)
for row in range(0, reference.height, block_size):
height = min(block_size, reference.height - row)
for col in range(0, reference.width, block_size):
width = min(block_size, reference.width - col)
window = Window(col, row, width, height)
arrays = {
role: source.read(1, window=window, masked=True, out_dtype="float64")
for role, source in sources.items()
}
values, valid = _calculate_window(
recipe,
arrays,
units,
effective_parameters,
)
aggregate.update(values[valid])
filled = np.full(values.shape, recipe.output_nodata, dtype=recipe.output_dtype)
filled[valid] = values[valid].astype(recipe.output_dtype, copy=False)
destination.write(filled, 1, window=window)
factors = _overview_factors(reference.width, reference.height)
if factors:
overview_resampling = Resampling.nearest if recipe.categorical else Resampling.average
destination.build_overviews(factors, overview_resampling)
destination.update_tags(
ns="rio_overview",
resampling=overview_resampling.name,
)
staged_output = output_path.with_name(
f".{output_path.name}.{uuid.uuid4().hex}.tmp"
)
try:
copy_raster(
temporary,
staged_output,
driver="COG",
compress="ZSTD",
blocksize=_valid_tiff_block(block_size),
overview_resampling="nearest" if recipe.categorical else "average",
BIGTIFF="IF_SAFER",
)
staged_output.replace(output_path)
finally:
staged_output.unlink(missing_ok=True)
statistics = aggregate.to_dict()
provenance_path = output_path.with_suffix(".provenance.json")
provenance = {
"schema": "ras2cng.raster-recipe/v1",
"recipe": asdict(recipe),
"generatedAt": datetime.now(timezone.utc).isoformat(),
"authority": "ras2cng controlled raster arithmetic",
"interpolationAuthority": "RASMapper/RasProcess source rasters",
"plan": plan,
"profile": profile,
"inputs": {
role: {
"file": path.name,
"sizeBytes": path.stat().st_size,
"modifiedNs": path.stat().st_mtime_ns,
"units": units[role][2],
**({"sha256": _sha256(path)} if hash_assets else {}),
}
for role, path in input_paths.items()
},
"parameters": effective_parameters,
"output": {
"file": output_path.name,
"sizeBytes": output_path.stat().st_size,
"modifiedNs": output_path.stat().st_mtime_ns,
**({"sha256": _sha256(output_path)} if hash_assets else {}),
"units": output_units,
"dtype": recipe.output_dtype,
"nodata": recipe.output_nodata,
"statistics": statistics,
"crs": reference.crs.to_string(),
"transform": list(reference.transform)[:6],
"width": reference.width,
"height": reference.height,
},
}
provenance_path.write_text(json.dumps(provenance, indent=2) + "\n", encoding="utf-8")
return RasterRecipeResult(
output_path=output_path,
provenance_path=provenance_path,
recipe_id=recipe.recipe_id,
recipe_version=recipe.version,
units=output_units,
statistics=statistics,
)
finally:
for source in sources.values():
source.close()