Function bodies 32 total

def parse_args(argv: Sequence[str] | None = None) -> argparse.Namespace:
    """Parse command-line arguments."""
    parser = argparse.ArgumentParser(
        prog="terrain-to-stl",
        description="Generate 3D-printable STL terrain tiles from real-world elevation data.",
    )

    # --- mutually exclusive region source ---
    source = parser.add_mutually_exclusive_group()
    source.add_argument(
        "--bbox",
        type=str,
        default=None,
        help="Bounding box as 'south,west,north,east' (decimal lat/lon)",
    )
    source.add_argument(
        "--landmark",
        type=str,
        default=None,
        help="Landmark name to geocode (e.g. 'Manhattan')",
    )
    source.add_argument(
        "--config",
        type=Path,
        default=None,
        help="Path to a JSON config file with all parameters",
    )

    # --- tunables ---
    parser.add_argument("--scale", type=int, default=3333, help="Horizontal scale factor (default: 3333)")
    parser.add_a

def build_config(args: argparse.Namespace) -> Config:
    """Convert parsed CLI args into a Config object."""
    if args.config is not None:
        return Config.from_json(args.config)

    kwargs: dict = dict(
        scale=args.scale,
        vert_exag=args.exag,
        buffer_m=args.buffer,
        tile_rows=args.tile_rows,
        tile_cols=args.tile_cols,
        build_plate_x_mm=args.plate_x,
        build_plate_y_mm=args.plate_y,
        max_print_height_mm=args.max_height,
        output_dir=args.output,
    )

    if args.bbox is not None:
        parts = [float(v) for v in args.bbox.split(",")]
        if len(parts) != 4:
            raise SystemExit("--bbox requires exactly 4 comma-separated values: south,west,north,east")
        south, west, north, east = parts
        kwargs.update(bbox_south=south, bbox_west=west, bbox_north=north, bbox_east=east)

    elif args.landmark is not None:
        # Lazy import so missing fetch_data doesn't break --help
        from terrain

def main(argv: Sequence[str] | None = None) -> None:
    """Run the full terrain-to-stl pipeline."""
    args = parse_args(argv)
    cfg = build_config(args)

    # ---- summary ----
    print("=" * 60)
    print("terrain-to-stl")
    print("=" * 60)
    print(f"  Scale          : 1:{cfg.scale}")
    print(f"  Vert exag      : {cfg.vert_exag}x")
    print(f"  Tile grid      : {cfg.tile_rows} rows x {cfg.tile_cols} cols")
    print(f"  Supertall      : >{cfg.supertall_threshold_m:.0f} m")
    print(f"  Output         : {cfg.output_dir}")
    print("=" * 60)

    # ---- pipeline stages (lazy imports) ----
    from terrain_to_stl.fetch_data import fetch_all  # type: ignore[import-not-found]

    fetch_all(cfg, force=args.force)

    from terrain_to_stl.process import process_all  # type: ignore[import-not-found]

    process_all(cfg)

    from terrain_to_stl.generate import generate_all  # type: ignore[import-not-found]

    generate_all(cfg)

    from terrain_to_stl.tile import split_all

class Config:
    """All tunable parameters for terrain-to-stl pipeline."""

    # Geographic bounding box
    bbox_north: float = 40.800
    bbox_south: float = 40.764
    bbox_east: float = -73.949
    bbox_west: float = -73.981
    buffer_m: float = 250.0

    # Scale & exaggeration
    scale: int = 3333
    vert_exag: float = 4.0

    # Build plate dimensions
    build_plate_x_mm: float = 230.0
    build_plate_y_mm: float = 190.0
    max_print_height_mm: float = 200.0

    # Tiling
    tile_cols: int = 2
    tile_rows: int = 7

    # Base
    base_thickness_mm: float = 5.0

    # Peg / socket alignment
    peg_diameter_mm: float = 4.0
    peg_depth_mm: float = 6.0
    socket_clearance_mm: float = 0.3
    pegs_per_edge: int = 2

    # Building setbacks
    setback_floor_1: int = 6
    setback_floor_2: int = 20
    setback_ratio: float = 0.15

    # Supertall building pegs
    supertall_peg_diameter_mm: float = 3.0
    supertall_peg_depth_mm: float = 8.0
    supertall_clearance_mm: f

    def real_to_model_mm(self, meters: float) -> float:
        """Convert real-world horizontal distance (m) to model mm."""
        return meters / self.scale * 1000

    def real_to_model_mm_vertical(self, meters: float) -> float:
        """Convert real-world height (m) to model mm with vertical exaggeration."""
        return meters / self.scale * 1000 * self.vert_exag

    def from_json(cls, path: str | Path) -> Config:
        """Load a JSON file and apply overrides on top of defaults."""
        with open(path) as f:
            overrides = json.load(f)
        # Convert output_dir string to Path if present
        if "output_dir" in overrides:
            overrides["output_dir"] = Path(overrides["output_dir"])
        return cls(**overrides)

def geocode_landmark(name: str) -> dict:
    """Geocode a landmark name to a bounding-box dict via Nominatim.

    Returns
    -------
    dict
        Keys: north, south, east, west (floats).

    Raises
    ------
    ValueError
        If Nominatim returns no results for *name*.
    """
    url = "https://nominatim.openstreetmap.org/search"
    params = {
        "q": name,
        "format": "json",
        "limit": 1,
    }
    headers = {"User-Agent": "terrain-to-stl/0.1"}
    resp = requests.get(url, params=params, headers=headers, timeout=10)
    resp.raise_for_status()

    results = resp.json()
    if not results:
        raise ValueError(f"No geocoding results for '{name}'")

    hit = results[0]
    bb = hit["boundingbox"]  # [south, north, west, east] as strings
    return {
        "south": float(bb[0]),
        "north": float(bb[1]),
        "west": float(bb[2]),
        "east": float(bb[3]),
    }

def fetch_dem(bbox: dict, output_dir: Path, force: bool = False) -> Path:
    """Download a USGS 3DEP DEM GeoTIFF covering *bbox*.

    Parameters
    ----------
    bbox : dict
        north / south / east / west.
    output_dir : Path
        Directory for cached files.
    force : bool
        Re-download even if the file already exists.

    Returns
    -------
    Path
        Path to the downloaded ``dem.tif``.
    """
    output_dir = Path(output_dir)
    output_dir.mkdir(parents=True, exist_ok=True)
    dem_path = output_dir / "dem.tif"

    if dem_path.exists() and not force:
        print(f"  DEM cached: {dem_path} ({dem_path.stat().st_size / 1024:.0f} KB)")
        return dem_path

    tnm_url = "https://tnmaccess.nationalmap.gov/api/v1/products"

    datasets = [
        "National Elevation Dataset (NED) 1/3 arc-second",
        "3DEP 1/3 arc-second",
    ]

    download_url: str | None = None

    for dataset in datasets:
        print(f"  Querying TNM for dataset: {datase

def fetch_nyc_buildings(
    bbox: dict,
    buffer_m: float,
    output_dir: Path,
    force: bool = False,
) -> Path:
    """Download building footprints from NYC OpenData (Socrata).

    Parameters
    ----------
    bbox : dict
        north / south / east / west.
    buffer_m : float
        Buffer in metres (used for context but not applied to the spatial
        query itself — pass the already-buffered bbox if needed).
    output_dir : Path
        Directory for cached files.
    force : bool
        Re-download even if the file already exists.

    Returns
    -------
    Path
        Path to the saved ``buildings.geojson``.
    """
    output_dir = Path(output_dir)
    output_dir.mkdir(parents=True, exist_ok=True)
    out_path = output_dir / "buildings.geojson"

    if out_path.exists() and not force:
        print(f"  Buildings cached: {out_path} ({out_path.stat().st_size / 1024:.0f} KB)")
        return out_path

    endpoint = "https://data.cityofnewyork.us/resource/qb5r-6d

def fetch_nyc_water(
    bbox: dict,
    output_dir: Path,
    force: bool = False,
) -> Path:
    """Download water-body polygons from NYC OpenData (Socrata).

    Parameters
    ----------
    bbox : dict
        north / south / east / west.
    output_dir : Path
        Directory for cached files.
    force : bool
        Re-download even if the file already exists.

    Returns
    -------
    Path
        Path to the saved ``water.geojson``.
    """
    output_dir = Path(output_dir)
    output_dir.mkdir(parents=True, exist_ok=True)
    out_path = output_dir / "water.geojson"

    if out_path.exists() and not force:
        print(f"  Water cached: {out_path} ({out_path.stat().st_size / 1024:.0f} KB)")
        return out_path

    endpoint = "https://data.cityofnewyork.us/resource/drh3-e2fd.geojson"
    where_clause = (
        f"within_box(the_geom, "
        f"{bbox['south']}, {bbox['west']}, "
        f"{bbox['north']}, {bbox['east']})"
    )
    params = {
        "$where": wher

def fetch_all(
    config: Config,
    output_dir: Path,
    force: bool = False,
) -> tuple[Path, Path, Path]:
    """Fetch DEM, buildings, and water data for the configured bbox.

    Parameters
    ----------
    config : Config
        Pipeline configuration (supplies bbox and buffer_m).
    output_dir : Path
        Root directory for cached downloads.
    force : bool
        Re-download everything even if cached.

    Returns
    -------
    tuple[Path, Path, Path]
        (dem_path, buildings_path, water_path)
    """
    print("Step 1/4: Fetching data...")
    output_dir = Path(output_dir)
    output_dir.mkdir(parents=True, exist_ok=True)

    bbox = config.bbox

    dem_path = fetch_dem(bbox, output_dir, force=force)
    buildings_path = fetch_nyc_buildings(bbox, config.buffer_m, output_dir, force=force)
    water_path = fetch_nyc_water(bbox, output_dir, force=force)

    print("  All data fetched.")
    return dem_path, buildings_path, water_path

def heightmap_to_mesh(
    heightmap: np.ndarray,
    extent_x_mm: float,
    extent_y_mm: float,
    base_thickness_mm: float = 5.0,
) -> trimesh.Trimesh:
    """Build a watertight solid mesh from a 2-D heightmap array.

    Top surface: triangulated grid (2 triangles per cell).
    Bottom surface: flat plane at z = -base_thickness_mm.
    Side walls: connect top and bottom perimeter edges.

    Parameters
    ----------
    heightmap : 2-D array of elevation values (mm).
    extent_x_mm, extent_y_mm : physical size of the model (mm).
    base_thickness_mm : thickness of the flat base below z = 0.

    Returns
    -------
    trimesh.Trimesh – watertight terrain solid.
    """
    rows, cols = heightmap.shape
    dx = extent_x_mm / (cols - 1)
    dy = extent_y_mm / (rows - 1)

    # --- Top surface vertices ---
    # Row 0 = north → y = (rows-1)*dy, row (rows-1) = south → y = 0
    top_verts = []
    for r in range(rows):
        for c in range(cols):
            x = c * dx
          

def punch_water_holes(
    mesh: trimesh.Trimesh,
    water_mask: np.ndarray,
    extent_x_mm: float,
    extent_y_mm: float,
) -> trimesh.Trimesh:
    """Boolean-subtract water cells from the terrain mesh.

    Parameters
    ----------
    mesh : terrain mesh.
    water_mask : boolean 2-D array (True == water).
    extent_x_mm, extent_y_mm : physical extents of the model.

    Returns
    -------
    trimesh.Trimesh with water regions removed.
    """
    rows, cols = water_mask.shape
    dx = extent_x_mm / (cols - 1) if cols > 1 else extent_x_mm
    dy = extent_y_mm / (rows - 1) if rows > 1 else extent_y_mm
    cell_area = dx * dy

    # Build shapely boxes for each True cell
    boxes = []
    for r in range(rows):
        for c in range(cols):
            if water_mask[r, c]:
                x0 = c * dx - dx / 2
                y0 = (rows - 1 - r) * dy - dy / 2
                boxes.append(box(x0, y0, x0 + dx, y0 + dy))

    if not boxes:
        return mesh

    merged = unary_un

def extrude_building(
    footprint_mm: list[tuple[float, float]],
    base_z_mm: float,
    height_mm: float,
    num_floors: int | None,
    setback_floor_1: int = 6,
    setback_floor_2: int = 20,
    setback_ratio: float = 0.15,
) -> trimesh.Trimesh:
    """Create a building mesh with optional setback tiers.

    Parameters
    ----------
    footprint_mm : exterior ring coordinates in model mm.
    base_z_mm : ground elevation in model mm.
    height_mm : total building height in model mm.
    num_floors : number of floors (None → simple extrusion).
    setback_floor_1 : floor at which the first setback occurs.
    setback_floor_2 : floor at which the second setback occurs.
    setback_ratio : fraction of perimeter used for setback offset.

    Returns
    -------
    trimesh.Trimesh – building solid.
    """
    poly = Polygon(footprint_mm)
    if not poly.is_valid:
        poly = poly.buffer(0)

    num_floors = num_floors or 1

    # Simple case: no setbacks
    if num_floors <

def generate_supertall_piece(
    footprint_mm: list[tuple[float, float]],
    base_z_mm: float,
    height_mm: float,
    num_floors: int | None,
    setback_floor_1: int = 6,
    setback_floor_2: int = 20,
    setback_ratio: float = 0.15,
    peg_diameter_mm: float = 3.0,
    peg_depth_mm: float = 8.0,
    clearance_mm: float = 0.25,
) -> tuple[trimesh.Trimesh, trimesh.Trimesh]:
    """Generate a supertall building as a removable piece with alignment peg.

    Parameters
    ----------
    footprint_mm : exterior ring coordinates in model mm.
    base_z_mm : ground elevation on the terrain (mm).
    height_mm : total building height (mm).
    num_floors : number of floors.
    setback_floor_1, setback_floor_2, setback_ratio : setback parameters.
    peg_diameter_mm : diameter of the alignment peg.
    peg_depth_mm : length of the peg extending downward.
    clearance_mm : radial clearance added to the socket hole.

    Returns
    -------
    (piece, socket_hole) — the printable buil

def generate_all(
    terrain_model: dict[str, Any],
    water_mask: np.ndarray,
    buildings: list[dict],
    config: Config,
    stl_dir: str | Path,
) -> tuple[trimesh.Trimesh, list[dict]]:
    """Generate all meshes: terrain, water cutouts, buildings, supertalls.

    Parameters
    ----------
    terrain_model : dict with keys heightmap, extent_x_mm, extent_y_mm.
    water_mask : boolean 2-D array (True == water).
    buildings : list of building dicts from process_terrain.
    config : project Config.
    stl_dir : directory for STL output files.

    Returns
    -------
    (combined_mesh, supertall_info_list).
    """
    print("Step 3/4: Generating meshes...")

    stl_dir = Path(stl_dir)
    stl_dir.mkdir(parents=True, exist_ok=True)

    heightmap = terrain_model["heightmap"]
    extent_x_mm = terrain_model["extent_x_mm"]
    extent_y_mm = terrain_model["extent_y_mm"]

    # --- Terrain mesh ---
    terrain_mesh = heightmap_to_mesh(
        heightmap,
        extent_x_mm,
 

def load_and_crop_dem(
    dem_path: str | Path,
    bbox: dict[str, float],
) -> tuple[np.ndarray, dict[str, Any]]:
    """Open a GeoTIFF DEM, crop to *bbox*, and reproject to EPSG:2263.

    Parameters
    ----------
    dem_path : path to a GeoTIFF elevation raster.
    bbox : dict with keys north, south, east, west (WGS-84 degrees).

    Returns
    -------
    (reprojected_array, metadata) where metadata contains:
        extent_x_m, extent_y_m, transform, crs, shape.
    """
    dem_path = Path(dem_path)

    with rasterio.open(dem_path) as src:
        # Crop to bounding box
        window = from_bounds(
            bbox["west"],
            bbox["south"],
            bbox["east"],
            bbox["north"],
            src.transform,
        )
        cropped = src.read(1, window=window)
        cropped_transform = rasterio.windows.transform(window, src.transform)
        src_crs = src.crs

    # Compute destination transform & shape for reprojection
    dst_transform, dst_widt

def apply_vertical_exaggeration(terrain: np.ndarray, exag: float) -> np.ndarray:
    """Zero-base the terrain and multiply by the exaggeration factor.

    Parameters
    ----------
    terrain : 2-D elevation array (real-world metres).
    exag : vertical exaggeration multiplier (e.g. 4.0).

    Returns
    -------
    Exaggerated elevation array with minimum == 0.
    """
    zeroed = terrain - terrain.min()
    return zeroed * exag

def terrain_to_model_coords(
    terrain: np.ndarray,
    real_extent_x_m: float,
    real_extent_y_m: float,
    config: Config,
) -> dict[str, Any]:
    """Scale a real-world heightmap into model (mm) coordinates.

    Parameters
    ----------
    terrain : 2-D elevation array (real-world metres, already exaggerated).
    real_extent_x_m, real_extent_y_m : footprint size in metres.
    config : project Config (supplies scale + vert_exag).

    Returns
    -------
    dict with keys: heightmap (ndarray in mm), extent_x_mm, extent_y_mm.
    """
    heightmap_mm = np.vectorize(config.real_to_model_mm_vertical)(terrain)
    extent_x_mm = config.real_to_model_mm(real_extent_x_m)
    extent_y_mm = config.real_to_model_mm(real_extent_y_m)

    return {
        "heightmap": heightmap_mm,
        "extent_x_mm": extent_x_mm,
        "extent_y_mm": extent_y_mm,
    }

def create_water_mask(
    water_path: str | Path,
    terrain_shape: tuple[int, int],
    transform: rasterio.transform.Affine,
    crs: str,
) -> np.ndarray:
    """Rasterize a water-polygon GeoJSON onto the terrain grid.

    Parameters
    ----------
    water_path : path to a GeoJSON file with water body polygons.
    terrain_shape : (rows, cols) of the target raster.
    transform : affine transform of the target raster.
    crs : CRS string of the target raster.

    Returns
    -------
    Boolean ndarray (True == water).
    """
    water_path = Path(water_path)

    gdf = gpd.read_file(water_path)
    gdf = gdf.to_crs(crs)

    # Build (geometry, value) pairs for rasterize
    shapes = [(mapping(geom), 1) for geom in gdf.geometry if geom is not None]

    if not shapes:
        return np.zeros(terrain_shape, dtype=bool)

    rasterized = rasterio.features.rasterize(
        shapes,
        out_shape=terrain_shape,
        transform=transform,
        fill=0,
        dtype="ui

def process_buildings(
    buildings_path: str | Path,
    terrain: np.ndarray,
    metadata: dict[str, Any],
    config: Config,
) -> list[dict]:
    """Load building footprints, extract heights, and sample base elevations.

    Parameters
    ----------
    buildings_path : path to a GeoJSON with building footprints.
    terrain : 2-D DEM array (reprojected, in target CRS).
    metadata : dict from ``load_and_crop_dem`` (must have transform, crs, shape).
    config : project Config.

    Returns
    -------
    List of dicts with keys:
        name, footprint_coords, base_elev_m, height_m, num_floors, is_supertall.
    """
    buildings_path = Path(buildings_path)
    gdf = gpd.read_file(buildings_path)
    gdf = gdf.to_crs(metadata["crs"])

    transform = metadata["transform"]
    inv_transform = ~transform  # inverse affine

    results: list[dict] = []

    for idx, row in gdf.iterrows():
        geom = row.geometry
        if geom is None:
            continue

        # Height:

def classify_buildings(
    buildings: list[dict],
    config: Config,
) -> list[dict]:
    """Mark buildings as supertall when height exceeds threshold.

    Parameters
    ----------
    buildings : list of building dicts (must contain ``height_m``).
    config : project Config (supplies ``supertall_threshold_m``).

    Returns
    -------
    Same list with ``is_supertall`` updated in-place (also returned).
    """
    threshold = config.supertall_threshold_m
    for bldg in buildings:
        bldg["is_supertall"] = bldg["height_m"] > threshold
    return buildings

def process_all(
    dem_path: str | Path,
    buildings_path: str | Path,
    water_path: str | Path,
    config: Config,
    output_dir: str | Path,
) -> tuple[dict, np.ndarray, list[dict]]:
    """Run the full terrain-processing pipeline.

    Parameters
    ----------
    dem_path : path to DEM GeoTIFF.
    buildings_path : path to buildings GeoJSON.
    water_path : path to water GeoJSON.
    config : project Config.
    output_dir : directory to write intermediate files.

    Returns
    -------
    (terrain_model_dict, water_mask, buildings_list).
    """
    output_dir = Path(output_dir)
    output_dir.mkdir(parents=True, exist_ok=True)

    # Step 1: Load & crop DEM
    print("Step 2/4: Processing terrain...")
    terrain_raw, metadata = load_and_crop_dem(dem_path, config.bbox)

    # Step 2: Vertical exaggeration
    terrain_exag = apply_vertical_exaggeration(terrain_raw, config.vert_exag)

    # Step 3: Convert to model coordinates
    terrain_model = terrain_to_model_coords

def compute_tile_bounds(
    total_x: float,
    total_y: float,
    rows: int,
    cols: int,
) -> list[dict]:
    """Divide a rectangular area into a grid of tile bounds.

    Parameters
    ----------
    total_x : total extent in the X (east-west) direction.
    total_y : total extent in the Y (north-south) direction.
    rows : number of rows (north-south divisions).
    cols : number of columns (east-west divisions).

    Returns
    -------
    List of dicts with keys: row, col, x_min, x_max, y_min, y_max.
    Row 0 = south, col 0 = west.
    """
    tile_w = total_x / cols
    tile_h = total_y / rows

    tiles: list[dict] = []
    for r in range(rows):
        for c in range(cols):
            tiles.append(
                {
                    "row": r,
                    "col": c,
                    "x_min": c * tile_w,
                    "x_max": (c + 1) * tile_w,
                    "y_min": r * tile_h,
                    "y_max": (r + 1) * tile_h,
                }
  

def slice_mesh_to_tile(
    mesh: trimesh.Trimesh,
    x_min: float,
    x_max: float,
    y_min: float,
    y_max: float,
) -> trimesh.Trimesh:
    """Clip a mesh to the given XY tile bounds.

    Creates a clipping box that covers the full Z range of the mesh
    (plus 10 mm margin) and uses boolean intersection to extract the tile.

    Parameters
    ----------
    mesh : source mesh to clip.
    x_min, x_max : tile bounds in X.
    y_min, y_max : tile bounds in Y.

    Returns
    -------
    Clipped trimesh.Trimesh.
    """
    z_min = mesh.bounds[0][2] - 10.0
    z_max = mesh.bounds[1][2] + 10.0

    center = [
        (x_min + x_max) / 2,
        (y_min + y_max) / 2,
        (z_min + z_max) / 2,
    ]
    extents = [
        x_max - x_min,
        y_max - y_min,
        z_max - z_min,
    ]

    clip_box = trimesh.creation.box(
        extents=extents,
        transform=trimesh.transformations.translation_matrix(center),
    )
    return mesh.intersection(clip_box)

def create_alignment_peg(
    x: float,
    y: float,
    z_bottom: float,
    z_top: float,
    diameter: float,
    depth: float,
) -> trimesh.Trimesh:
    """Create a cylindrical alignment peg at (x, y).

    Parameters
    ----------
    x, y : horizontal position of the peg centre.
    z_bottom : Z coordinate of the base of the peg.
    z_top : unused (kept for interface symmetry).
    diameter : peg diameter.
    depth : peg height (length along Z).

    Returns
    -------
    trimesh.Trimesh cylinder.
    """
    radius = diameter / 2.0
    peg = trimesh.creation.cylinder(radius=radius, height=depth)
    # Centre the peg so its base sits at z_bottom
    peg.apply_translation([x, y, z_bottom + depth / 2.0])
    return peg

def create_alignment_hole(
    x: float,
    y: float,
    z_bottom: float,
    diameter: float,
    depth: float,
    clearance: float,
) -> trimesh.Trimesh:
    """Create a cylinder for boolean-subtracting an alignment hole.

    The hole is slightly wider (clearance on each side) and slightly
    deeper than the corresponding peg.

    Parameters
    ----------
    x, y : horizontal position.
    z_bottom : Z coordinate for the base of the hole.
    diameter : nominal peg diameter.
    depth : nominal peg depth.
    clearance : radial clearance added to the radius.

    Returns
    -------
    trimesh.Trimesh cylinder (to be subtracted from the tile mesh).
    """
    radius = diameter / 2.0 + clearance
    hole_depth = depth + 1.0  # slightly deeper for easy insertion
    hole = trimesh.creation.cylinder(radius=radius, height=hole_depth)
    hole.apply_translation([x, y, z_bottom + hole_depth / 2.0])
    return hole

def add_alignment_features(
    tiles: list[dict],
    config: Config,
) -> list[dict]:
    """Add pegs and holes to tile meshes for assembly alignment.

    For each tile:
    - North/east edges that border another tile get **pegs** (concatenated).
    - South/west edges that border another tile get **holes** (boolean subtracted).

    Pegs are placed at 1/3 and 2/3 of the edge length.

    Parameters
    ----------
    tiles : list of tile dicts; each must have keys
        ``row``, ``col``, ``mesh``, ``x_min``, ``x_max``, ``y_min``, ``y_max``.
    config : project Config (peg specs).

    Returns
    -------
    Same list with ``mesh`` values modified in-place.
    """
    max_row = max(t["row"] for t in tiles)
    max_col = max(t["col"] for t in tiles)

    diameter = config.peg_diameter_mm
    depth = config.peg_depth_mm
    clearance = config.socket_clearance_mm
    n_pegs = config.pegs_per_edge

    for tile in tiles:
        mesh = tile["mesh"]
        z_bottom = mesh.bounds[0]

def generate_assembly_guide(
    tiles: list[dict],
    supertall_info: list[dict],
    config: Config,
) -> str:
    """Generate a plain-text assembly guide.

    Parameters
    ----------
    tiles : list of tile dicts with row, col, x_min, x_max, y_min, y_max.
    supertall_info : list of dicts with ``name`` and ``file`` keys.
    config : project Config.

    Returns
    -------
    Multi-line string with grid layout, specs, and printing tips.
    """
    max_row = max(t["row"] for t in tiles)
    max_col = max(t["col"] for t in tiles)
    rows = max_row + 1
    cols = max_col + 1

    # Tile dimensions (from first tile)
    t0 = tiles[0]
    tile_w = t0["x_max"] - t0["x_min"]
    tile_h = t0["y_max"] - t0["y_min"]

    lines: list[str] = []
    lines.append("=" * 60)
    lines.append("TERRAIN-TO-STL ASSEMBLY GUIDE")
    lines.append("=" * 60)
    lines.append("")

    # Grid layout (north at top)
    lines.append("GRID LAYOUT (north at top):")
    lines.append("-" * 40)
    for r 

def render_previews(
    full_mesh: trimesh.Trimesh,
    tiles: list[dict],
    config: Config,
    output_dir: str | Path,
) -> None:
    """Render a top-down tile grid preview and save as PNG.

    Parameters
    ----------
    full_mesh : the complete (unsplit) mesh.
    tiles : list of tile dicts with row, col, and bounds.
    config : project Config.
    output_dir : directory for the output image.
    """
    import matplotlib

    matplotlib.use("Agg")
    import matplotlib.patches as mpatches
    import matplotlib.pyplot as plt

    output_dir = Path(output_dir)
    output_dir.mkdir(parents=True, exist_ok=True)

    fig, ax = plt.subplots(1, 1, figsize=(10, 8))
    ax.set_aspect("equal")
    ax.set_title("Tile Grid Layout (top-down view)")
    ax.set_xlabel("X (mm)")
    ax.set_ylabel("Y (mm)")

    for tile in tiles:
        x0 = tile["x_min"]
        y0 = tile["y_min"]
        w = tile["x_max"] - tile["x_min"]
        h = tile["y_max"] - tile["y_min"]

        rect = mpatches

def split_all(
    full_mesh: trimesh.Trimesh,
    supertall_info: list[dict],
    buildings: list[dict],
    config: Config,
    stl_dir: str | Path,
    preview_dir: str | Path,
) -> None:
    """Run the full tile-splitting pipeline.

    Parameters
    ----------
    full_mesh : complete terrain mesh.
    supertall_info : list of supertall building dicts.
    buildings : list of all building dicts.
    config : project Config.
    stl_dir : directory for STL output files.
    preview_dir : directory for preview images.
    """
    print("Step 4/4: Splitting into tiles...")

    stl_dir = Path(stl_dir)
    preview_dir = Path(preview_dir)
    stl_dir.mkdir(parents=True, exist_ok=True)
    preview_dir.mkdir(parents=True, exist_ok=True)

    # Compute total extents from mesh bounds
    total_x = full_mesh.bounds[1][0] - full_mesh.bounds[0][0]
    total_y = full_mesh.bounds[1][1] - full_mesh.bounds[0][1]

    # Compute tile bounds
    tile_bounds = compute_tile_bounds(
        total_x, t