Dataset Overview¶

Notebook that:

1. Downloads ml-data artifacts from S3 (split IDs, class_weights.json, optional OSM counts)
2. Computes dataset stats: unique HUCs, train/val/test counts, OSM bridge counts, total points
3. Produces class distribution horizontal bar chart (from class_weights.json)
4. Optionally produces per-bridge point count histogram (requires local SILVER_NORMALIZED_DIR)
5. Downloads HUC8 boundaries from S3 (WBDHU8_webproj.gpkg) and produces a map of HUC8 distribution

Configure S3 bucket/prefix and optional AWS profile below; all data is read from S3.

Config¶

Set S3 bucket/prefix for ml-data, optional AWS profile, public HUC8 boundaries URL or local path, and cache directory.

In [1]:

import os

# S3: ml-data prefix (split files, class_weights.json, optional osm_bridge_counts.json)
S3_BUCKET = os.environ.get("BRIDGE_S3_BUCKET", "fimc-data")
S3_ML_PREFIX = os.environ.get("BRIDGE_S3_ML_PREFIX", "bridge-classification/ml-data")
# Optional AWS profile (e.g. "Data"); leave None for default credentials
# AWS_PROFILE = os.environ.get("AWS_PROFILE", None)
AWS_PROFILE = os.environ.get("AWS_PROFILE", "Data")

# Local cache for downloaded files (run notebook from repo root or from notebooks/)
try:
    NOTEBOOK_DIR = os.path.dirname(os.path.abspath(__file__))
except NameError:
    NOTEBOOK_DIR = os.getcwd()
if not os.path.basename(NOTEBOOK_DIR) == "notebooks":
    NOTEBOOK_DIR = os.path.join(NOTEBOOK_DIR, "notebooks")
CACHE_DIR = os.path.join(NOTEBOOK_DIR, "data")
OUTPUT_DIR = os.path.join(NOTEBOOK_DIR, "outputs")
os.makedirs(CACHE_DIR, exist_ok=True)
os.makedirs(OUTPUT_DIR, exist_ok=True)

# HUC8 boundaries: S3 key (downloaded to cache and used for map)
S3_HUC8_KEY = os.environ.get("BRIDGE_S3_HUC8_KEY", "bridge-classification/hucs/WBDHU8_webproj.gpkg")
# OSM gpkg prefix for fallback OSM counts (when osm_bridge_counts.json is not on S3)
S3_OSM_HUCS_PREFIX = os.environ.get("BRIDGE_S3_OSM_HUCS_PREFIX", "bridge-classification/osm/hucs")

# Optional: local path to silver_training_normalized (per-bridge JSONs) for point-count histogram. Set to None to skip.
SILVER_NORMALIZED_DIR = "data/ml-data/silver_training_normalized"  # e.g. None, Path("data/ml-data/silver_training_normalized")

# Class names and colors for class distribution bar chart (same order as class_weights total_counts 0-3)
CLASS_NAMES = ["Background", "Ground/Water", "Bridge Deck", "Obstacles"]
CLASS_COLORS = ["#9E9E9E", "#8D6E63", "#2196F3", "#FF5722"]

import os # S3: ml-data prefix (split files, class_weights.json, optional osm_bridge_counts.json) S3_BUCKET = os.environ.get("BRIDGE_S3_BUCKET", "fimc-data") S3_ML_PREFIX = os.environ.get("BRIDGE_S3_ML_PREFIX", "bridge-classification/ml-data") # Optional AWS profile (e.g. "Data"); leave None for default credentials # AWS_PROFILE = os.environ.get("AWS_PROFILE", None) AWS_PROFILE = os.environ.get("AWS_PROFILE", "Data") # Local cache for downloaded files (run notebook from repo root or from notebooks/) try: NOTEBOOK_DIR = os.path.dirname(os.path.abspath(__file__)) except NameError: NOTEBOOK_DIR = os.getcwd() if not os.path.basename(NOTEBOOK_DIR) == "notebooks": NOTEBOOK_DIR = os.path.join(NOTEBOOK_DIR, "notebooks") CACHE_DIR = os.path.join(NOTEBOOK_DIR, "data") OUTPUT_DIR = os.path.join(NOTEBOOK_DIR, "outputs") os.makedirs(CACHE_DIR, exist_ok=True) os.makedirs(OUTPUT_DIR, exist_ok=True) # HUC8 boundaries: S3 key (downloaded to cache and used for map) S3_HUC8_KEY = os.environ.get("BRIDGE_S3_HUC8_KEY", "bridge-classification/hucs/WBDHU8_webproj.gpkg") # OSM gpkg prefix for fallback OSM counts (when osm_bridge_counts.json is not on S3) S3_OSM_HUCS_PREFIX = os.environ.get("BRIDGE_S3_OSM_HUCS_PREFIX", "bridge-classification/osm/hucs") # Optional: local path to silver_training_normalized (per-bridge JSONs) for point-count histogram. Set to None to skip. SILVER_NORMALIZED_DIR = "data/ml-data/silver_training_normalized" # e.g. None, Path("data/ml-data/silver_training_normalized") # Class names and colors for class distribution bar chart (same order as class_weights total_counts 0-3) CLASS_NAMES = ["Background", "Ground/Water", "Bridge Deck", "Obstacles"] CLASS_COLORS = ["#9E9E9E", "#8D6E63", "#2196F3", "#FF5722"]

In [2]:

# Run after Config. All imports used elsewhere in the notebook.
import os
from pathlib import Path
import json
import geopandas as gpd
import boto3
import numpy as np
from concurrent.futures import ThreadPoolExecutor, as_completed
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from matplotlib.patches import Patch
from typing import Optional

# Run after Config. All imports used elsewhere in the notebook. import os from pathlib import Path import json import geopandas as gpd import boto3 import numpy as np from concurrent.futures import ThreadPoolExecutor, as_completed import matplotlib.pyplot as plt import matplotlib.ticker as ticker from matplotlib.patches import Patch from typing import Optional

Download from S3¶

Download split ID files and class_weights.json from the ml-data prefix.

In [ ]:

def get_s3_client():
    if AWS_PROFILE:
        session = boto3.Session(profile_name=AWS_PROFILE)
    else:
        session = boto3.Session()
    return session.client("s3")

s3 = get_s3_client()

ML_FILES = [
    "split_train_ids.txt",
    "split_val_ids.txt",
    "split_test_ids.txt",
    "class_weights.json",
]

for key in ML_FILES:
    s3_key = f"{S3_ML_PREFIX.rstrip('/')}/{key}"
    local_path = Path(CACHE_DIR) / key
    try:
        s3.download_file(S3_BUCKET, s3_key, str(local_path))
        print(f"Downloaded {key}")
    except Exception as e:
        print(f"Failed {key}: {e}")

def get_s3_client(): if AWS_PROFILE: session = boto3.Session(profile_name=AWS_PROFILE) else: session = boto3.Session() return session.client("s3") s3 = get_s3_client() ML_FILES = [ "split_train_ids.txt", "split_val_ids.txt", "split_test_ids.txt", "class_weights.json", ] for key in ML_FILES: s3_key = f"{S3_ML_PREFIX.rstrip('/')}/{key}" local_path = Path(CACHE_DIR) / key try: s3.download_file(S3_BUCKET, s3_key, str(local_path)) print(f"Downloaded {key}") except Exception as e: print(f"Failed {key}: {e}")

Downloaded split_train_ids.txt
Downloaded split_val_ids.txt
Downloaded split_test_ids.txt
Downloaded class_weights.json

Stats 1 – Split (line counts and unique HUCs)¶

Read each split file; report line counts and unique HUCs (first segment of each line huc_id/bridge_stem).

In [ ]:

cache = Path(CACHE_DIR)
train_path = cache / "split_train_ids.txt"
val_path = cache / "split_val_ids.txt"
test_path = cache / "split_test_ids.txt"

def load_ids(path):
    lines = []
    with open(path) as f:
        for line in f:
            line = line.strip()
            if line:
                lines.append(line)
    return lines

def unique_hucs(lines):
    return set(line.split("/", 1)[0] for line in lines)

train_ids = load_ids(train_path)
val_ids = load_ids(val_path)
test_ids = load_ids(test_path)

n_train, n_val, n_test = len(train_ids), len(val_ids), len(test_ids)
n_total = n_train + n_val + n_test

hucs_train = unique_hucs(train_ids)
hucs_val = unique_hucs(val_ids)
hucs_test = unique_hucs(test_ids)
hucs_all = hucs_train | hucs_val | hucs_test

print("Dataset split (line counts):")
print(f"  Train:      {n_train:,}")
print(f"  Validation: {n_val:,}")
print(f"  Test:      {n_test:,}")
print(f"  Total:     {n_total:,}")
print()
print("Unique HUCs in split (union of train/val/test):")
print(f"  {len(hucs_all):,}")

cache = Path(CACHE_DIR) train_path = cache / "split_train_ids.txt" val_path = cache / "split_val_ids.txt" test_path = cache / "split_test_ids.txt" def load_ids(path): lines = [] with open(path) as f: for line in f: line = line.strip() if line: lines.append(line) return lines def unique_hucs(lines): return set(line.split("/", 1)[0] for line in lines) train_ids = load_ids(train_path) val_ids = load_ids(val_path) test_ids = load_ids(test_path) n_train, n_val, n_test = len(train_ids), len(val_ids), len(test_ids) n_total = n_train + n_val + n_test hucs_train = unique_hucs(train_ids) hucs_val = unique_hucs(val_ids) hucs_test = unique_hucs(test_ids) hucs_all = hucs_train | hucs_val | hucs_test print("Dataset split (line counts):") print(f" Train: {n_train:,}") print(f" Validation: {n_val:,}") print(f" Test: {n_test:,}") print(f" Total: {n_total:,}") print() print("Unique HUCs in split (union of train/val/test):") print(f" {len(hucs_all):,}")

Dataset split (line counts):
  Train:      476,980
  Validation: 120,266
  Test:      180
  Total:     597,426

Unique HUCs in split (union of train/val/test):
  1,955

Stats 2 – Total points (from class_weights.json)¶

Sum total_counts in class_weights.json (per-class point counts from training set).

In [ ]:

weights_path = Path(CACHE_DIR) / "class_weights.json"
with open(weights_path) as f:
    cw = json.load(f)

total_counts = cw.get("total_counts", {})
total_points = sum(int(v) for v in total_counts.values())

print("Total points (from class_weights.json):")
print(f"  {total_points:,}")
print("Per-class counts:")
for k in sorted(total_counts.keys(), key=int):
    print(f"  Class {k}: {int(total_counts[k]):,}")

weights_path = Path(CACHE_DIR) / "class_weights.json" with open(weights_path) as f: cw = json.load(f) total_counts = cw.get("total_counts", {}) total_points = sum(int(v) for v in total_counts.values()) print("Total points (from class_weights.json):") print(f" {total_points:,}") print("Per-class counts:") for k in sorted(total_counts.keys(), key=int): print(f" Class {k}: {int(total_counts[k]):,}")

Total points (from class_weights.json):
  3,846,139,057
Per-class counts:
  Class 0: 154,663,102
  Class 1: 645,015,460
  Class 2: 2,636,395,892
  Class 3: 410,064,603

Stats 3 – OSM bridge counts¶

Try to download precomputed osm_bridge_counts.json from S3. If missing, list HUCs under the OSM prefix, download each HUC's gpkg files to notebooks/data, count features with geopandas, and write osm_bridge_counts.json to notebooks/outputs.

In [ ]:

OSM_COUNTS_KEY = f"{S3_ML_PREFIX.rstrip('/')}/osm_bridge_counts.json"
osm_counts_path = Path(CACHE_DIR) / "osm_bridge_counts.json"
osm_output_path = Path(OUTPUT_DIR) / "osm_bridge_counts.json"

all_osm_bridges = None
osm_bridge_with_lidar = None

try:
    s3.download_file(S3_BUCKET, OSM_COUNTS_KEY, str(osm_counts_path))
    with open(osm_counts_path) as f:
        counts = json.load(f)
    all_osm_bridges = counts.get("all_osm_bridges")
    osm_bridge_with_lidar = counts.get("osm_bridge_with_lidar")
    print("Loaded osm_bridge_counts.json from S3.")
except Exception as e:
    print(f"osm_bridge_counts.json not found on S3 ({e}). Computing from OSM gpkg files (parallel)...")
    prefix = S3_OSM_HUCS_PREFIX.rstrip("/") + "/"
    paginator = s3.get_paginator("list_objects_v2")
    huc_prefixes = []
    for page in paginator.paginate(Bucket=S3_BUCKET, Prefix=prefix, Delimiter="/"):
        for p in page.get("CommonPrefixes", []):
            huc_prefixes.append(p["Prefix"].rstrip("/").split("/")[-1])

    def _osm_count_one_huc(args):
        import boto3
        import geopandas as gpd
        from pathlib import Path
        bucket, prefix_arg, huc_id, cache_dir, profile_name = args
        session = boto3.Session(profile_name=profile_name) if profile_name else boto3.Session()
        s3_client = session.client("s3")
        huc_dir = Path(cache_dir) / "osm" / "hucs" / huc_id
        huc_dir.mkdir(parents=True, exist_ok=True)
        n_all, n_lidar = 0, 0
        for name, key_suffix in [
            ("osm_bridges_subset", "osm_bridges_subset__{}.gpkg"),
            ("osm_bridges_lidar_subset", "osm_bridges_lidar_subset__{}.gpkg"),
        ]:
            s3_key = f"{prefix_arg}{huc_id}/{key_suffix.format(huc_id)}"
            local_gpkg = huc_dir / key_suffix.format(huc_id)
            try:
                s3_client.download_file(bucket, s3_key, str(local_gpkg))
                gdf = gpd.read_file(str(local_gpkg))
                n = len(gdf)
                if name == "osm_bridges_subset":
                    n_all = n
                else:
                    n_lidar = n
            except Exception:
                pass
        return (n_all, n_lidar)

    from concurrent.futures import ThreadPoolExecutor, as_completed
    workers = min(32, len(huc_prefixes))
    work_items = [(S3_BUCKET, prefix, huc_id, CACHE_DIR, AWS_PROFILE) for huc_id in huc_prefixes]
    total_all = 0
    total_lidar = 0
    with ThreadPoolExecutor(max_workers=workers) as executor:
        futures = {executor.submit(_osm_count_one_huc, item): item for item in work_items}
        for i, future in enumerate(as_completed(futures)):
            n_all, n_lidar = future.result()
            total_all += n_all
            total_lidar += n_lidar
            if (i + 1) % 200 == 0:
                print(f"  Processed {i + 1}/{len(huc_prefixes)} HUCs...")
    all_osm_bridges = total_all
    osm_bridge_with_lidar = total_lidar
    counts = {"all_osm_bridges": all_osm_bridges, "osm_bridge_with_lidar": osm_bridge_with_lidar}
    with open(osm_output_path, "w") as f:
        json.dump(counts, f, indent=2)
    print(f"Wrote {osm_output_path}")

if all_osm_bridges is not None and osm_bridge_with_lidar is not None:
    print(f"  all_osm_bridges:       {all_osm_bridges:,}")
    print(f"  osm_bridge_with_lidar: {osm_bridge_with_lidar:,}")

OSM_COUNTS_KEY = f"{S3_ML_PREFIX.rstrip('/')}/osm_bridge_counts.json" osm_counts_path = Path(CACHE_DIR) / "osm_bridge_counts.json" osm_output_path = Path(OUTPUT_DIR) / "osm_bridge_counts.json" all_osm_bridges = None osm_bridge_with_lidar = None try: s3.download_file(S3_BUCKET, OSM_COUNTS_KEY, str(osm_counts_path)) with open(osm_counts_path) as f: counts = json.load(f) all_osm_bridges = counts.get("all_osm_bridges") osm_bridge_with_lidar = counts.get("osm_bridge_with_lidar") print("Loaded osm_bridge_counts.json from S3.") except Exception as e: print(f"osm_bridge_counts.json not found on S3 ({e}). Computing from OSM gpkg files (parallel)...") prefix = S3_OSM_HUCS_PREFIX.rstrip("/") + "/" paginator = s3.get_paginator("list_objects_v2") huc_prefixes = [] for page in paginator.paginate(Bucket=S3_BUCKET, Prefix=prefix, Delimiter="/"): for p in page.get("CommonPrefixes", []): huc_prefixes.append(p["Prefix"].rstrip("/").split("/")[-1]) def _osm_count_one_huc(args): import boto3 import geopandas as gpd from pathlib import Path bucket, prefix_arg, huc_id, cache_dir, profile_name = args session = boto3.Session(profile_name=profile_name) if profile_name else boto3.Session() s3_client = session.client("s3") huc_dir = Path(cache_dir) / "osm" / "hucs" / huc_id huc_dir.mkdir(parents=True, exist_ok=True) n_all, n_lidar = 0, 0 for name, key_suffix in [ ("osm_bridges_subset", "osm_bridges_subset__{}.gpkg"), ("osm_bridges_lidar_subset", "osm_bridges_lidar_subset__{}.gpkg"), ]: s3_key = f"{prefix_arg}{huc_id}/{key_suffix.format(huc_id)}" local_gpkg = huc_dir / key_suffix.format(huc_id) try: s3_client.download_file(bucket, s3_key, str(local_gpkg)) gdf = gpd.read_file(str(local_gpkg)) n = len(gdf) if name == "osm_bridges_subset": n_all = n else: n_lidar = n except Exception: pass return (n_all, n_lidar) from concurrent.futures import ThreadPoolExecutor, as_completed workers = min(32, len(huc_prefixes)) work_items = [(S3_BUCKET, prefix, huc_id, CACHE_DIR, AWS_PROFILE) for huc_id in huc_prefixes] total_all = 0 total_lidar = 0 with ThreadPoolExecutor(max_workers=workers) as executor: futures = {executor.submit(_osm_count_one_huc, item): item for item in work_items} for i, future in enumerate(as_completed(futures)): n_all, n_lidar = future.result() total_all += n_all total_lidar += n_lidar if (i + 1) % 200 == 0: print(f" Processed {i + 1}/{len(huc_prefixes)} HUCs...") all_osm_bridges = total_all osm_bridge_with_lidar = total_lidar counts = {"all_osm_bridges": all_osm_bridges, "osm_bridge_with_lidar": osm_bridge_with_lidar} with open(osm_output_path, "w") as f: json.dump(counts, f, indent=2) print(f"Wrote {osm_output_path}") if all_osm_bridges is not None and osm_bridge_with_lidar is not None: print(f" all_osm_bridges: {all_osm_bridges:,}") print(f" osm_bridge_with_lidar: {osm_bridge_with_lidar:,}")

osm_bridge_counts.json not found on S3 (An error occurred (404) when calling the HeadObject operation: Not Found). Computing from OSM gpkg files (parallel)...
  Processed 200/2139 HUCs...
  Processed 400/2139 HUCs...
  Processed 600/2139 HUCs...
  Processed 800/2139 HUCs...
  Processed 1000/2139 HUCs...
  Processed 1200/2139 HUCs...
  Processed 1400/2139 HUCs...
  Processed 1600/2139 HUCs...
  Processed 1800/2139 HUCs...
  Processed 2000/2139 HUCs...
Wrote /Users/bbhandar/Work/ERT/git/bridge_classification/notebooks/outputs/osm_bridge_counts.json
  all_osm_bridges:       1,073,996
  osm_bridge_with_lidar: 514,125

Stats 3b – OSM bridges: LiDAR vs non-LiDAR (local gpkg)¶

Counts below are from local GeoPackages under data/osm/hucs/{huc_id}/: all subset (osm_bridges_subset__{huc_id}.gpkg), LiDAR subset (osm_bridges_lidar_subset__{huc_id}.gpkg), and non-LiDAR subset (osm_bridges_not_lidar_subset__{huc_id}.gpkg). When data/usgs_entwine/lidar_resources.geojson is present, total lidar source counts (catalog and per LiDAR / non-LiDAR bridges) are also shown; a single bridge can have multiple intersecting lidar sources.

In [3]:

def _count_intersecting_sources(lidar_gdf, bridge_geom_3857, buffer_meters=10.0):
    """Count lidar sources intersecting the buffered bridge geometry (both in 3857)."""
    if lidar_gdf is None or lidar_gdf.empty or bridge_geom_3857 is None:
        return 0
    try:
        buffered = bridge_geom_3857.buffer(buffer_meters)
        possible = list(lidar_gdf.sindex.intersection(buffered.bounds))
        if not possible:
            return 0
        candidates = lidar_gdf.iloc[possible]
        return int(candidates.intersects(buffered).sum())
    except Exception:
        return 0

def _count_one_huc(osm_hucs_dir: str, huc_id: str, lidar_gdf=None, buffer_meters=10.0):
    """Return (n_all, n_lidar, n_not_lidar, total_src_lidar, total_src_not_lidar,
              n_intersecting_lidar, n_intersecting_not_lidar) for one HUC."""
    from pathlib import Path
    import geopandas as gpd
    base = Path(osm_hucs_dir) / huc_id
    n_all, n_lidar, n_not_lidar = 0, 0, 0
    total_src_lidar, total_src_not_lidar = 0, 0
    n_intersecting_lidar, n_intersecting_not_lidar = 0, 0
    for name, suffix in [
        ("all", f"osm_bridges_subset__{huc_id}.gpkg"),
        ("lidar", f"osm_bridges_lidar_subset__{huc_id}.gpkg"),
        ("not_lidar", f"osm_bridges_not_lidar_subset__{huc_id}.gpkg"),
    ]:
        p = base / suffix
        if p.exists():
            gdf = gpd.read_file(str(p))
            n = len(gdf)
            if name == "all":
                n_all = n
            elif name == "lidar":
                n_lidar = n
                if lidar_gdf is not None and not lidar_gdf.empty and n > 0 and gdf.crs is not None:
                    gdf_3857 = gdf.to_crs(epsg=3857)
                    for geom in gdf_3857.geometry:
                        src_count = _count_intersecting_sources(lidar_gdf, geom, buffer_meters)
                        total_src_lidar += src_count
                        if src_count > 0:
                            n_intersecting_lidar += 1
            else:
                n_not_lidar = n
                if lidar_gdf is not None and not lidar_gdf.empty and n > 0 and gdf.crs is not None:
                    gdf_3857 = gdf.to_crs(epsg=3857)
                    for geom in gdf_3857.geometry:
                        src_count = _count_intersecting_sources(lidar_gdf, geom, buffer_meters)
                        total_src_not_lidar += src_count
                        if src_count > 0:
                            n_intersecting_not_lidar += 1
    return (n_all, n_lidar, n_not_lidar, total_src_lidar, total_src_not_lidar,
            n_intersecting_lidar, n_intersecting_not_lidar)

# Resolve local OSM HUCs dir and lidar resources path
_repo_root = Path(NOTEBOOK_DIR).parent
_osm_hucs = _repo_root / "data" / "osm" / "hucs"
if not _osm_hucs.exists():
    _osm_hucs = Path(CACHE_DIR) / "osm" / "hucs"
# _lidar_resources_path = _repo_root / "data" / "usgs_entwine" / "lidar_resources_jan_22_2026.geojson"
_lidar_resources_path = _repo_root / "data" / "usgs_entwine" / "lidar_resources.geojson"
_lidar_gdf = None
if _lidar_resources_path.exists():
    try:
        _lidar_gdf = gpd.read_file(str(_lidar_resources_path)).to_crs(epsg=3857)
    except Exception:
        pass
if not _osm_hucs.exists():
    print("Local OSM HUCs dir not found (tried repo data/osm/hucs and cache osm/hucs). Skip.")
else:
    huc_ids = sorted(d.name for d in _osm_hucs.iterdir() if d.is_dir())
    total_all = total_lidar = total_not_lidar = 0
    total_sources_lidar = total_sources_not_lidar = 0
    total_intersecting_lidar = total_intersecting_not_lidar = 0
    # Sequential loop: pyproj's PROJ SQLite database is not thread-safe,
    # so ThreadPoolExecutor causes UnicodeDecodeError during CRS operations.
    for i, huc_id in enumerate(huc_ids):
        a, b, c, s_li, s_not, i_li, i_not = _count_one_huc(str(_osm_hucs), huc_id, _lidar_gdf)
        total_all += a
        total_lidar += b
        total_not_lidar += c
        total_sources_lidar += s_li
        total_sources_not_lidar += s_not
        total_intersecting_lidar += i_li
        total_intersecting_not_lidar += i_not
        if (i + 1) % 500 == 0:
            print(f"  Processed {i + 1}/{len(huc_ids)} HUCs...")
    print("From local gpkg files (data/osm/hucs/{huc_id}/):")
    print(f"  total bridges (all):  {total_all:,}")
    print(f"  LiDAR bridges:        {total_lidar:,}")
    print(f"  non-LiDAR bridges:   {total_not_lidar:,}")
    if _lidar_gdf is not None:
        print(f"  total lidar sources (catalog):                    {len(_lidar_gdf):,}")
        print(f"  total lidar source pairs (LiDAR bridges):         {total_sources_lidar:,}")
        print(f"  total lidar source pairs (non-LiDAR bridges):     {total_sources_not_lidar:,}")
        print(f"  LiDAR bridges intersecting sources:               {total_intersecting_lidar:,}  (of {total_lidar:,})")
        print(f"  LiDAR bridges NOT intersecting sources:           {total_lidar - total_intersecting_lidar:,}")
        print(f"  non-LiDAR bridges intersecting sources:           {total_intersecting_not_lidar:,}  (of {total_not_lidar:,})")
        print(f"  non-LiDAR bridges NOT intersecting sources:       {total_not_lidar - total_intersecting_not_lidar:,}")
    else:
        print("  (lidar_resources.geojson not found; source counts skipped.)")

def _count_intersecting_sources(lidar_gdf, bridge_geom_3857, buffer_meters=10.0): """Count lidar sources intersecting the buffered bridge geometry (both in 3857).""" if lidar_gdf is None or lidar_gdf.empty or bridge_geom_3857 is None: return 0 try: buffered = bridge_geom_3857.buffer(buffer_meters) possible = list(lidar_gdf.sindex.intersection(buffered.bounds)) if not possible: return 0 candidates = lidar_gdf.iloc[possible] return int(candidates.intersects(buffered).sum()) except Exception: return 0 def _count_one_huc(osm_hucs_dir: str, huc_id: str, lidar_gdf=None, buffer_meters=10.0): """Return (n_all, n_lidar, n_not_lidar, total_src_lidar, total_src_not_lidar, n_intersecting_lidar, n_intersecting_not_lidar) for one HUC.""" from pathlib import Path import geopandas as gpd base = Path(osm_hucs_dir) / huc_id n_all, n_lidar, n_not_lidar = 0, 0, 0 total_src_lidar, total_src_not_lidar = 0, 0 n_intersecting_lidar, n_intersecting_not_lidar = 0, 0 for name, suffix in [ ("all", f"osm_bridges_subset__{huc_id}.gpkg"), ("lidar", f"osm_bridges_lidar_subset__{huc_id}.gpkg"), ("not_lidar", f"osm_bridges_not_lidar_subset__{huc_id}.gpkg"), ]: p = base / suffix if p.exists(): gdf = gpd.read_file(str(p)) n = len(gdf) if name == "all": n_all = n elif name == "lidar": n_lidar = n if lidar_gdf is not None and not lidar_gdf.empty and n > 0 and gdf.crs is not None: gdf_3857 = gdf.to_crs(epsg=3857) for geom in gdf_3857.geometry: src_count = _count_intersecting_sources(lidar_gdf, geom, buffer_meters) total_src_lidar += src_count if src_count > 0: n_intersecting_lidar += 1 else: n_not_lidar = n if lidar_gdf is not None and not lidar_gdf.empty and n > 0 and gdf.crs is not None: gdf_3857 = gdf.to_crs(epsg=3857) for geom in gdf_3857.geometry: src_count = _count_intersecting_sources(lidar_gdf, geom, buffer_meters) total_src_not_lidar += src_count if src_count > 0: n_intersecting_not_lidar += 1 return (n_all, n_lidar, n_not_lidar, total_src_lidar, total_src_not_lidar, n_intersecting_lidar, n_intersecting_not_lidar) # Resolve local OSM HUCs dir and lidar resources path _repo_root = Path(NOTEBOOK_DIR).parent _osm_hucs = _repo_root / "data" / "osm" / "hucs" if not _osm_hucs.exists(): _osm_hucs = Path(CACHE_DIR) / "osm" / "hucs" # _lidar_resources_path = _repo_root / "data" / "usgs_entwine" / "lidar_resources_jan_22_2026.geojson" _lidar_resources_path = _repo_root / "data" / "usgs_entwine" / "lidar_resources.geojson" _lidar_gdf = None if _lidar_resources_path.exists(): try: _lidar_gdf = gpd.read_file(str(_lidar_resources_path)).to_crs(epsg=3857) except Exception: pass if not _osm_hucs.exists(): print("Local OSM HUCs dir not found (tried repo data/osm/hucs and cache osm/hucs). Skip.") else: huc_ids = sorted(d.name for d in _osm_hucs.iterdir() if d.is_dir()) total_all = total_lidar = total_not_lidar = 0 total_sources_lidar = total_sources_not_lidar = 0 total_intersecting_lidar = total_intersecting_not_lidar = 0 # Sequential loop: pyproj's PROJ SQLite database is not thread-safe, # so ThreadPoolExecutor causes UnicodeDecodeError during CRS operations. for i, huc_id in enumerate(huc_ids): a, b, c, s_li, s_not, i_li, i_not = _count_one_huc(str(_osm_hucs), huc_id, _lidar_gdf) total_all += a total_lidar += b total_not_lidar += c total_sources_lidar += s_li total_sources_not_lidar += s_not total_intersecting_lidar += i_li total_intersecting_not_lidar += i_not if (i + 1) % 500 == 0: print(f" Processed {i + 1}/{len(huc_ids)} HUCs...") print("From local gpkg files (data/osm/hucs/{huc_id}/):") print(f" total bridges (all): {total_all:,}") print(f" LiDAR bridges: {total_lidar:,}") print(f" non-LiDAR bridges: {total_not_lidar:,}") if _lidar_gdf is not None: print(f" total lidar sources (catalog): {len(_lidar_gdf):,}") print(f" total lidar source pairs (LiDAR bridges): {total_sources_lidar:,}") print(f" total lidar source pairs (non-LiDAR bridges): {total_sources_not_lidar:,}") print(f" LiDAR bridges intersecting sources: {total_intersecting_lidar:,} (of {total_lidar:,})") print(f" LiDAR bridges NOT intersecting sources: {total_lidar - total_intersecting_lidar:,}") print(f" non-LiDAR bridges intersecting sources: {total_intersecting_not_lidar:,} (of {total_not_lidar:,})") print(f" non-LiDAR bridges NOT intersecting sources: {total_not_lidar - total_intersecting_not_lidar:,}") else: print(" (lidar_resources.geojson not found; source counts skipped.)")

  Processed 500/2139 HUCs...
  Processed 1000/2139 HUCs...
  Processed 1500/2139 HUCs...
  Processed 2000/2139 HUCs...
From local gpkg files (data/osm/hucs/{huc_id}/):
  total bridges (all):  1,073,996
  LiDAR bridges:        514,125
  non-LiDAR bridges:   559,365
  total lidar sources (catalog):                    2,243
  total lidar source pairs (LiDAR bridges):         770,703
  total lidar source pairs (non-LiDAR bridges):     769,484
  LiDAR bridges intersecting sources:               514,125  (of 514,125)
  LiDAR bridges NOT intersecting sources:           0
  non-LiDAR bridges intersecting sources:           469,450  (of 559,365)
  non-LiDAR bridges NOT intersecting sources:       89,915

HUC8 boundaries (S3)¶

Download WBD HUC8 boundaries from S3 (e.g. WBDHU8_webproj.gpkg) to the cache dir and load with geopandas for the map.

In [ ]:

huc8_gdf = None
huc8_local = Path(CACHE_DIR) / "WBDHU8_webproj.gpkg"
try:
    s3.download_file(S3_BUCKET, S3_HUC8_KEY, str(huc8_local))
    huc8_gdf = gpd.read_file(str(huc8_local))
    print(f"Loaded {len(huc8_gdf)} HUC8 features from S3 ({S3_HUC8_KEY}).")
except Exception as e:
    print(f"HUC8 boundaries download failed: {e}")

huc8_gdf = None huc8_local = Path(CACHE_DIR) / "WBDHU8_webproj.gpkg" try: s3.download_file(S3_BUCKET, S3_HUC8_KEY, str(huc8_local)) huc8_gdf = gpd.read_file(str(huc8_local)) print(f"Loaded {len(huc8_gdf)} HUC8 features from S3 ({S3_HUC8_KEY}).") except Exception as e: print(f"HUC8 boundaries download failed: {e}")

Loaded 2399 HUC8 features from S3 (bridge-classification/hucs/WBDHU8_webproj.gpkg).

Class distribution (horizontal bar chart)¶

Uses class_weights.json total_counts; class names and colors from Config.

In [ ]:

weights_path = Path(CACHE_DIR) / "class_weights.json"
with open(weights_path) as f:
    cw = json.load(f)
counts = cw["total_counts"]
values = [counts[str(i)] for i in range(len(CLASS_NAMES))]
total = sum(values)

fig, ax = plt.subplots(figsize=(10, 4))
y_pos = range(len(CLASS_NAMES))
bars = ax.barh(y_pos, values, color=CLASS_COLORS[: len(CLASS_NAMES)], edgecolor="white", height=0.6)
ax.set_yticks(y_pos)
ax.set_yticklabels(CLASS_NAMES, fontsize=13)
ax.invert_yaxis()
ax.set_xlabel("Total Points", fontsize=13)
ax.set_title("Class Distribution Across Training Set", fontsize=16, fontweight="bold")
ax.tick_params(labelsize=11)
ax.grid(True, alpha=0.3, axis="x")
for bar, val in zip(bars, values):
    pct = val / total * 100
    label = f"{val / 1e6:.1f}M ({pct:.1f}%)"
    ax.text(bar.get_width() + total * 0.01, bar.get_y() + bar.get_height() / 2,
            label, va="center", fontsize=12, fontweight="bold")
ax.set_xlim(0, max(values) * 1.25)
ax.xaxis.set_major_formatter(ticker.FuncFormatter(lambda x, _: f"{x/1e6:.0f}M"))
plt.tight_layout()
out_path = Path(OUTPUT_DIR) / "class_distribution.png"
fig.savefig(out_path, dpi=300, bbox_inches="tight", facecolor="white")
plt.show()
print(f"Saved: {out_path}")

weights_path = Path(CACHE_DIR) / "class_weights.json" with open(weights_path) as f: cw = json.load(f) counts = cw["total_counts"] values = [counts[str(i)] for i in range(len(CLASS_NAMES))] total = sum(values) fig, ax = plt.subplots(figsize=(10, 4)) y_pos = range(len(CLASS_NAMES)) bars = ax.barh(y_pos, values, color=CLASS_COLORS[: len(CLASS_NAMES)], edgecolor="white", height=0.6) ax.set_yticks(y_pos) ax.set_yticklabels(CLASS_NAMES, fontsize=13) ax.invert_yaxis() ax.set_xlabel("Total Points", fontsize=13) ax.set_title("Class Distribution Across Training Set", fontsize=16, fontweight="bold") ax.tick_params(labelsize=11) ax.grid(True, alpha=0.3, axis="x") for bar, val in zip(bars, values): pct = val / total * 100 label = f"{val / 1e6:.1f}M ({pct:.1f}%)" ax.text(bar.get_width() + total * 0.01, bar.get_y() + bar.get_height() / 2, label, va="center", fontsize=12, fontweight="bold") ax.set_xlim(0, max(values) * 1.25) ax.xaxis.set_major_formatter(ticker.FuncFormatter(lambda x, _: f"{x/1e6:.0f}M")) plt.tight_layout() out_path = Path(OUTPUT_DIR) / "class_distribution.png" fig.savefig(out_path, dpi=300, bbox_inches="tight", facecolor="white") plt.show() print(f"Saved: {out_path}")

Saved: /Users/bbhandar/Work/ERT/git/bridge_classification/notebooks/outputs/class_distribution.png

Per-bridge point count histogram¶

Requires SILVER_NORMALIZED_DIR set in Config to a local path containing per-bridge JSONs (stats.point_count). Skip if not set.

In [ ]:

_repo_root = Path.cwd()
if (_repo_root / "utils").is_dir():
    pass
elif (_repo_root.parent / "utils").is_dir():
    _repo_root = _repo_root.parent

def get_point_count_from_json(path: str | Path) -> Optional[int]:
    try:
        with open(path) as f:
            data = json.load(f)
        if "stats" in data and "point_count" in data["stats"]:
            return data["stats"]["point_count"]
    except (OSError, json.JSONDecodeError, KeyError, TypeError):
        pass
    return None


if SILVER_NORMALIZED_DIR is None or not Path(_repo_root / SILVER_NORMALIZED_DIR).exists():
    print("SILVER_NORMALIZED_DIR not set or path does not exist. Set it in Config to a local silver_training_normalized directory to generate the per-bridge point count histogram.")
else:
    silver_path = Path(_repo_root / SILVER_NORMALIZED_DIR)
    json_files = list(silver_path.rglob("*.json"))
    max_workers = min(32, (os.cpu_count() or 4) * 2)
    print(f"Reading point_count from {len(json_files):,} JSONs (max_workers={max_workers})...")
    with ThreadPoolExecutor(max_workers=max_workers) as executor:
        results = list(executor.map(get_point_count_from_json, json_files, chunksize=2000))
    point_counts = [c for c in results if c is not None]
    if not point_counts:
        print("No point_count found in JSONs under SILVER_NORMALIZED_DIR.")
    else:
        point_counts = np.array(point_counts)
        median_pc = np.median(point_counts)
        mean_pc = np.mean(point_counts)
        print(f"Bridges: {len(point_counts):,}, Median: {median_pc:,.0f}, Mean: {mean_pc:,.0f}")
        print(f"Range: [{point_counts.min():,} — {point_counts.max():,}]")
        fig, ax = plt.subplots(figsize=(10, 5))
        log_bins = np.logspace(np.log10(max(1, point_counts.min())), np.log10(point_counts.max()), 40)
        ax.hist(point_counts, bins=log_bins, color="#2196F3", edgecolor="white", linewidth=0.5, alpha=0.85)
        ax.set_xscale("log")
        ax.axvline(median_pc, color="#E91E63", linestyle="--", linewidth=2, label=f"Median: {median_pc:,.0f}")
        ax.axvline(mean_pc, color="#FF9800", linestyle="--", linewidth=2, label=f"Mean: {mean_pc:,.0f}")
        ax.set_title(f"Per-bridge point count (training set, n={len(point_counts):,})", fontsize=16, fontweight="bold")
        ax.set_xlabel("Points per Bridge (log scale)", fontsize=13)
        ax.set_ylabel("Count", fontsize=13)
        ax.tick_params(labelsize=11)
        ax.legend(fontsize=12)
        ax.grid(True, alpha=0.3, axis="y")
        plt.tight_layout()
        out_path = Path(OUTPUT_DIR) / "bridge_point_distribution.png"
        fig.savefig(out_path, dpi=300, bbox_inches="tight", facecolor="white")
        plt.show()
        print(f"Saved: {out_path}")

_repo_root = Path.cwd() if (_repo_root / "utils").is_dir(): pass elif (_repo_root.parent / "utils").is_dir(): _repo_root = _repo_root.parent def get_point_count_from_json(path: str | Path) -> Optional[int]: try: with open(path) as f: data = json.load(f) if "stats" in data and "point_count" in data["stats"]: return data["stats"]["point_count"] except (OSError, json.JSONDecodeError, KeyError, TypeError): pass return None if SILVER_NORMALIZED_DIR is None or not Path(_repo_root / SILVER_NORMALIZED_DIR).exists(): print("SILVER_NORMALIZED_DIR not set or path does not exist. Set it in Config to a local silver_training_normalized directory to generate the per-bridge point count histogram.") else: silver_path = Path(_repo_root / SILVER_NORMALIZED_DIR) json_files = list(silver_path.rglob("*.json")) max_workers = min(32, (os.cpu_count() or 4) * 2) print(f"Reading point_count from {len(json_files):,} JSONs (max_workers={max_workers})...") with ThreadPoolExecutor(max_workers=max_workers) as executor: results = list(executor.map(get_point_count_from_json, json_files, chunksize=2000)) point_counts = [c for c in results if c is not None] if not point_counts: print("No point_count found in JSONs under SILVER_NORMALIZED_DIR.") else: point_counts = np.array(point_counts) median_pc = np.median(point_counts) mean_pc = np.mean(point_counts) print(f"Bridges: {len(point_counts):,}, Median: {median_pc:,.0f}, Mean: {mean_pc:,.0f}") print(f"Range: [{point_counts.min():,} — {point_counts.max():,}]") fig, ax = plt.subplots(figsize=(10, 5)) log_bins = np.logspace(np.log10(max(1, point_counts.min())), np.log10(point_counts.max()), 40) ax.hist(point_counts, bins=log_bins, color="#2196F3", edgecolor="white", linewidth=0.5, alpha=0.85) ax.set_xscale("log") ax.axvline(median_pc, color="#E91E63", linestyle="--", linewidth=2, label=f"Median: {median_pc:,.0f}") ax.axvline(mean_pc, color="#FF9800", linestyle="--", linewidth=2, label=f"Mean: {mean_pc:,.0f}") ax.set_title(f"Per-bridge point count (training set, n={len(point_counts):,})", fontsize=16, fontweight="bold") ax.set_xlabel("Points per Bridge (log scale)", fontsize=13) ax.set_ylabel("Count", fontsize=13) ax.tick_params(labelsize=11) ax.legend(fontsize=12) ax.grid(True, alpha=0.3, axis="y") plt.tight_layout() out_path = Path(OUTPUT_DIR) / "bridge_point_distribution.png" fig.savefig(out_path, dpi=300, bbox_inches="tight", facecolor="white") plt.show() print(f"Saved: {out_path}")

Reading point_count from 597,426 JSONs (max_workers=28)...
Bridges: 597,426, Median: 4,113, Mean: 8,081
Range: [20 — 5,667,854]

Saved: /Users/bbhandar/Work/ERT/git/bridge_classification/notebooks/outputs/bridge_point_distribution.png

Map: HUC8 distribution¶

Join split HUC8 IDs to the boundaries (match by HUC code). Plot; HUC8s in the dataset are colored, others gray. Save to notebooks/outputs/huc8_distribution.png.

In [ ]:

if huc8_gdf is None or len(huc8_gdf) == 0:
    print("No HUC8 boundaries loaded; skipping map.")
else:
    # Find HUC code column (WBD often uses huc8, HUC8, or similar)
    cols = [c for c in huc8_gdf.columns if c and c.upper().replace(" ", "") == "HUC8"]
    if not cols:
        cols = [c for c in huc8_gdf.columns if "huc" in c.lower() and "8" in c]
    if not cols:
        cols = [c for c in huc8_gdf.columns if "huc" in c.lower()]
    huc_col = cols[0] if cols else huc8_gdf.columns[0]

    # Normalize to 8-digit string for join (leading zeros)
    def norm(s):
        s = str(s).strip()
        if s.isdigit():
            return s.zfill(8)
        return s

    huc8_gdf["_huc8_norm"] = huc8_gdf[huc_col].apply(norm)
    in_split = huc8_gdf["_huc8_norm"].isin(hucs_all)

    # Compute tight bounds using per-feature quantiles to exclude remote outlier
    # territories (e.g. Pacific islands far from the main US landmass) that would
    # otherwise create large blank whitespace to the right/bottom of the figure.
    geom_bounds = huc8_gdf.geometry.bounds
    minx = geom_bounds["minx"].quantile(0.01)
    miny = geom_bounds["miny"].quantile(0.01)
    maxx = geom_bounds["maxx"].quantile(0.99)
    maxy = geom_bounds["maxy"].quantile(0.99)

    data_aspect = (maxx - minx) / (maxy - miny)
    fig_height = 22 / data_aspect
    fig, ax = plt.subplots(1, 1, figsize=(22, fig_height))
    huc8_gdf[~in_split].plot(ax=ax, color="#e8e8e8", edgecolor="#cccccc", linewidth=0.15, label="_nolegend_")
    huc8_gdf[in_split].plot(ax=ax, color="#1e7bb8", edgecolor="#0d3d5c", linewidth=0.25, label="_nolegend_")
    ax.set_xlim(minx, maxx)
    ax.set_ylim(miny, maxy)
    ax.set_title("HUC8 distribution in dataset split", fontsize=20, fontweight="bold")
    ax.legend(
        handles=[
            Patch(facecolor="#1e7bb8", edgecolor="#0d3d5c", label="In dataset"),
            Patch(facecolor="#e8e8e8", edgecolor="#cccccc", label="Not in dataset"),
        ],
        loc="lower left",
        fontsize=12,
        framealpha=0.95,
    )
    ax.set_axis_off()
    out_path = Path(OUTPUT_DIR) / "huc8_distribution.png"
    plt.savefig(out_path, dpi=200, bbox_inches="tight", pad_inches=0.1)
    plt.show()
    print(f"Saved: {out_path}")

if huc8_gdf is None or len(huc8_gdf) == 0: print("No HUC8 boundaries loaded; skipping map.") else: # Find HUC code column (WBD often uses huc8, HUC8, or similar) cols = [c for c in huc8_gdf.columns if c and c.upper().replace(" ", "") == "HUC8"] if not cols: cols = [c for c in huc8_gdf.columns if "huc" in c.lower() and "8" in c] if not cols: cols = [c for c in huc8_gdf.columns if "huc" in c.lower()] huc_col = cols[0] if cols else huc8_gdf.columns[0] # Normalize to 8-digit string for join (leading zeros) def norm(s): s = str(s).strip() if s.isdigit(): return s.zfill(8) return s huc8_gdf["_huc8_norm"] = huc8_gdf[huc_col].apply(norm) in_split = huc8_gdf["_huc8_norm"].isin(hucs_all) # Compute tight bounds using per-feature quantiles to exclude remote outlier # territories (e.g. Pacific islands far from the main US landmass) that would # otherwise create large blank whitespace to the right/bottom of the figure. geom_bounds = huc8_gdf.geometry.bounds minx = geom_bounds["minx"].quantile(0.01) miny = geom_bounds["miny"].quantile(0.01) maxx = geom_bounds["maxx"].quantile(0.99) maxy = geom_bounds["maxy"].quantile(0.99) data_aspect = (maxx - minx) / (maxy - miny) fig_height = 22 / data_aspect fig, ax = plt.subplots(1, 1, figsize=(22, fig_height)) huc8_gdf[~in_split].plot(ax=ax, color="#e8e8e8", edgecolor="#cccccc", linewidth=0.15, label="_nolegend_") huc8_gdf[in_split].plot(ax=ax, color="#1e7bb8", edgecolor="#0d3d5c", linewidth=0.25, label="_nolegend_") ax.set_xlim(minx, maxx) ax.set_ylim(miny, maxy) ax.set_title("HUC8 distribution in dataset split", fontsize=20, fontweight="bold") ax.legend( handles=[ Patch(facecolor="#1e7bb8", edgecolor="#0d3d5c", label="In dataset"), Patch(facecolor="#e8e8e8", edgecolor="#cccccc", label="Not in dataset"), ], loc="lower left", fontsize=12, framealpha=0.95, ) ax.set_axis_off() out_path = Path(OUTPUT_DIR) / "huc8_distribution.png" plt.savefig(out_path, dpi=200, bbox_inches="tight", pad_inches=0.1) plt.show() print(f"Saved: {out_path}")

Saved: /Users/bbhandar/Work/ERT/git/bridge_classification/notebooks/outputs/huc8_distribution.png