gig-line/Tools/blosm_blockout.py

"""Headless Blosm blockout: import an OSM extract into a .blend, no GUI.

Run with Blender in background mode. The add-on path has to point at the folder
that contains `blosm/` (for the Flatpak install that is
`~/.var/app/org.blender.Blender/config/blender/5.2/scripts`):

    B=/var/lib/flatpak/app/org.blender.Blender/x86_64/stable/active/files
    LD_LIBRARY_PATH="$B/blender/lib" \
    BLENDER_USER_SCRIPTS="$HOME/.var/app/org.blender.Blender/config/blender/5.2/scripts" \
    "$B/blender/blender" -b --factory-startup \
        --python Tools/blosm_blockout.py -- \
        --osm Geo/raw/osm-api-map-2026-11-06.osm \
        --bounds 34.1006,-118.3422,34.1020,-118.3370 \
        --out Art/blender/blockout/blockout-2026-11-06.blend

What it does:

1. enables Blosm and checks that `scene.blosm` exists;
2. imports OSM data from a local `.osm` file (no network, no Overpass);
3. reports what arrived, per import category and in metres, so a scale or
   missing-category failure is visible in the log rather than at first render;
4. saves the .blend.

No source assets are modified in place: the output path is the only write.
"""

import argparse
import math
import os
import sys

import addon_utils
import bpy


def parse_args():
    argv = sys.argv
    argv = argv[argv.index("--") + 1 :] if "--" in argv else []
    ap = argparse.ArgumentParser()
    ap.add_argument("--osm", required=True, help="local .osm file to import")
    ap.add_argument("--bounds", required=True, help="south,west,north,east in WGS84 degrees")
    ap.add_argument("--out", required=True, help="output .blend path")
    ap.add_argument("--level-height", type=float, default=3.2, help="metres per building level")
    ap.add_argument("--origin", default=None,
                    help="lat,lon to use as the Blender origin (project origin). Default: Blosm centres on the extract")
    ap.add_argument("--no-subdivide", action="store_true", help="skip OSM area subdivision")
    ap.add_argument("--load-missing", action="store_true",
                    help="let Blosm fetch relation members missing from the extract (needs network, breaks offline reproducibility)")
    ap.add_argument("--mode", default="3Dsimple", choices=("3Dsimple", "3Drealistic", "2D"),
                    help="Blosm import mode; 3Drealistic needs the paid material asset pack")
    ap.add_argument("--assets-dir", default=None,
                    help="directory with building_materials.blend/vegetation.blend, only needed for 3Drealistic")
    ap.add_argument("--data-dir", default="Art/blender/.cache/blosm",
                    help="where Blosm may cache downloads")
    return ap.parse_args(argv)


def enable_blosm():
    present = [m.__name__ for m in addon_utils.modules() if m.__name__ == "blosm"]
    if not present:
        raise SystemExit(
            "Blosm is not on Blender's script path. Set BLENDER_USER_SCRIPTS to the "
            "scripts folder that contains blosm/, e.g. "
            "$HOME/.var/app/org.blender.Blender/config/blender/5.2/scripts"
        )
    addon_utils.enable("blosm", default_set=True, persistent=False)
    if not hasattr(bpy.context.scene, "blosm"):
        raise SystemExit("Blosm enabled but scene.blosm is missing; add-on registration failed")


def configure_preferences(args):
    """Blosm refuses to import until an asset directory is set.

    The check only requires the directory to exist; the paid asset pack
    (building_materials.blend, vegetation.blend) supplies richer materials, and
    without it the blockout still gets geometry and placeholder materials.
    """
    import blosm

    prefs = bpy.context.preferences.addons["blosm"].preferences
    assets = args.assets_dir or os.path.join(os.path.dirname(blosm.__file__), "assets")
    if os.path.isdir(assets):
        prefs.assetsDir = assets
    elif args.mode == "3Drealistic":
        raise SystemExit(f"3Drealistic needs an asset directory; {assets} does not exist")
    data_dir = os.path.abspath(args.data_dir)
    os.makedirs(data_dir, exist_ok=True)
    prefs.dataDir = data_dir
    # In background mode the GUI never populates the default level table, and the
    # add-on then takes a code path that imports a module relative to the wrong
    # package ("from gui import addDefaultLevels") and dies. One entry avoids it.
    scene = bpy.context.scene
    levels = scene.blosm.defaultLevels
    if len(levels) == 0:
        item = levels.add()
        item.levels = 4
        item.weight = 1
        print("seeded one default level entry (headless guard)")
    print(f"blosm assets: {assets}")
    print(f"blosm data:   {data_dir}")


def import_osm(args):
    south, west, north, east = (float(v) for v in args.bounds.split(","))
    scene = bpy.context.scene
    osmp = scene.blosm

    if args.origin:
        lat, lon = (float(v) for v in args.origin.split(","))
        scene["lat"] = lat
        scene["lon"] = lon
        osmp.relativeToInitialImport = True
        print(f"origin set to {lat}, {lon} (project origin)")

    osmp.commandLineMode = True
    osmp.mode = args.mode
    assert osmp.mode == args.mode
    osmp.dataType = "osm"
    osmp.osmSource = "file"
    osmp.osmFilepath = os.path.abspath(args.osm)
    osmp.coordinatesAsFilter = True
    osmp.minLat = south
    osmp.maxLat = north
    osmp.minLon = west
    osmp.maxLon = east

    # what to bring in: buildings and the street network are what 0.1.0 needs
    osmp.buildings = True
    osmp.highways = True
    osmp.water = False
    osmp.forests = False
    osmp.vegetation = False
    osmp.railways = False
    osmp.singleObject = False
    osmp.levelHeight = args.level_height
    # default: stay offline and reproducible from the committed extract alone
    osmp.loadMissingMembers = bool(args.load_missing)
    osmp.subdivide = not args.no_subdivide
    osmp.subdivisionSize = 1.0

    print(f"importing {args.osm} bounds {args.bounds}")
    bpy.ops.blosm.import_data()
    print("import finished")


def proj_batch(points, proj_bin):
    """WGS84 (lat, lon) -> EPSG:26911 (E, N) in one PROJ call."""
    import subprocess

    app_files = os.path.dirname(os.path.dirname(proj_bin))
    env = dict(os.environ)
    env["LD_LIBRARY_PATH"] = f"{app_files}/lib:{env.get('LD_LIBRARY_PATH', '')}"
    payload = "\n".join(f"{lon:.9f} {lat:.9f}" for lat, lon in points) + "\n"
    out = subprocess.run(
        [proj_bin, "+proj=utm", "+zone=11", "+datum=NAD83", "-f", "%.5f"],
        input=payload, capture_output=True, text=True, env=env, check=True,
    ).stdout.strip().splitlines()
    return [tuple(float(v) for v in line.split()) for line in out]


def find_proj():
    import glob
    import shutil

    env = os.environ.get("PROJ_BIN")
    if env and os.path.exists(env):
        return env
    which = shutil.which("proj")
    if which:
        return which
    pattern = os.path.expanduser(
        "~/.local/share/flatpak/app/org.qgis.qgis/x86_64/stable/*/files/bin/proj"
    )
    matches = sorted(glob.glob(pattern))
    if not matches:
        raise SystemExit("no PROJ binary; set PROJ_BIN")
    return matches[0]


def correct_frame(args):
    """Map Blosm's spherical Transverse Mercator onto true EPSG:26911.

    Blosm projects with a sphere of radius 6378137 about the origin, which is a
    conformal but ~0.3 % oversized frame relative to ellipsoidal UTM. The
    difference between the two over a 500 m block is an affine map to well under
    a centimetre, so fit a similarity transform on four spread control points and
    apply it to the imported objects. After this, the blockout uses the same
    metric frame as the DEM and the origin documented in Geo/manifests.
    """
    import math

    import bpy
    from mathutils import Matrix
    from util.transverse_mercator import TransverseMercator

    lat, lon = (float(v) for v in args.origin.split(","))
    tm = TransverseMercator(lat=lat, lon=lon)

    south, west, north, east = (float(v) for v in args.bounds.split(","))
    corners = [(south, west), (south, east), (north, west), (north, east),
               ((south + north) / 2, (west + east) / 2)]
    blosm_xy = [tm.fromGeographic(la, lo)[:2] for la, lo in corners]
    utm_xy = proj_batch(corners, find_proj())

    # full affine fit (6 parameters): absorbs the small shear the spherical
    # projection introduces across the block, which a similarity leaves as ~30 cm
    n = len(blosm_xy)
    bxy = [p for p in blosm_xy]
    uxy = [q for q in utm_xy]
    # solve [x, y, 1] * [a, b, c]^T = u  for both components, least squares
    def solve(target):
        ata = [[0.0] * 3 for _ in range(3)]
        atb = [0.0] * 3
        for (x, y), t in zip(bxy, target):
            row = (x, y, 1.0)
            for i in range(3):
                for j in range(3):
                    ata[i][j] += row[i] * row[j]
                atb[i] += row[i] * t
        # Gaussian elimination
        m = [ata[i][:] + [atb[i]] for i in range(3)]
        for i in range(3):
            piv = max(range(i, 3), key=lambda r: abs(m[r][i]))
            m[i], m[piv] = m[piv], m[i]
            for r in range(3):
                if r != i:
                    f = m[r][i] / m[i][i]
                    for c in range(i, 4):
                        m[r][c] -= f * m[i][c]
        return [m[i][3] / m[i][i] for i in range(3)]

    a, b, tx = solve([q[0] for q in uxy])
    d, e, ty = solve([q[1] for q in uxy])

    residuals = [
        math.hypot(a * p[0] + b * p[1] + tx - q[0], d * p[0] + e * p[1] + ty - q[1])
        for p, q in zip(blosm_xy, utm_xy)
    ]
    print(f"frame correction affine: X = {a:.6f}*x {b:+.6f}*y {tx:+.2f}, "
          f"Y = {d:+.6f}*x {e:.6f}*y {ty:+.2f}")
    print(f"frame correction residual: max {max(residuals)*100:.1f} cm")

    # local origin for the UTM frame: the project origin
    origin_lat_lon = (lat, lon)
    origin_xy = proj_batch([origin_lat_lon], find_proj())[0]

    m = Matrix(((a, b, 0.0, tx - origin_xy[0]), (d, e, 0.0, ty - origin_xy[1]),
                (0.0, 0.0, 1.0, 0.0), (0.0, 0.0, 0.0, 1.0)))
    # only transform roots: children inherit their parent's transform, and
    # transforming both would apply the correction twice
    roots = [o for o in bpy.data.objects if o.parent is None]
    for obj in roots:
        obj.matrix_world = m @ obj.matrix_world
    print(f"frame correction applied to {len(roots)} root objects "
          f"({len(bpy.data.objects) - len(roots)} children inherit)")


def report():
    """Per-object measurements in metres, so a scale bug is obvious in the log."""
    collections = {c.name: c for c in bpy.data.collections}
    meshes = [o for o in bpy.data.objects if o.type == "MESH"]
    print(f"objects: {len(bpy.data.objects)} total, {len(meshes)} mesh, "
          f"collections: {sorted(collections)}")
    if not meshes:
        print("WARNING: no mesh objects were created")
        return

    def dimension(obj):
        return max(obj.dimensions)

    meshes.sort(key=dimension, reverse=True)
    print("largest meshes (name, dimensions in metres):")
    for obj in meshes[:12]:
        dx, dy, dz = obj.dimensions
        print(f"  {obj.name[:40]:40} {dx:7.2f} x {dy:7.2f} x {dz:7.2f}")

    widths = [dimension(o) for o in meshes]
    print(f"mesh extent: max {max(widths):.2f} m, min {min(widths):.3f} m, "
          f"median-ish {sorted(widths)[len(widths)//2]:.2f} m")
    if max(widths) > 1000:
        print("WARNING: something is over 1 km across; check units and bounds")
    if max(widths) < 5:
        print("WARNING: nothing larger than 5 m; the filter probably matched nothing")


def main():
    args = parse_args()
    enable_blosm()
    configure_preferences(args)
    import_osm(args)
    if args.origin:
        correct_frame(args)
    report()

    out = os.path.abspath(args.out)
    os.makedirs(os.path.dirname(out), exist_ok=True)
    bpy.ops.wm.save_as_mainfile(filepath=out)
    print(f"saved {out}")


if __name__ == "__main__":
    main()