#!/usr/bin/env python3 """Terrain-grid and movement-cost functions used by the complete-matrix analysis. The complete-matrix script imports this module to calculate all 22 base paths and 54 principal sensitivity runs. This module also retains the command-line routine for the 14 base and 42 sensitivity comparison baselines. Those stored baselines verify the common grid, costs and display geometry; their display-line lengths are not the analytical lengths reported in the paper. Dijkstra's algorithm finds a minimum-cost path on the constructed raster graph. The base grid uses zoom 8 / downsample 4. Sensitivity configurations vary grid resolution, movement neighbourhood, search margin and non-Tobler terrain terms. Source-endpoint placement sensitivity is not tested. """ from __future__ import annotations import argparse import heapq import json import math import platform import urllib.request from dataclasses import dataclass from pathlib import Path from typing import Iterable import numpy as np from PIL import Image TERRARIUM_URL = "https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png" USER_AGENT = "Obsidian-least-cost-paths/4.5 (archaeological reproducibility)" MARGIN_DEG = 0.9 PRIMARY_ZOOM = 8 PRIMARY_DOWNSAMPLE = 4 SENSITIVITY_KEYS = ( "harb-central", "harb-bingol", "harb-nemrut", "karahan-central", "karahan-bingol", "gobekli-central", "gobekli-bingol", ) SENSITIVITY_VARIANTS = ( { "key": "higher_resolution", "zoom": 8, "downsample": 3, "margin_deg": 0.9, "neighbourhood_radius": 2, "terrain_weight_scale": 1.0, }, { "key": "three_by_three_neighbourhood", "zoom": 8, "downsample": 4, "margin_deg": 0.9, "neighbourhood_radius": 1, "terrain_weight_scale": 1.0, }, { "key": "narrower_search_margin", "zoom": 8, "downsample": 4, "margin_deg": 0.45, "neighbourhood_radius": 2, "terrain_weight_scale": 1.0, }, { "key": "wider_search_margin", "zoom": 8, "downsample": 4, "margin_deg": 1.35, "neighbourhood_radius": 2, "terrain_weight_scale": 1.0, }, { "key": "weaker_non_tobler_terms", "zoom": 8, "downsample": 4, "margin_deg": 0.9, "neighbourhood_radius": 2, "terrain_weight_scale": 0.75, }, { "key": "stronger_non_tobler_terms", "zoom": 8, "downsample": 4, "margin_deg": 0.9, "neighbourhood_radius": 2, "terrain_weight_scale": 1.25, }, ) @dataclass(frozen=True) class Point: key: str name: str lat: float lon: float @dataclass class CostGrid: elev: np.ndarray lats: np.ndarray lons: np.ndarray east_west: np.ndarray north_south: np.ndarray ruggedness: np.ndarray topographic_position: np.ndarray flow_accumulation: np.ndarray tile_ids: list[str] def deg2num(lat_deg: float, lon_deg: float, zoom: int) -> tuple[int, int]: lat_rad = math.radians(lat_deg) n = 2**zoom x = int((lon_deg + 180.0) / 360.0 * n) y = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n) return x, y def pixel_to_lon(global_x: float, zoom: int) -> float: return global_x / (256 * (2**zoom)) * 360.0 - 180.0 def pixel_to_lat(global_y: float, zoom: int) -> float: n = math.pi - 2.0 * math.pi * global_y / (256 * (2**zoom)) return math.degrees(math.atan(math.sinh(n))) def haversine_km(lat1: float, lon1: float, lat2: float, lon2: float) -> float: radius_km = 6371.0 dlat = math.radians(lat2 - lat1) dlon = math.radians(lon2 - lon1) a = ( math.sin(dlat / 2.0) ** 2 + math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) * math.sin(dlon / 2.0) ** 2 ) return 2.0 * radius_km * math.atan2(math.sqrt(a), math.sqrt(1.0 - a)) def fetch_tile(cache_dir: Path, zoom: int, x: int, y: int) -> np.ndarray: cache_dir.mkdir(parents=True, exist_ok=True) path = cache_dir / f"z{zoom}_{x}_{y}.png" if not path.exists(): url = TERRARIUM_URL.format(z=zoom, x=x, y=y) request = urllib.request.Request(url, headers={"User-Agent": USER_AGENT}) with urllib.request.urlopen(request, timeout=60) as response: path.write_bytes(response.read()) with Image.open(path) as image: return np.array(image.convert("RGB"), dtype=np.uint8) def decode_terrarium(rgb: np.ndarray) -> np.ndarray: red = rgb[:, :, 0].astype(np.float32) green = rgb[:, :, 1].astype(np.float32) blue = rgb[:, :, 2].astype(np.float32) return red * 256.0 + green + blue / 256.0 - 32768.0 def block_mean(array: np.ndarray, factor: int) -> np.ndarray: usable_rows = (array.shape[0] // factor) * factor usable_cols = (array.shape[1] // factor) * factor cropped = array[:usable_rows, :usable_cols] reshaped = cropped.reshape( usable_rows // factor, factor, usable_cols // factor, factor, ) return reshaped.mean(axis=(1, 3)) def tile_window( points: Iterable[Point], zoom: int, margin_deg: float, ) -> tuple[tuple[int, ...], tuple[int, ...]]: points = tuple(points) min_lat = min(point.lat for point in points) - margin_deg max_lat = max(point.lat for point in points) + margin_deg min_lon = min(point.lon for point in points) - margin_deg max_lon = max(point.lon for point in points) + margin_deg x0, y1 = deg2num(min_lat, min_lon, zoom) x1, y0 = deg2num(max_lat, max_lon, zoom) xs = tuple(range(min(x0, x1), max(x0, x1) + 1)) ys = tuple(range(min(y0, y1), max(y0, y1) + 1)) return xs, ys def build_distance_tables(lats: np.ndarray, lons: np.ndarray) -> tuple[np.ndarray, np.ndarray]: east_west = np.array( [haversine_km(lat, lons[0], lat, lons[1]) for lat in lats], dtype=np.float32, ) north_south = np.array( [ haversine_km(lats[index], lons[0], lats[index + 1], lons[0]) for index in range(len(lats) - 1) ], dtype=np.float32, ) return east_west, north_south def compute_ruggedness(elev: np.ndarray) -> np.ndarray: padded = np.pad(elev, 1, mode="edge") local_max = padded[1:-1, 1:-1].copy() local_min = padded[1:-1, 1:-1].copy() for dr in (-1, 0, 1): for dc in (-1, 0, 1): window = padded[ 1 + dr : 1 + dr + elev.shape[0], 1 + dc : 1 + dc + elev.shape[1], ] local_max = np.maximum(local_max, window) local_min = np.minimum(local_min, window) return (local_max - local_min).astype(np.float32) def box_mean(array: np.ndarray, radius: int) -> np.ndarray: if radius <= 0: return array.astype(np.float32) padded = np.pad(array.astype(np.float64), radius, mode="edge") integral = np.pad(padded, ((1, 0), (1, 0)), mode="constant").cumsum(axis=0).cumsum(axis=1) window = radius * 2 + 1 summed = ( integral[window:, window:] - integral[:-window, window:] - integral[window:, :-window] + integral[:-window, :-window] ) return (summed / (window * window)).astype(np.float32) def compute_topographic_position(elev: np.ndarray, radius: int = 4) -> np.ndarray: return elev.astype(np.float32) - box_mean(elev, radius) def compute_flow_accumulation( elev: np.ndarray, east_west: np.ndarray, north_south: np.ndarray, ) -> np.ndarray: rows, cols = elev.shape downstream = np.full(rows * cols, -1, dtype=np.int32) offsets = ((-1, 0), (1, 0), (0, -1), (0, 1), (-1, -1), (-1, 1), (1, -1), (1, 1)) for row in range(rows): for col in range(cols): best_gradient = 0.0 best_index = -1 for dr, dc in offsets: next_row = row + dr next_col = col + dc if next_row < 0 or next_row >= rows or next_col < 0 or next_col >= cols: continue if row == next_row: distance_km = float(east_west[row]) elif col == next_col: distance_km = float(north_south[min(row, next_row)]) else: horizontal = float(east_west[min(row, next_row)]) vertical = float(north_south[min(row, next_row)]) distance_km = math.hypot(horizontal, vertical) drop = float(elev[row, col] - elev[next_row, next_col]) if drop <= 0.0: continue gradient = drop / max(distance_km * 1000.0, 1.0) if gradient > best_gradient: best_gradient = gradient best_index = next_row * cols + next_col downstream[row * cols + col] = best_index accumulation = np.ones(rows * cols, dtype=np.float32) for index in np.argsort(elev.ravel())[::-1]: downstream_index = int(downstream[int(index)]) if downstream_index != -1: accumulation[downstream_index] += accumulation[int(index)] return accumulation.reshape(rows, cols) def build_cost_grid( points: Iterable[Point], zoom: int, factor: int, margin_deg: float, cache_dir: Path, ) -> CostGrid: xs, ys = tile_window(points, zoom, margin_deg) rows = [] tile_ids: list[str] = [] for y in ys: row_tiles = [] for x in xs: tile_ids.append(f"z{zoom}_{x}_{y}.png") row_tiles.append(fetch_tile(cache_dir, zoom, x, y)) rows.append(np.concatenate(row_tiles, axis=1)) mosaic = np.concatenate(rows, axis=0) elev = block_mean(decode_terrarium(mosaic), factor) base_x = xs[0] * 256 base_y = ys[0] * 256 col_indices = np.arange(elev.shape[1], dtype=float) * factor + factor / 2.0 row_indices = np.arange(elev.shape[0], dtype=float) * factor + factor / 2.0 lons = np.array([pixel_to_lon(base_x + index, zoom) for index in col_indices], dtype=np.float64) lats = np.array([pixel_to_lat(base_y + index, zoom) for index in row_indices], dtype=np.float64) east_west, north_south = build_distance_tables(lats, lons) ruggedness = compute_ruggedness(elev) topographic_position = compute_topographic_position(elev) flow_accumulation = compute_flow_accumulation(elev, east_west, north_south) return CostGrid( elev=elev, lats=lats, lons=lons, east_west=east_west, north_south=north_south, ruggedness=ruggedness, topographic_position=topographic_position, flow_accumulation=flow_accumulation, tile_ids=tile_ids, ) def nearest_cell(lat: float, lon: float, lats: np.ndarray, lons: np.ndarray) -> tuple[int, int]: return int(np.argmin(np.abs(lats - lat))), int(np.argmin(np.abs(lons - lon))) def bresenham_cells(row0: int, col0: int, row1: int, col1: int) -> list[tuple[int, int]]: cells: list[tuple[int, int]] = [] x0, y0, x1, y1 = col0, row0, col1, row1 dx = abs(x1 - x0) dy = -abs(y1 - y0) sx = 1 if x0 < x1 else -1 sy = 1 if y0 < y1 else -1 error = dx + dy while True: cells.append((y0, x0)) if x0 == x1 and y0 == y1: break doubled = 2 * error if doubled >= dy: error += dy x0 += sx if doubled <= dx: error += dx y0 += sy return cells def scale_modifier(multiplier: float, terrain_weight_scale: float) -> float: """Scale a multiplicative non-Tobler modifier around the neutral value 1.""" return 1.0 + (multiplier - 1.0) * terrain_weight_scale def edge_cost( grid: CostGrid, row: int, col: int, next_row: int, next_col: int, terrain_weight_scale: float, ) -> float: if row == next_row: distance_km = float(grid.east_west[row]) elif col == next_col: distance_km = float(grid.north_south[min(row, next_row)]) else: horizontal = float(grid.east_west[min(row, next_row)]) vertical = float(grid.north_south[min(row, next_row)]) distance_km = math.hypot(horizontal, vertical) rise_m = float(grid.elev[next_row, next_col] - grid.elev[row, col]) slope = rise_m / max(distance_km * 1000.0, 1.0) speed_kmh = 6.0 * math.exp(-3.5 * abs(slope + 0.05)) speed_kmh = max(speed_kmh, 0.35) average_elevation = (float(grid.elev[next_row, next_col]) + float(grid.elev[row, col])) / 2.0 local_relief = ( float(grid.ruggedness[next_row, next_col]) + float(grid.ruggedness[row, col]) ) / 2.0 altitude_penalty = 1.0 if average_elevation > 2000.0: altitude_penalty += min((average_elevation - 2000.0) / 2200.0, 0.85) if abs(slope) > 0.55: altitude_penalty *= 2.0 if average_elevation < -10.0: return math.inf ruggedness_penalty = 1.0 + min(local_relief / 420.0, 1.35) average_tpi = ( float(grid.topographic_position[next_row, next_col]) + float(grid.topographic_position[row, col]) ) / 2.0 landform_penalty = 1.0 if average_tpi < -25.0: landform_penalty *= max(0.84, 1.0 + average_tpi / 650.0) elif average_tpi > 45.0: landform_penalty *= 1.0 + min(average_tpi / 700.0, 0.35) average_flow = ( float(grid.flow_accumulation[next_row, next_col]) + float(grid.flow_accumulation[row, col]) ) / 2.0 log_flow = math.log1p(average_flow) valley_bonus = 1.0 if average_tpi <= 5.0: valley_bonus -= min(max(log_flow - 1.5, 0.0) * 0.04, 0.12) channel_penalty = 1.0 if log_flow > 5.5: channel_penalty += min((log_flow - 5.5) * 0.18, 0.7) return ( distance_km / speed_kmh * scale_modifier(altitude_penalty, terrain_weight_scale) * scale_modifier(ruggedness_penalty, terrain_weight_scale) * scale_modifier(landform_penalty, terrain_weight_scale) * scale_modifier(channel_penalty, terrain_weight_scale) * scale_modifier(valley_bonus, terrain_weight_scale) ) def dijkstra_path( grid: CostGrid, start: tuple[int, int], goal: tuple[int, int], neighbourhood_radius: int, terrain_weight_scale: float, ) -> tuple[list[tuple[int, int]], float]: """Return the exact minimum-cost path on the configured movement graph.""" rows, cols = grid.elev.shape total = rows * cols start_index = start[0] * cols + start[1] goal_index = goal[0] * cols + goal[1] scores = np.full(total, np.inf, dtype=np.float64) parents = np.full(total, -1, dtype=np.int32) closed = np.zeros(total, dtype=bool) queue: list[tuple[float, int]] = [(0.0, start_index)] scores[start_index] = 0.0 neighbours = [ (dr, dc, bresenham_cells(0, 0, dr, dc)) for dr in range(-neighbourhood_radius, neighbourhood_radius + 1) for dc in range(-neighbourhood_radius, neighbourhood_radius + 1) if not (dr == 0 and dc == 0) ] while queue: current_cost, current = heapq.heappop(queue) if closed[current]: continue if current == goal_index: break closed[current] = True row, col = divmod(current, cols) for dr, dc, template in neighbours: next_row = row + dr next_col = col + dc if next_row < 0 or next_row >= rows or next_col < 0 or next_col >= cols: continue next_index = next_row * cols + next_col if closed[next_index]: continue translated = [(row + local_row, col + local_col) for local_row, local_col in template] step_cost = 0.0 valid = True for (segment_row, segment_col), (target_row, target_col) in zip( translated, translated[1:], ): local_cost = edge_cost( grid, segment_row, segment_col, target_row, target_col, terrain_weight_scale, ) if not math.isfinite(local_cost): valid = False break step_cost += local_cost if not valid: continue tentative = current_cost + step_cost if tentative < scores[next_index]: scores[next_index] = tentative parents[next_index] = current heapq.heappush(queue, (tentative, next_index)) if goal_index != start_index and parents[goal_index] == -1: raise RuntimeError(f"No cost-weighted path found between cells {start} and {goal}.") path_indices = [goal_index] while path_indices[-1] != start_index: previous = int(parents[path_indices[-1]]) if previous == -1: raise RuntimeError("Broken parent chain during path reconstruction.") path_indices.append(previous) path_indices.reverse() return [divmod(index, cols) for index in path_indices], float(scores[goal_index]) def rdp(points: list[tuple[float, float]], epsilon: float) -> list[tuple[float, float]]: if len(points) < 3: return points start = np.array(points[0], dtype=float) end = np.array(points[-1], dtype=float) segment = end - start segment_norm = np.linalg.norm(segment) if segment_norm == 0.0: index = 0 maximum = 0.0 else: index = 0 maximum = 0.0 for candidate_index in range(1, len(points) - 1): point = np.array(points[candidate_index], dtype=float) relative = point - start distance = abs(segment[0] * relative[1] - segment[1] * relative[0]) / segment_norm if distance > maximum: maximum = distance index = candidate_index if maximum > epsilon: left = rdp(points[: index + 1], epsilon) right = rdp(points[index:], epsilon) return left[:-1] + right return [points[0], points[-1]] def simplify_polyline( coordinates: list[tuple[float, float]], epsilon: float, ) -> list[tuple[float, float]]: if len(coordinates) <= 2: return coordinates stride = 1 if len(coordinates) > 4000: stride = 4 elif len(coordinates) > 2000: stride = 3 elif len(coordinates) > 1000: stride = 2 reduced = coordinates[::stride] if reduced[-1] != coordinates[-1]: reduced.append(coordinates[-1]) return rdp(reduced, epsilon) def path_to_latlon( path: list[tuple[int, int]], lats: np.ndarray, lons: np.ndarray, ) -> list[tuple[float, float]]: coordinates = [(float(lats[row]), float(lons[col])) for row, col in path] return simplify_polyline(coordinates, 0.0035 if len(coordinates) < 900 else 0.005) def force_exact_endpoints( coordinates: list[tuple[float, float]], source: Point, target: Point, ) -> list[tuple[float, float]]: start = (source.lat, source.lon) end = (target.lat, target.lon) if not coordinates: return [start, end] fixed = list(coordinates) if haversine_km(start[0], start[1], fixed[0][0], fixed[0][1]) > 0.01: fixed.insert(0, start) else: fixed[0] = start if haversine_km(end[0], end[1], fixed[-1][0], fixed[-1][1]) > 0.01: fixed.append(end) else: fixed[-1] = end return fixed def polyline_length(coordinates: list[tuple[float, float]]) -> float: return sum( haversine_km(lat1, lon1, lat2, lon2) for (lat1, lon1), (lat2, lon2) in zip(coordinates, coordinates[1:]) ) def spatial_scenario_status(route: dict) -> str: layer = route.get("layer") if layer == "source_zone_comparison": return "preliminary local source-zone test" if layer == "regional_direct": return "source-attested comparator scenario" return "unsourced source-site test" def run_route( route: dict, point_index: dict[str, Point], zoom: int, factor: int, grid: CostGrid, neighbourhood_radius: int, terrain_weight_scale: float, ) -> dict: source = point_index[route["source_key"]] target = point_index[route["target_key"]] start = nearest_cell(source.lat, source.lon, grid.lats, grid.lons) goal = nearest_cell(target.lat, target.lon, grid.lats, grid.lons) print(f"routing {route['key']} with Dijkstra", flush=True) cell_path, accumulated_cost = dijkstra_path( grid, start, goal, neighbourhood_radius, terrain_weight_scale, ) path = force_exact_endpoints(path_to_latlon(cell_path, grid.lats, grid.lons), source, target) return { "key": route["key"], "name": route["name"], "source_key": route["source_key"], "target_key": route["target_key"], "source_name": route["source_name"], "target_name": route["target_name"], "layer": route.get("layer"), "evidence_filter_key": route.get("evidence_filter_key"), "spatial_scenario_status": spatial_scenario_status(route), "path_status": "deterministic least-cost path scenario; not a reconstructed journey", "path_kind": "anisotropic terrain-conditioned least-cost path", "method_note": ( f"Dijkstra minimum-cost path on a {neighbourhood_radius * 2 + 1} x " f"{neighbourhood_radius * 2 + 1} movement graph, z{zoom}/downsample {factor}; " f"non-Tobler terrain-weight scale {terrain_weight_scale:g}" ), "zoom": zoom, "downsample": factor, "length_km": polyline_length(path), "composite_cost_index": accumulated_cost, "path": [[lat, lon] for lat, lon in path], "tile_ids": grid.tile_ids, } def build_output( source_data: dict, selected_routes: list[dict], zoom: int, factor: int, margin_deg: float, neighbourhood_radius: int, terrain_weight_scale: float, cache_dir: Path, run_label: str, ) -> dict: point_index = { point["key"]: Point( key=point["key"], name=point["name"], lat=float(point["coords"][0]), lon=float(point["coords"][1]), ) for point in source_data["points"] } used_keys = { key for route in selected_routes for key in (route["source_key"], route["target_key"]) } points = [point for point in source_data["points"] if point["key"] in used_keys] model_points = [point_index[point["key"]] for point in points] xs, ys = tile_window(model_points, zoom, margin_deg) print( f"preparing shared z{zoom}/downsample {factor} grid " f"x{xs[0]}-{xs[-1]} y{ys[0]}-{ys[-1]}", flush=True, ) grid = build_cost_grid(model_points, zoom, factor, margin_deg, cache_dir) routes = [ run_route( route, point_index, zoom, factor, grid, neighbourhood_radius, terrain_weight_scale, ) for route in selected_routes ] effective_cell_width_km = float(np.median(grid.east_west)) effective_cell_height_km = float(np.median(grid.north_south)) return { "dataset": "Taş Tepeler obsidian deterministic least-cost path scenarios", "run_label": run_label, "input_route_selection": source_data.get("route_selection"), "algorithm": "Dijkstra's algorithm (zero heuristic)", "optimality_scope": "minimum accumulated cost on the constructed raster movement graph", "interpretive_scope": ( "comparative spatial scenarios only; paths are not evidence of actual Neolithic journeys, " "transfer mechanisms or route use" ), "resolution": { "zoom": zoom, "downsample": factor, "effective_cell_width_km_median": effective_cell_width_km, "effective_cell_height_km_median": effective_cell_height_km, }, "grid_extent": { "tile_x_min": xs[0], "tile_x_max": xs[-1], "tile_y_min": ys[0], "tile_y_max": ys[-1], "shared_by_all_routes_in_run": True, "search_margin_degrees": margin_deg, }, "cost_function": { "walking_speed": "6 × exp[-3.5 × |slope + 0.05|] km/h; floor 0.35 km/h", "altitude": ( "above 2000 m: 1 + min((mean_elevation - 2000) / 2200, 0.85); " "multiplied by 2 where |slope| > 0.55; cells below -10 m are impassable" ), "ruggedness": "1 + min(local_relief / 420, 1.35)", "topographic_position": ( "valley (mean TPI < -25): max(0.84, 1 + mean_TPI / 650); " "ridge (mean TPI > 45): 1 + min(mean_TPI / 700, 0.35)" ), "flow_accumulation": ( "valley bonus where TPI <= 5: 1 - min(max(log1p(flow) - 1.5, 0) × 0.04, 0.12); " "channel penalty where log1p(flow) > 5.5: 1 + min((log1p(flow) - 5.5) × 0.18, 0.7)" ), "non_tobler_terrain_weight_scale": terrain_weight_scale, "neighbourhood": ( f"{neighbourhood_radius * 2 + 1} × {neighbourhood_radius * 2 + 1}, " "with multi-cell moves costed across intervening Bresenham cells" ), "output_units": ( "length_km is geographic path length; composite_cost_index is a relative " "scenario-comparison value, not elapsed time" ), }, "raster_processing": { "source_scheme": "Mapzen Terrarium Web Mercator tiles", "elevation_decode": "R × 256 + G + B / 256 - 32768 metres", "downsampling": "non-overlapping arithmetic block mean", "distance_calculation": "haversine distance between raster-cell centres", "no_data_handling": ( "Mapzen Terrarium RGB tiles have no separate NoData mask in this implementation; " "decoded mean elevations below -10 m are treated as impassable" ), }, "runtime": { "python": platform.python_version(), "numpy": np.__version__, "pillow": Image.__version__, }, "terrain_source": TERRARIUM_URL, "points": points, "routes": routes, } def parse_args() -> argparse.Namespace: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--input", required=True, type=Path, help="Pure point and scenario input JSON") parser.add_argument("--output", required=True, type=Path, help="Primary 14-route output JSON") parser.add_argument( "--sensitivity-output", required=True, type=Path, help="Seven-route multi-parameter sensitivity suite JSON", ) parser.add_argument( "--cache-dir", type=Path, default=Path(__file__).resolve().parent / "terrain_cache", help="Mapzen Terrarium tile cache (downloaded if absent)", ) return parser.parse_args() def main() -> None: args = parse_args() with args.input.open(encoding="utf-8") as handle: source_data = json.load(handle) routes = source_data.get("routes", []) if len(routes) != 14: raise ValueError(f"Expected 14 article routes, found {len(routes)}") primary = build_output( source_data, routes, PRIMARY_ZOOM, PRIMARY_DOWNSAMPLE, MARGIN_DEG, 2, 1.0, args.cache_dir, "primary common-resolution rerun", ) sensitivity_routes = [route for route in routes if route["key"] in SENSITIVITY_KEYS] if [route["key"] for route in sensitivity_routes] != list(SENSITIVITY_KEYS): missing = sorted(set(SENSITIVITY_KEYS) - {route["key"] for route in sensitivity_routes}) raise ValueError(f"Missing headline sensitivity routes: {missing}") sensitivity_runs = [] for variant in SENSITIVITY_VARIANTS: sensitivity_runs.append( build_output( source_data, sensitivity_routes, variant["zoom"], variant["downsample"], variant["margin_deg"], variant["neighbourhood_radius"], variant["terrain_weight_scale"], args.cache_dir, variant["key"], ) ) sensitivity = { "dataset": "Taş Tepeler obsidian least-cost path sensitivity suite", "baseline_run": { "zoom": PRIMARY_ZOOM, "downsample": PRIMARY_DOWNSAMPLE, "margin_deg": MARGIN_DEG, "neighbourhood_radius": 2, "terrain_weight_scale": 1.0, }, "scope": ( "seven headline paths; tests grid resolution, search margin, movement neighbourhood " "and the strength of non-Tobler terrain modifiers; endpoint placement remains untested" ), "runs": sensitivity_runs, } args.output.parent.mkdir(parents=True, exist_ok=True) args.sensitivity_output.parent.mkdir(parents=True, exist_ok=True) args.output.write_text(json.dumps(primary, indent=2, ensure_ascii=False) + "\n", encoding="utf-8") args.sensitivity_output.write_text( json.dumps(sensitivity, indent=2, ensure_ascii=False) + "\n", encoding="utf-8", ) print(f"wrote {args.output}") print(f"wrote {args.sensitivity_output}") if __name__ == "__main__": main()