sightlines.py

"""
Citywide 360° “visibility polygons” per intersection or any point/node, blocked by building footprints.
"""

from __future__ import annotations

import math
import os
from dataclasses import dataclass
from typing import Optional, Literal

import numpy as np
import geopandas as gpd
import pandas as pd
import requests
from shapely.geometry import Point, LineString, Polygon
from shapely.strtree import STRtree

import osmnx as ox


# ----------------------------
# Utilities
# ----------------------------
STATE_ABBR_TO_NAME = {
    "AL": "Alabama", "AK": "Alaska", "AZ": "Arizona", "AR": "Arkansas",
    "CA": "California", "CO": "Colorado", "CT": "Connecticut", "DE": "Delaware",
    "DC": "DistrictofColumbia", "FL": "Florida", "GA": "Georgia", "HI": "Hawaii",
    "ID": "Idaho", "IL": "Illinois", "IN": "Indiana", "IA": "Iowa",
    "KS": "Kansas", "KY": "Kentucky", "LA": "Louisiana", "ME": "Maine",
    "MD": "Maryland", "MA": "Massachusetts", "MI": "Michigan", "MN": "Minnesota",
    "MS": "Mississippi", "MO": "Missouri", "MT": "Montana", "NE": "Nebraska",
    "NV": "Nevada", "NH": "NewHampshire", "NJ": "NewJersey", "NM": "NewMexico",
    "NY": "NewYork", "NC": "NorthCarolina", "ND": "NorthDakota", "OH": "Ohio",
    "OK": "Oklahoma", "OR": "Oregon", "PA": "Pennsylvania", "RI": "RhodeIsland",
    "SC": "SouthCarolina", "SD": "SouthDakota", "TN": "Tennessee", "TX": "Texas",
    "UT": "Utah", "VT": "Vermont", "VA": "Virginia", "WA": "Washington",
    "WV": "WestVirginia", "WI": "Wisconsin", "WY": "Wyoming",
}

OSM_PRESET_FILTERS = {
    "walk_basic": (
        '["highway"~"primary|primary_link|secondary|secondary_link|'
        'tertiary|tertiary_link|residential|path|living_street|pedestrian"]'
    ),
    # Reasonable “drive-ish” variant (still excludes motorways by default)
    "drive_basic": (
        '["highway"~"primary|primary_link|secondary|secondary_link|'
        'tertiary|tertiary_link|residential|living_street|unclassified"]'
    ),
    # Include service roads too
    "drive_plus_service": (
        '["highway"~"primary|primary_link|secondary|secondary_link|'
        'tertiary|tertiary_link|residential|living_street|unclassified|service"]'
    ),
}

def _place_string(city: str, state: str | None = None, country: str | None = "USA") -> str:
    parts = [p for p in [city, state, country] if p and str(p).strip()]
    return ", ".join(parts)

def _county_place_string(county: str, state: str | None = None, country: str | None = "USA") -> str:
    parts = [p for p in [f"{county} County", state, country] if p and str(p).strip()]
    return ", ".join(parts)

def _download_file(url: str, dst: str, chunk_size: int = 1024 * 1024, verbose: bool = True) -> None:
    os.makedirs(os.path.dirname(dst) or ".", exist_ok=True)
    with requests.get(url, stream=True, timeout=300) as r:
        r.raise_for_status()
        total = int(r.headers.get("Content-Length", "0") or 0)
        done = 0
        with open(dst, "wb") as f:
            for chunk in r.iter_content(chunk_size=chunk_size):
                if not chunk:
                    continue
                f.write(chunk)
                done += len(chunk)
                if verbose and total > 0:
                    pct = 100.0 * done / total
                    print(f"\r  downloading… {pct:5.1f}%", end="")
    if verbose:
        if total > 0:
            print("\r  download complete.")
        else:
            print("  download complete.")


# ----------------------------
# Configuration
# ----------------------------
@dataclass
class SightlineConfig:
    target_crs: str = "EPSG:2285"  # Seattle-area feet (WA North ftUS)
    view_dist: float = 200.0       # feet if target_crs is feet
    step_deg: float = 3.0          # 120 rays per intersection at 3 degrees
    consolidate_intersections: bool = True
    consolidate_tolerance: float = 40.0  # feet in target_crs
    # OSMnx download controls
    network_type: str = "drive"    # ignored if custom_filter is used
    osm_preset: Optional[str] = None         # e.g., "walk_basic"
    osm_custom_filter: Optional[str] = None  # overrides preset if provided
    min_area: float = 5.0          # ft^2: filters tiny/degenerate results
    drop_if_too_small: bool = True
    verbose: bool = True


# ----------------------------
# Main class
# ----------------------------
class CitySightlines:
    def __init__(
        self,
        city: str,
        state: Optional[str] = None,
        country: str = "USA",
        cfg: SightlineConfig = SightlineConfig(),
        county: Optional[str] = None,
    ):
        self.city = city
        self.state = state
        self.country = country
        self.cfg = cfg
        self.county = county

        self.place = _place_string(city, state, country)
        self.county_place = _county_place_string(county, state, country) if county else None

        self.aoi_city: Optional[gpd.GeoDataFrame] = None
        self.aoi_county: Optional[gpd.GeoDataFrame] = None
        self.aoi: Optional[gpd.GeoDataFrame] = None
        self.intersections: Optional[gpd.GeoDataFrame] = None
        self.buildings: Optional[gpd.GeoDataFrame] = None
        self.osm_nodes_full: Optional[gpd.GeoDataFrame] = None
        self.osm_edges_full: Optional[gpd.GeoDataFrame] = None

    # ----------------------------
    # AOI
    # ----------------------------
    def fetch_aoi(self) -> gpd.GeoDataFrame:
        if self.cfg.verbose:
            print(f"[AOI] Geocoding city boundary: {self.place}")
        aoi_city = ox.geocode_to_gdf(self.place).to_crs(self.cfg.target_crs)
        self.aoi_city = aoi_city

        if self.county:
            if self.cfg.verbose:
                print(f"[AOI] Geocoding county boundary: {self.county_place}")
            aoi_county = ox.geocode_to_gdf(self.county_place).to_crs(self.cfg.target_crs)
            self.aoi_county = aoi_county

            inter = aoi_city.unary_union.intersection(aoi_county.unary_union)
            self.aoi = gpd.GeoDataFrame({"geometry": [inter]}, crs=self.cfg.target_crs)
            if self.cfg.verbose:
                print("[AOI] Using AOI = city ∩ county")
        else:
            self.aoi = aoi_city
            if self.cfg.verbose:
                print("[AOI] Using AOI = city boundary")

        return self.aoi

    # ----------------------------
    # OSM intersections (with custom_filter)
    # ----------------------------
    def _resolve_osm_custom_filter(self) -> Optional[str]:
        """
        Mirror your logic:
        - If osm_custom_filter provided, use it
        - Else if osm_preset provided, use that preset
        - Else None => use network_type
        """
        custom_filter = None
        if self.cfg.osm_custom_filter is not None:
            custom_filter = self.cfg.osm_custom_filter
        elif self.cfg.osm_preset is not None:
            if self.cfg.osm_preset not in OSM_PRESET_FILTERS:
                raise ValueError(
                    f"Unknown osm_preset='{self.cfg.osm_preset}'. "
                    f"Valid: {sorted(OSM_PRESET_FILTERS.keys())}"
                )
            custom_filter = OSM_PRESET_FILTERS[self.cfg.osm_preset]
        return custom_filter

    def fetch_intersections(
        self,
        limit_n: Optional[int] = None,
        limit_mode: Literal["first", "random"] = "first",
        random_seed: int = 42,
    ) -> gpd.GeoDataFrame:
        if self.aoi is None:
            self.fetch_aoi()

        custom_filter = self._resolve_osm_custom_filter()

        if custom_filter is not None:
            if self.cfg.verbose:
                print(f"[OSM] Downloading graph with custom_filter (preset={self.cfg.osm_preset}, override={'yes' if self.cfg.osm_custom_filter else 'no'})")
            G = ox.graph_from_place(self.place, custom_filter=custom_filter, simplify=True)
        else:
            if self.cfg.verbose:
                print(f"[OSM] Downloading graph (network_type={self.cfg.network_type}) for: {self.place}")
            G = ox.graph_from_place(self.place, network_type=self.cfg.network_type, simplify=True)

        # ---- Store FULL nodes & ways (edges)
        nodes_full, edges_full = ox.graph_to_gdfs(G, nodes=True, edges=True)
        # Project full nodes/edges to target CRS
        # (edges_full.crs is usually set; if not, project_graph first)
        try:
            nodes_full = nodes_full.to_crs(self.cfg.target_crs)
            edges_full = edges_full.to_crs(self.cfg.target_crs)
        except Exception:
            Gp_tmp = ox.project_graph(G, to_crs=self.cfg.target_crs)
            nodes_full, edges_full = ox.graph_to_gdfs(Gp_tmp, nodes=True, edges=True)
            nodes_full = nodes_full.to_crs(self.cfg.target_crs)
            edges_full = edges_full.to_crs(self.cfg.target_crs)

        # Clip full nodes/edges to AOI (optional but usually what you want)
        aoi_geom = self.aoi.unary_union
        nodes_full = nodes_full[nodes_full.geometry.intersects(aoi_geom)].copy()
        edges_full = edges_full[edges_full.geometry.intersects(aoi_geom)].copy()

        self.osm_nodes_full = nodes_full.reset_index(drop=False)  # keep original node ids as columns
        self.osm_edges_full = edges_full.reset_index(drop=False)  # keep u/v/key as columns

        if self.cfg.verbose:
            print(f"[OSM] Full nodes (AOI-clipped): {len(self.osm_nodes_full):,}")
            print(f"[OSM] Full ways/edges (AOI-clipped): {len(self.osm_edges_full):,}")

        # ---- Build the INTERSECTIONS layer used for sightlines 
        if self.cfg.consolidate_intersections:
            if self.cfg.verbose:
                print(f"[OSM] Consolidating intersections (tolerance={self.cfg.consolidate_tolerance} in {self.cfg.target_crs})")
            Gp = ox.project_graph(G, to_crs=self.cfg.target_crs)
            Gc = ox.consolidate_intersections(
                Gp,
                tolerance=self.cfg.consolidate_tolerance,
                rebuild_graph=True,
                dead_ends=False,
            )
            nodes_int, _ = ox.graph_to_gdfs(Gc, nodes=True, edges=True)
            inter = gpd.GeoDataFrame(geometry=nodes_int.geometry, crs=self.cfg.target_crs).reset_index(drop=True)
        else:
            # Use full nodes as base, optionally filter to "true intersections"
            nodes_tmp = self.osm_nodes_full.copy()
            if "street_count" in nodes_tmp.columns:
                nodes_tmp = nodes_tmp[nodes_tmp["street_count"] >= 3].copy()
            inter = gpd.GeoDataFrame(geometry=nodes_tmp.geometry, crs=self.cfg.target_crs).reset_index(drop=True)

        # Clip to AOI
        inter = inter[inter.geometry.intersects(aoi_geom)].copy().reset_index(drop=True)

        if self.cfg.verbose:
            print(f"[OSM] Intersections after AOI clip: {len(inter):,}")

        # Testing limiter
        if limit_n is not None and limit_n < len(inter):
            if limit_mode == "first":
                inter = inter.iloc[:limit_n].copy()
            else:
                inter = inter.sample(n=limit_n, random_state=random_seed).copy()
            inter = inter.reset_index(drop=True)
            if self.cfg.verbose:
                print(f"[OSM] Intersections after limit ({limit_mode} {limit_n:,}): {len(inter):,}")

        self.intersections = inter
        return inter

    # ----------------------------
    # Microsoft US Buildings (state ZIP)
    # ----------------------------


    def load_buildings_if_exists(
        self,
        path: str,
        layer: str | None = None,
        target_crs: str | None = None,
    ) -> bool:
        """
        If a buildings file exists locally, load it into self.buildings and return True.
        Supports GPKG (with optional layer), GeoJSON, Shapefile, etc.

        This is the fastest path because you avoid reading the WA-wide GeoJSON zip again.
        """
        if not os.path.exists(path):
            if self.cfg.verbose:
                print(f"[MS-US] Buildings cache not found: {path}")
            return False

        if self.cfg.verbose:
            msg = f"[MS-US] Loading cached buildings: {path}"
            if layer:
                msg += f" (layer={layer})"
            print(msg)

        gdf = gpd.read_file(path, layer=layer) if layer else gpd.read_file(path)

        # Ensure CRS
        if target_crs is None:
            target_crs = self.cfg.target_crs

        if gdf.crs is None:
            # assume target CRS if missing
            gdf = gdf.set_crs(target_crs, allow_override=True)
        elif str(gdf.crs) != str(target_crs):
            gdf = gdf.to_crs(target_crs)

        self.buildings = gdf.reset_index(drop=True)

        if self.cfg.verbose:
            print(f"[MS-US] Loaded cached buildings: {len(self.buildings):,}")
        return True


    def fetch_buildings_ms_us_cached(
        self,
        state: str | None = None,
        cache_dir: str = "./_cache_ms_buildings",
        force_download: bool = False,
        clip_to: Literal["aoi", "city", "county", "none"] = "aoi",
        try_bbox_read: bool = True,
        # clipped-cache options
        clipped_cache_path: str | None = None,
        clipped_cache_layer: str | None = "ms_buildings",
        reuse_clipped_cache: bool = True,
        # if we create buildings, optionally write to clipped cache
        write_clipped_cache: bool = True,
    ) -> gpd.GeoDataFrame:
        """
        Wrapper around fetch_buildings_ms_us() that:
        1) Reuses a pre-clipped cache file (gpkg/geojson/etc) if it exists (fastest)
        2) Otherwise downloads/reads the state ZIP (zip is cached), clips, and optionally writes the clipped cache.

        """
        # 1) Try clipped cache first
        if reuse_clipped_cache and clipped_cache_path:
            loaded = self.load_buildings_if_exists(
                clipped_cache_path,
                layer=clipped_cache_layer,
                target_crs=self.cfg.target_crs,
            )
            if loaded:
                return self.buildings

        # 2) Otherwise build it (this uses your existing fetch_buildings_ms_us)
        bld = self.fetch_buildings_ms_us(
            state=state,
            cache_dir=cache_dir,
            force=force_download,
            clip_to=clip_to,
            try_bbox_read=try_bbox_read,
        )

        # 3) Optionally write clipped cache for next time
        if write_clipped_cache and clipped_cache_path:
            os.makedirs(os.path.dirname(clipped_cache_path) or ".", exist_ok=True)
            ext = os.path.splitext(clipped_cache_path)[1].lower()

            if self.cfg.verbose:
                print(f"[MS-US] Writing clipped buildings cache: {clipped_cache_path}")

            if ext in [".gpkg", ".geopackage"]:
                bld.to_file(clipped_cache_path, layer=(clipped_cache_layer or "ms_buildings"), driver="GPKG", mode="w")
            elif ext in [".geojson", ".json"]:
                bld.to_file(clipped_cache_path, driver="GeoJSON")
            else:
                # default to gpkg
                bld.to_file(clipped_cache_path, layer=(clipped_cache_layer or "ms_buildings"), driver="GPKG", mode="w")

        return bld
        
    def fetch_buildings_ms_us(
        self,
        state: Optional[str] = None,
        cache_dir: str = "./_cache_ms_buildings",
        force: bool = False,
        clip_to: Literal["aoi", "city", "county", "none"] = "aoi",
        try_bbox_read: bool = True,
    ) -> gpd.GeoDataFrame:
        if self.aoi is None:
            self.fetch_aoi()

        if state is None:
            state = self.state
        if state is None:
            raise ValueError("State is required for US-only MS buildings (e.g., 'WA').")

        state_key = state.strip()
        if len(state_key) == 2:
            state_key = state_key.upper()
            if state_key not in STATE_ABBR_TO_NAME:
                raise ValueError(f"Unknown state abbreviation: {state_key}")
            state_name = STATE_ABBR_TO_NAME[state_key]
        else:
            state_name = state_key.replace(" ", "")

        base = "https://minedbuildings.z5.web.core.windows.net/legacy/usbuildings-v2"
        url = f"{base}/{state_name}.geojson.zip"

        zip_path = os.path.join(cache_dir, f"{state_name}.geojson.zip")
        if force or (not os.path.exists(zip_path)):
            if self.cfg.verbose:
                print(f"[MS-US] Downloading: {url}")
            _download_file(url, zip_path, verbose=self.cfg.verbose)
        else:
            if self.cfg.verbose:
                print(f"[MS-US] Using cached: {zip_path}")

        # Choose clip geometry
        if clip_to == "none":
            clip_geom_wgs = None
        elif clip_to == "city":
            if self.aoi_city is None:
                self.fetch_aoi()
            clip_geom_wgs = self.aoi_city.to_crs("EPSG:4326").unary_union
        elif clip_to == "county":
            if not self.county:
                raise ValueError("clip_to='county' requires county=... in CitySightlines().")
            if self.aoi_county is None:
                self.fetch_aoi()
            clip_geom_wgs = self.aoi_county.to_crs("EPSG:4326").unary_union
        else:
            clip_geom_wgs = self.aoi.to_crs("EPSG:4326").unary_union

        gdf = None
        if try_bbox_read and clip_geom_wgs is not None:
            bbox = gpd.GeoSeries([clip_geom_wgs], crs="EPSG:4326").total_bounds
            try:
                if self.cfg.verbose:
                    print("[MS-US] Attempting bbox read…")
                gdf = gpd.read_file(zip_path, bbox=tuple(bbox))
            except Exception:
                gdf = None

        if gdf is None:
            if self.cfg.verbose:
                print("[MS-US] Reading full state file (then clipping)…")
            gdf = gpd.read_file(zip_path)

        if gdf.crs is None:
            gdf = gdf.set_crs("EPSG:4326", allow_override=True)
        else:
            gdf = gdf.to_crs("EPSG:4326")

        if clip_geom_wgs is not None:
            gdf = gdf[gdf.geometry.notna() & gdf.geometry.is_valid].copy()
            gdf = gdf[gdf.geometry.intersects(clip_geom_wgs)].copy()

        gdf = gdf.to_crs(self.cfg.target_crs).reset_index(drop=True)

        self.buildings = gdf
        if self.cfg.verbose:
            print(f"[MS-US] Buildings loaded: {len(gdf):,}")
        return gdf

    # ----------------------------
    # Visibility polygon builder
    # ----------------------------
    def build_visibility_polygons(
        self,
        intersection_id_col: Optional[str] = None,
        chunk_size: int = 2500,
    ) -> gpd.GeoDataFrame:
        if self.intersections is None:
            raise ValueError("Intersections not loaded. Run fetch_intersections() first.")
        if self.buildings is None:
            raise ValueError("Buildings not loaded. Run fetch_buildings_ms_us() first.")

        inter = self.intersections
        bld = self.buildings

        if inter.crs != bld.crs:
            raise ValueError("CRS mismatch between intersections and buildings.")

        if self.cfg.verbose:
            print(f"[VIS] Building STRtree for {len(bld):,} buildings...")

        b_geoms = list(bld.geometry)
        b_tree = STRtree(b_geoms)

        def _as_geoms(tree_results):
            if len(tree_results) == 0:
                return []
            first = tree_results[0]
            if isinstance(first, (int, np.integer)):
                return [b_geoms[i] for i in tree_results]
            return list(tree_results)

        angles = np.arange(0.0, 360.0, self.cfg.step_deg, dtype=float)

        records = []
        n = len(inter)

        for start in range(0, n, chunk_size):
            end = min(start + chunk_size, n)
            if self.cfg.verbose:
                print(f"[VIS] Processing intersections {start:,}–{end-1:,} of {n:,}...")

            for i, row in inter.iloc[start:end].iterrows():
                geom = row.geometry
                if geom is None or geom.is_empty:
                    continue

                pt = geom if geom.geom_type == "Point" else geom.centroid
                x0, y0 = pt.x, pt.y
                inter_id = row[intersection_id_col] if intersection_id_col else i

                candidates = b_tree.query(pt.buffer(self.cfg.view_dist).envelope)
                candidates = _as_geoms(candidates)

                boundary_pts = []
                origin = pt

                for deg in angles:
                    x1, y1 = self._ray_endpoint_xy(x0, y0, self.cfg.view_dist, float(deg))
                    ray = LineString([(x0, y0), (x1, y1)])

                    nearest_hit = None
                    nearest_d = self.cfg.view_dist

                    for bg in candidates:
                        if bg is None or bg.is_empty:
                            continue
                        if not ray.intersects(bg):
                            continue

                        inter_geom = ray.intersection(bg)
                        hit_pt = self._nearest_hit_point(inter_geom, origin)
                        if hit_pt is None:
                            continue

                        d = origin.distance(hit_pt)
                        if d < 1e-6:
                            continue

                        if d < nearest_d:
                            nearest_d = d
                            nearest_hit = hit_pt

                    boundary_pts.append((x1, y1) if nearest_hit is None else (nearest_hit.x, nearest_hit.y))

                poly = Polygon(boundary_pts)
                if not poly.is_valid:
                    poly = poly.buffer(0)

                if poly.is_empty:
                    continue

                if poly.geom_type == "MultiPolygon":
                    poly = max(list(poly.geoms), key=lambda p: p.area)

                if self.cfg.drop_if_too_small and poly.area < self.cfg.min_area:
                    continue

                records.append(
                    {
                        "intersection_id": inter_id,
                        "view_dist": float(self.cfg.view_dist),
                        "step_deg": float(self.cfg.step_deg),
                        "geometry": poly,
                    }
                )

        out = gpd.GeoDataFrame(records, crs=inter.crs)
        if self.cfg.verbose:
            print(f"[VIS] Output polygons: {len(out):,}")
        return out


    # ----------------------------
    # Export
    # ----------------------------
    def export_all(
        self,
        vis_gdf: gpd.GeoDataFrame,
        path: str,
        layer_vis: str = "visibility",
        layer_intersections: str = "intersections",
        layer_buildings: str = "ms_buildings",
        layer_osm_nodes: str = "osm_nodes_full",
        layer_osm_ways: str = "osm_ways_full",
        buildings_simplify_tolerance: float | None = None,
        overwrite_layers: bool = True,
    ) -> None:
        """
        Export to ONE GeoPackage:
        - visibility polygons (vis_gdf)
        - intersections used for analysis (self.intersections)
        - Microsoft buildings (self.buildings)
        - full OSM nodes (self.osm_nodes_full)
        - full OSM ways/edges (self.osm_edges_full)
        """
        if self.cfg.verbose:
            print(f"[OUT] Writing bundle to: {path}")

        if self.intersections is None:
            raise ValueError("intersections is None. Run fetch_intersections() first.")
        if self.buildings is None:
            raise ValueError("buildings is None. Run fetch_buildings_ms_us() first.")
        if self.osm_nodes_full is None or self.osm_edges_full is None:
            raise ValueError("OSM full nodes/ways not available. Run fetch_intersections() first (it populates them).")

        os.makedirs(os.path.dirname(path) or ".", exist_ok=True)

        # Optionally simplify buildings for output size
        bld_out = self.buildings
        if buildings_simplify_tolerance is not None and buildings_simplify_tolerance > 0:
            if self.cfg.verbose:
                print(f"[OUT] Simplifying buildings (tolerance={buildings_simplify_tolerance}) for export...")
            bld_out = bld_out.copy()
            bld_out["geometry"] = bld_out.geometry.simplify(buildings_simplify_tolerance, preserve_topology=True)

        mode_first = "w" if overwrite_layers else "a"

        # Write layers (first write uses mode_first, rest append)
        vis_gdf.to_file(path, layer=layer_vis, driver="GPKG", mode=mode_first)
        self.intersections.to_file(path, layer=layer_intersections, driver="GPKG", mode="a")
        bld_out.to_file(path, layer=layer_buildings, driver="GPKG", mode="a")
        self.osm_nodes_full.to_file(path, layer=layer_osm_nodes, driver="GPKG", mode="a")
        self.osm_edges_full.to_file(path, layer=layer_osm_ways, driver="GPKG", mode="a")

        if self.cfg.verbose:
            print("[OUT] Done.")


    def export_split(
        self,
        vis_gdf: gpd.GeoDataFrame,
        out_dir: str,
        vis_name: str = "visibility.gpkg",
        nodes_name: str = "osm_nodes.gpkg",
        buildings_name: str = "ms_buildings.gpkg",
        buildings_simplify_tolerance: float | None = None,
    ) -> None:
        """
        Export as THREE separate GeoPackages (sometimes nicer for quick debugging).
        """
        os.makedirs(out_dir, exist_ok=True)

        if self.intersections is None:
            raise ValueError("intersections is None. Run fetch_intersections() first.")
        if self.buildings is None:
            raise ValueError("buildings is None. Run fetch_buildings_ms_us() first.")

        vis_path = os.path.join(out_dir, vis_name)
        nodes_path = os.path.join(out_dir, nodes_name)
        bld_path = os.path.join(out_dir, buildings_name)

        if self.cfg.verbose:
            print(f"[OUT] Writing: {vis_path}")
        vis_gdf.to_file(vis_path, layer="visibility", driver="GPKG", mode="w")

        if self.cfg.verbose:
            print(f"[OUT] Writing: {nodes_path}")
        self.intersections.to_file(nodes_path, layer="osm_nodes", driver="GPKG", mode="w")

        bld_out = self.buildings
        if buildings_simplify_tolerance is not None and buildings_simplify_tolerance > 0:
            if self.cfg.verbose:
                print(f"[OUT] Simplifying buildings (tolerance={buildings_simplify_tolerance}) for export...")
            bld_out = bld_out.copy()
            bld_out["geometry"] = bld_out.geometry.simplify(buildings_simplify_tolerance, preserve_topology=True)

        if self.cfg.verbose:
            print(f"[OUT] Writing: {bld_path}")
        bld_out.to_file(bld_path, layer="ms_buildings", driver="GPKG", mode="w")

        if self.cfg.verbose:
            print("[OUT] Done.")

    # ----------------------------
    # Geometry helpers
    # ----------------------------
    @staticmethod
    def _unit_vector_from_deg(deg: float):
        rad = math.radians(deg)
        return math.cos(rad), math.sin(rad)

    def _ray_endpoint_xy(self, x: float, y: float, dist: float, deg: float):
        ux, uy = self._unit_vector_from_deg(deg)
        return (x + ux * dist, y + uy * dist)

    @staticmethod
    def _nearest_hit_point(intersection_geom, origin: Point):
        if intersection_geom.is_empty:
            return None

        gt = intersection_geom.geom_type

        if gt == "Point":
            return intersection_geom

        if gt == "MultiPoint":
            pts = list(intersection_geom.geoms)
            return min(pts, key=lambda p: origin.distance(p)) if pts else None

        if gt in ("LineString", "LinearRing"):
            return intersection_geom.interpolate(intersection_geom.project(origin))

        if gt == "MultiLineString":
            candidates = [ln.interpolate(ln.project(origin)) for ln in intersection_geom.geoms]
            return min(candidates, key=lambda p: origin.distance(p)) if candidates else None

        if gt == "GeometryCollection":
            pts = []
            for g in intersection_geom.geoms:
                p = CitySightlines._nearest_hit_point(g, origin)
                if p is not None:
                    pts.append(p)
            return min(pts, key=lambda p: origin.distance(p)) if pts else None

        if gt == "Polygon":
            b = intersection_geom.boundary
            return b.interpolate(b.project(origin))

        if gt == "MultiPolygon":
            candidates = []
            for poly in intersection_geom.geoms:
                b = poly.boundary
                candidates.append(b.interpolate(b.project(origin)))
            return min(candidates, key=lambda p: origin.distance(p)) if candidates else None

        return None


# ----------------------------
# Example usage (Seattle)
# ----------------------------

cfg = SightlineConfig(
    target_crs="EPSG:2285",
    view_dist=500.0,
    step_deg=5.0,                    # faster for first run
    consolidate_intersections=True,
    consolidate_tolerance=40.0,
    # Choose one:
    # 1) Preset
    osm_preset="walk_basic",
    # 2) Or override with your own:
    # osm_custom_filter='["highway"~"residential|tertiary|secondary"]',
    # If no preset/custom_filter, network_type is used instead
    network_type="drive",
    verbose=True,
)

cs = CitySightlines(city="Seattle", state="WA", country="USA", county="King", cfg=cfg)

cs.fetch_aoi()

# Quick test: sample 500 intersections
cs.fetch_intersections(limit_n=500, limit_mode="random", random_seed=7)

# cs.fetch_buildings_ms_us(state="WA", cache_dir="./_cache_ms", clip_to="aoi", try_bbox_read=True)
cs.fetch_buildings_ms_us_cached(
    state="WA",
    cache_dir="./_cache_ms",                    # stores Washington.geojson.zip
    clip_to="aoi",
    clipped_cache_path="./_cache_ms/seattle_msbuildings_aoi.gpkg",  # stores *clipped* buildings
    clipped_cache_layer="ms_buildings",
    reuse_clipped_cache=True,
    write_clipped_cache=True,
)

vis = cs.build_visibility_polygons(chunk_size=1000)

cs.export_all(vis, 
    "seattle_bundle.gpkg",
    layer_vis="vis_poly",
    layer_intersections="intersections_used",
    layer_buildings="ms_buildings",
    layer_osm_nodes="osm_nodes",
    layer_osm_ways="osm_ways",
    buildings_simplify_tolerance=1.0
)

# cs.__dict__