From 2f52e953c2924742a81d50551a718eea98588743 Mon Sep 17 00:00:00 2001 From: WagnerJ-Dev Date: Sat, 27 Jun 2026 13:16:56 -0400 Subject: [PATCH 01/13] =?UTF-8?q?Working=20LiDAR=E2=86=92solar=20pipeline:?= =?UTF-8?q?=20CRS=20fix,=20in-memory=20clip,=20panel-packing=20energy=20mo?= =?UTF-8?q?del?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- memory/tech_stack.md | 8 +- requirements.txt | 2 +- test_pipeline.py | 291 +++++++++++++++++++++++++++++-------------- 3 files changed, 206 insertions(+), 95 deletions(-) diff --git a/memory/tech_stack.md b/memory/tech_stack.md index 504d80d..76f4c19 100644 --- a/memory/tech_stack.md +++ b/memory/tech_stack.md @@ -17,8 +17,14 @@ a bottleneck. **Frontend:** FastAPI + Leaflet.js + OpenStreetMap tiles. **Environment:** -- Mac: `brew install pdal` then `pip install -r requirements.txt` (no conda needed) +- Mac: `brew install pdal` then `pip install -r requirements.txt` in a venv (no conda needed) - Windows: conda-forge environment via environment.yml (PDAL native deps need conda on Windows) +**PDAL pip gotcha (hit 2026-06-25):** The PyPI package is `pdal`, NOT `python-pdal`. +`python-pdal` is the *conda* name; on PyPI it's an empty 0.0.1 stub that installs +nothing. requirements.txt has been corrected to `pdal`. The pip `pdal` wheel builds +against the Homebrew libpdal via pdal-config, so `brew install pdal` must come first. +Verified working: Homebrew PDAL 2.10.2 + pip pdal 3.5.3 + Python 3.9.6 venv. + **How to apply:** User prefers simple over complex. Don't suggest adding servers, Docker, or extra dependencies unless absolutely necessary. diff --git a/requirements.txt b/requirements.txt index ac96e1d..a7b6f69 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,4 @@ -python-pdal +pdal pvlib open3d duckdb diff --git a/test_pipeline.py b/test_pipeline.py index 1960139..2190008 100644 --- a/test_pipeline.py +++ b/test_pipeline.py @@ -19,22 +19,45 @@ import pdal import pvlib import requests +import shapely +from pyproj import Transformer from scipy.spatial import ConvexHull from shapely.geometry import Polygon +from shapely.ops import transform as shapely_transform warnings.filterwarnings("ignore") # ── Configuration ───────────────────────────────────────────────────────────── -# Small test area: suburban Chester County, PA (near Exton) -# Adjust this bbox to pick a different area — keep it small (0.02° × 0.02° max) -TEST_BBOX = (-75.640, 40.015, -75.620, 40.030) # (west, south, east, north) -TEST_LAT = 40.022 -TEST_LON = -75.630 +# Small test area: residential single-family streets, West Chester borough, PA. +# Picked because the goal is scoring HOMES — the old Exton bbox (below) was all +# commercial/retail. Adjust these 4 numbers to target a different area; keep it +# small (0.02° × 0.02° max). (west, south, east, north) +TEST_BBOX = (-75.615, 39.955, -75.600, 39.970) +TEST_LAT = 39.962 +TEST_LON = -75.607 + +# Previous test area (Exton commercial strip — big-box roofs, no houses): +# TEST_BBOX = (-75.640, 40.015, -75.620, 40.030); TEST_LAT/LON = 40.022, -75.630 TILE_CACHE_DIR = "tile_cache" # LiDAR tiles cached here (re-used on re-runs) DB_PATH = "solar_results.duckdb" +# ── Panel / system model ────────────────────────────────────────────────────── +# All physical: change these to model different panels/configurations. The energy +# model packs WHOLE panels onto each roof plane after a fire-code edge setback, +# rather than assuming a fixed % of the roof is usable. +PANEL_WIDTH_M = 1.05 # standard residential module footprint (~1.05 × 1.74 m) +PANEL_HEIGHT_M = 1.74 +PANEL_WATTS = 400 # module rated DC power at STC (1000 W/m²) +ROOF_SETBACK_M = 0.46 # fire-code clear pathway from plane edges (~18 inches) +PANEL_PACKING = 0.90 # fraction of the setback area that packs into whole modules +SYSTEM_LOSSES = 0.86 # inverter + wiring + soiling (~14% total loss) +MAX_AZIMUTH_OFFSET = 120 # skip planes facing >this many degrees from due south + +PANEL_AREA_M2 = PANEL_WIDTH_M * PANEL_HEIGHT_M # m² per module +MODULE_EFF = PANEL_WATTS / (PANEL_AREA_M2 * 1000.0) # derived STC efficiency + MIN_ROOF_POINTS = 20 # skip buildings with fewer LiDAR points than this RANSAC_ITERATIONS = 300 # lower = faster, higher = more accurate RANSAC_DISTANCE_THRESHOLD = 0.20 # meters; how far a point can be from a plane @@ -129,6 +152,7 @@ def get_buildings_from_osm(bbox): r = requests.post( "https://overpass-api.de/api/interpreter", data={"data": query}, + headers={"User-Agent": "SolarGrader/1.0 (joshuawagnersd@gmail.com)"}, timeout=90, ) r.raise_for_status() @@ -172,7 +196,7 @@ def get_tmy_data(lat, lon): Returns a DataFrame with hourly irradiance columns: ghi, dni, dhi. """ print(f" Fetching TMY data from PVGIS for ({lat:.3f}, {lon:.3f})...") - tmy, _, _, _ = pvlib.iotools.get_pvgis_tmy( + tmy, _ = pvlib.iotools.get_pvgis_tmy( # pvlib >=0.11 returns (data, metadata) latitude=lat, longitude=lon, outputformat="json", @@ -185,49 +209,68 @@ def get_tmy_data(lat, lon): # ── Step 4: LiDAR clipping ──────────────────────────────────────────────────── -def clip_lidar_to_building(tile_path, footprint_poly, buffer_deg=0.00002): +def get_tile_srs(tile_path): + """Read the LiDAR tile's spatial reference (WKT) from its LAS header.""" + p = pdal.Pipeline(json.dumps( + {"pipeline": [{"type": "readers.las", "filename": tile_path.replace("\\", "/"), "count": 1}]} + )) + p.execute() + md = p.metadata["metadata"]["readers.las"] + return md.get("comp_spatialreference") or md.get("spatialreference") + + +def load_tile_points(tile_path): + """ + Read the entire LiDAR tile into memory ONCE as a structured numpy array. + Returns (points, X, Y) where X/Y are float views for fast bbox masking. + Reading once and filtering per-building in RAM is ~100x faster than + re-opening the 145 MB file for every building. + """ + p = pdal.Pipeline(json.dumps( + {"pipeline": [{"type": "readers.las", "filename": tile_path.replace("\\", "/")}]} + )) + p.execute() + pts = p.arrays[0] + return pts, np.asarray(pts["X"], dtype=float), np.asarray(pts["Y"], dtype=float) + + +def clip_lidar_to_building(tile_points, footprint_poly, transformer, buffer_m=0.5): """ - Use PDAL to extract LiDAR points that fall within a building footprint. + Extract LiDAR points inside a building footprint from already-loaded tile + points. The footprint comes in lon/lat (EPSG:4326); the tile is in projected + meters, so we reproject the polygon into the tile's CRS before clipping. Returns an (N, 3) numpy array of [X, Y, Z], or None if too few points. """ - buffered_wkt = footprint_poly.buffer(buffer_deg).wkt - - pipeline_json = json.dumps( - { - "pipeline": [ - {"type": "readers.las", "filename": tile_path.replace("\\", "/")}, - {"type": "filters.crop", "polygon": buffered_wkt}, - ] - } - ) + pts, X, Y = tile_points + proj_poly = shapely_transform(transformer.transform, footprint_poly).buffer(buffer_m) - try: - pipeline = pdal.Pipeline(pipeline_json) - pipeline.execute() - arrays = pipeline.arrays - if not arrays or len(arrays[0]) == 0: - return None - - pts = arrays[0] - - # Prefer building-classified points (class 6); fall back to height filter - if "Classification" in pts.dtype.names: - roof_pts = pts[pts["Classification"] == 6] - if len(roof_pts) < MIN_ROOF_POINTS: - z_min = pts["Z"].min() - roof_pts = pts[pts["Z"] > z_min + 2.0] - else: - z_min = pts["Z"].min() - roof_pts = pts[pts["Z"] > z_min + 2.0] + # Fast bounding-box prefilter, then exact point-in-polygon on the small subset + minx, miny, maxx, maxy = proj_poly.bounds + bbox_mask = (X >= minx) & (X <= maxx) & (Y >= miny) & (Y <= maxy) + if not bbox_mask.any(): + return None - if len(roof_pts) < MIN_ROOF_POINTS: - return None + sub = pts[bbox_mask] + inside = shapely.contains_xy(proj_poly, sub["X"], sub["Y"]) + sub = sub[inside] + if len(sub) == 0: + return None - return np.column_stack([roof_pts["X"], roof_pts["Y"], roof_pts["Z"]]) + # Prefer building-classified points (class 6); fall back to height filter + if "Classification" in sub.dtype.names: + roof_pts = sub[sub["Classification"] == 6] + if len(roof_pts) < MIN_ROOF_POINTS: + z_min = sub["Z"].min() + roof_pts = sub[sub["Z"] > z_min + 2.0] + else: + z_min = sub["Z"].min() + roof_pts = sub[sub["Z"] > z_min + 2.0] - except Exception: + if len(roof_pts) < MIN_ROOF_POINTS: return None + return np.column_stack([roof_pts["X"], roof_pts["Y"], roof_pts["Z"]]) + # ── Step 5: Roof plane extraction (RANSAC) ──────────────────────────────────── @@ -239,6 +282,17 @@ def extract_roof_planes(points_xyz): pcd = o3d.geometry.PointCloud() pcd.points = o3d.utility.Vector3dVector(points_xyz) + # Ground point density (points per m² of horizontal area) from ONE convex + # hull over all points. Used to size each plane by its point count — this + # avoids the inflation you get from summing overlapping per-plane hulls. + # Coordinates are in projected meters (tile CRS), so areas are already m². + all_xy = np.asarray(points_xyz)[:, :2] + try: + covered_m2 = float(ConvexHull(all_xy).volume) if len(all_xy) >= 3 else 0.0 + except Exception: + covered_m2 = 0.0 + ground_density = (len(points_xyz) / covered_m2) if covered_m2 > 0 else 0.0 + planes = [] remaining = pcd @@ -279,19 +333,15 @@ def extract_roof_planes(points_xyz): # In projected coords, +X=east, +Y=north → atan2(nx, ny) gives bearing azimuth_deg = float((np.degrees(np.arctan2(nx, ny)) + 360) % 360) - # Estimate plane area via convex hull of projected inlier points - inlier_pts = np.asarray(remaining.points)[inliers] - try: - if len(inlier_pts) >= 3: - hull = ConvexHull(inlier_pts[:, :2]) - area_m2 = float(hull.volume) # in 2D, .volume = area in coordinate units² - # If coords are in degrees (geographic), convert to m² - if abs(inlier_pts[0, 0]) < 180: # likely lat/lon, not projected - area_m2 *= (111320 ** 2) * np.cos(np.radians(np.mean(inlier_pts[:, 1]))) - else: - area_m2 = 0.0 - except Exception: - area_m2 = float(len(inlier_pts)) # rough fallback: 1 pt ≈ 1 m² + # Plane area = (inlier point count ÷ ground density) gives the plane's + # horizontal (ground-projected) area; dividing by cos(tilt) converts it + # to true tilted surface area. Each point belongs to exactly one plane, + # so summing these areas can never exceed the building's real footprint. + if ground_density > 0: + ground_area_m2 = len(inliers) / ground_density + area_m2 = ground_area_m2 / max(np.cos(np.radians(tilt_deg)), 0.1) + else: + area_m2 = 0.0 planes.append( { @@ -309,31 +359,45 @@ def extract_roof_planes(points_xyz): # ── Step 6: Solar potential calculation (pvlib) ─────────────────────────────── +def panels_on_plane(plane_area_m2): + """ + How many whole panels physically fit on a roof plane of this area. + Approximate the plane as an equal-area square, erode it inward by the + fire-code setback on all sides, apply a packing factor, then divide by the + module footprint. Small/cut-up roofs lose proportionally more area — the + realistic behaviour a flat usable-% can't capture. + """ + side = np.sqrt(plane_area_m2) + usable_side = max(0.0, side - 2 * ROOF_SETBACK_M) + usable_area = (usable_side ** 2) * PANEL_PACKING + return int(usable_area // PANEL_AREA_M2) + + def calculate_annual_kwh(lat, lon, tmy, planes): """ - Calculate total annual kWh potential across all viable roof planes. - Returns (total_kwh, best_plane). + Estimate realistic annual production by packing whole panels onto each viable + roof plane and running pvlib's irradiance model for that plane's orientation. + Returns (total_kwh, best_plane, system_kw_dc, total_panels). """ location = pvlib.location.Location(latitude=lat, longitude=lon, tz="Etc/GMT+5") solar_pos = location.get_solarposition(tmy.index) - PANEL_EFFICIENCY = 0.20 # typical monocrystalline panel - SYSTEM_LOSSES = 0.86 # inverter + wiring + soiling (~14% total loss) - MAX_PANEL_AREA = 150.0 # sanity cap in m² — no residential roof has more - total_kwh = 0.0 + total_panels = 0 best_plane = None best_kwh = 0.0 for plane in planes: - # Skip near-north-facing planes (>120° from south = not worth paneling) - if abs(plane["azimuth_deg"] - 180) > 120: + # Skip planes facing too far from due south to be worth paneling + if abs(plane["azimuth_deg"] - 180) > MAX_AZIMUTH_OFFSET: continue - area = min(plane["area_m2"], MAX_PANEL_AREA) - if area <= 0: + n_panels = panels_on_plane(plane["area_m2"]) + if n_panels <= 0: continue + panel_area = n_panels * PANEL_AREA_M2 + poa = pvlib.irradiance.get_total_irradiance( surface_tilt=plane["tilt_deg"], surface_azimuth=plane["azimuth_deg"], @@ -344,17 +408,18 @@ def calculate_annual_kwh(lat, lon, tmy, planes): dhi=tmy["dhi"], ) - plane_kwh = float( - (poa["poa_global"].fillna(0) * area * PANEL_EFFICIENCY * SYSTEM_LOSSES).sum() - / 1000 - ) + # Annual plane-of-array irradiance (kWh/m²/yr) × panel area × efficiency + poa_annual = float(poa["poa_global"].fillna(0).sum() / 1000) + plane_kwh = poa_annual * panel_area * MODULE_EFF * SYSTEM_LOSSES total_kwh += plane_kwh + total_panels += n_panels if plane_kwh > best_kwh: best_kwh = plane_kwh best_plane = plane - return total_kwh, best_plane + system_kw_dc = total_panels * PANEL_WATTS / 1000.0 + return total_kwh, best_plane, system_kw_dc, total_panels # ── Step 7: Grading ─────────────────────────────────────────────────────────── @@ -391,6 +456,9 @@ def grade_home(annual_kwh, total_area_m2, primary_azimuth_deg): def setup_db(path): con = duckdb.connect(path) + # Fresh table each run so test results reflect only the current tile/bbox + # (the real multi-tile pipeline would accumulate instead of dropping). + con.execute("DROP TABLE IF EXISTS homes") con.execute(""" CREATE TABLE IF NOT EXISTS homes ( osm_id BIGINT PRIMARY KEY, @@ -402,6 +470,8 @@ def setup_db(path): primary_azimuth_deg DOUBLE, roof_plane_count INTEGER, annual_kwh_potential DOUBLE, + system_kw_dc DOUBLE, + panel_count INTEGER, solar_score INTEGER, solar_grade VARCHAR(2), processed_at TIMESTAMP @@ -413,7 +483,7 @@ def setup_db(path): def save_result(con, building, result): con.execute( """ - INSERT OR REPLACE INTO homes VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) + INSERT OR REPLACE INTO homes VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) """, [ building["osm_id"], @@ -425,6 +495,8 @@ def save_result(con, building, result): result["primary_azimuth"], result["plane_count"], result["annual_kwh"], + result["system_kw_dc"], + result["panel_count"], result["score"], result["grade"], pd.Timestamp.now(), @@ -448,6 +520,16 @@ def main(): print("\n[2/5] Downloading tile...") tile_path = download_tile(tile_info["downloadURL"], TILE_CACHE_DIR) + # Build a transformer from lon/lat (OSM) into the tile's projected CRS (meters) + tile_srs = get_tile_srs(tile_path) + transformer = Transformer.from_crs("EPSG:4326", tile_srs, always_xy=True) + print(f" Tile CRS: {tile_srs.split(',')[0].split('[')[-1].strip(chr(34))}") + + # Load the whole tile into memory once (filtered per-building in RAM below) + print(" Loading tile points into memory...") + tile_points = load_tile_points(tile_path) + print(f" Loaded {len(tile_points[0]):,} points.") + # 3. Get buildings print("\n[3/5] Fetching building footprints from OpenStreetMap...") buildings = get_buildings_from_osm(TEST_BBOX) @@ -466,6 +548,7 @@ def main(): counts = {"scored": 0, "skipped_no_points": 0, "skipped_no_planes": 0} grade_counts = {} timings = [] + density_samples = [] # (n_points, footprint_m2) for buildings that had LiDAR t_total = time.time() for i, building in enumerate(buildings): @@ -477,7 +560,9 @@ def main(): ) # Clip LiDAR - points = clip_lidar_to_building(tile_path, building["geometry"]) + points = clip_lidar_to_building(tile_points, building["geometry"], transformer) + if points is not None: + density_samples.append((len(points), building["footprint_m2"])) if points is None: counts["skipped_no_points"] += 1 continue @@ -489,10 +574,10 @@ def main(): continue # Solar calc - annual_kwh, best_plane = calculate_annual_kwh( + annual_kwh, best_plane, system_kw, n_panels = calculate_annual_kwh( building["lat"], building["lon"], tmy, planes ) - if annual_kwh == 0 or best_plane is None: + if annual_kwh == 0 or best_plane is None or n_panels == 0: counts["skipped_no_planes"] += 1 continue @@ -509,6 +594,8 @@ def main(): "primary_azimuth": best_plane["azimuth_deg"], "plane_count": len(planes), "annual_kwh": annual_kwh, + "system_kw_dc": system_kw, + "panel_count": n_panels, "grade": grade, "score": score, }, @@ -541,11 +628,14 @@ def main(): print(f"\n Top 10 A/A+ leads:") df = con.execute(""" SELECT lat, lon, - ROUND(annual_kwh_potential) AS kwh_yr, - ROUND(usable_roof_area_m2) AS roof_m2, - ROUND(primary_tilt_deg) AS tilt, - ROUND(primary_azimuth_deg) AS azimuth, - solar_grade AS grade + ROUND(annual_kwh_potential) AS kwh_yr, + ROUND(system_kw_dc, 1) AS kw_dc, + panel_count AS panels, + ROUND(footprint_m2) AS osm_m2, + ROUND(usable_roof_area_m2) AS roof_m2, + ROUND(primary_tilt_deg) AS tilt, + ROUND(primary_azimuth_deg) AS azimuth, + solar_grade AS grade FROM homes WHERE solar_grade IN ('A+', 'A') ORDER BY annual_kwh_potential DESC @@ -558,23 +648,38 @@ def main(): else: print(df.to_string(index=False)) - # Diagnostic: LiDAR point density check + # Realism check: specific yield (kWh per kW). PA should land ~1,100–1,300. + print(f"\n Realism check (kWh per kW installed — PA expected ~1,100–1,300):") + yields = con.execute(""" + SELECT annual_kwh_potential / system_kw_dc AS specific_yield, system_kw_dc + FROM homes WHERE system_kw_dc > 0 + """).fetchdf() + if not yields.empty: + med_yield = float(yields["specific_yield"].median()) + med_kw = float(yields["system_kw_dc"].median()) + print(f" Median specific yield: ~{med_yield:,.0f} kWh/kW (median system: {med_kw:.1f} kW DC)") + if 950 <= med_yield <= 1450: + print(" OK — production is in the physically realistic range.") + else: + print(" WARNING — outside the expected band; the panel/energy model needs review.") + + # Diagnostic: LiDAR point density check (averaged over buildings that had LiDAR) print(f"\n LiDAR density check:") - sample = buildings[0] if buildings else None - if sample: - pts = clip_lidar_to_building(tile_path, sample["geometry"]) - if pts is not None: - footprint_m2 = sample["footprint_m2"] - density = len(pts) / max(footprint_m2, 1) - print(f" Sample building: {len(pts)} points over ~{footprint_m2:.0f} m²") - print(f" Point density: ~{density:.2f} pts/m²") - if density < 0.5: - print(" WARNING: Low density (<0.5 pts/m²) — RANSAC results may be unreliable.") - print(" Consider finding a higher-density tile for this area.") - elif density >= 4: - print(" Density is good — RANSAC results should be reliable.") - else: - print(" Density is acceptable.") + if density_samples: + densities = [n / max(fp, 1) for n, fp in density_samples] + density = float(np.median(densities)) + print(f" Buildings with LiDAR: {len(density_samples)} of {len(buildings)}") + print(f" Median point density: ~{density:.2f} pts/m²") + if density < 0.5: + print(" WARNING: Low density (<0.5 pts/m²) — RANSAC results may be unreliable.") + print(" Consider finding a higher-density tile for this area.") + elif density >= 4: + print(" Density is good — RANSAC results should be reliable.") + else: + print(" Density is acceptable.") + else: + print(" No buildings overlapped the tile — density unknown.") + print(" The test bbox is likely larger than this single tile's footprint.") print(f"\n Results saved to: {DB_PATH}") print(f" Open with: python -c \"import duckdb; print(duckdb.connect('{DB_PATH}').execute('SELECT * FROM homes LIMIT 5').fetchdf())\"") From 9c1b5b2fbda55c5577c54721014e2562e124849d Mon Sep 17 00:00:00 2001 From: WagnerJ-Dev Date: Sat, 27 Jun 2026 13:17:22 -0400 Subject: [PATCH 02/13] Stubbing --- .claude/hooks/block-secret-commit.sh | 33 +++++++++++ .env.example | 10 ++++ .gitignore | 6 ++ CLAUDE.md | 84 ++++++++++++++++++++++++++++ 4 files changed, 133 insertions(+) create mode 100755 .claude/hooks/block-secret-commit.sh create mode 100644 .env.example create mode 100644 CLAUDE.md diff --git a/.claude/hooks/block-secret-commit.sh b/.claude/hooks/block-secret-commit.sh new file mode 100755 index 0000000..b64426e --- /dev/null +++ b/.claude/hooks/block-secret-commit.sh @@ -0,0 +1,33 @@ +#!/usr/bin/env bash +# PreToolUse(Bash) guard: block `git commit` when about-to-be-committed changes +# contain hardcoded secrets. Enforces CLAUDE.md principle #2 (env vars only). +# +# Scans `git diff HEAD` (staged + unstaged tracked changes — covers both a normal +# staged commit and `git commit -a`). Only inspects ADDED lines. Exit 2 blocks the +# tool call and surfaces stderr back to Claude. Exit 0 allows the commit. +set -uo pipefail + +# Diff of everything tracked vs the last commit; exclude the secrets template, +# which intentionally lists key NAMES with empty values. +diff=$(git diff HEAD --no-color -- . ':(exclude).env.example' 2>/dev/null) || exit 0 + +# Added lines only (drop the +++ file headers). +added=$(printf '%s\n' "$diff" | grep '^+' | grep -v '^+++') || true +[ -z "$added" ] && exit 0 + +# secret-name = "literal-8+chars" | AWS key id | PEM private-key header +patterns='(api[_-]?key|secret|passwd|password|token|access[_-]?key|client[_-]?secret|auth_token)["'"'"' ]*[:=][[:space:]]*["'"'"'][^"'"'"' ]{8,}|AKIA[0-9A-Z]{16}|-----BEGIN [A-Z ]*PRIVATE KEY-----' + +hits=$(printf '%s\n' "$added" | grep -niE "$patterns") || true + +if [ -n "$hits" ]; then + { + echo "COMMIT BLOCKED — staged changes look like they contain a hardcoded secret." + echo "Per CLAUDE.md principle #2, secrets must come from env vars / .env, never inline." + echo "Offending added lines:" + printf '%s\n' "$hits" + echo "If this is a false positive, the user can commit manually or adjust the hook." + } >&2 + exit 2 +fi +exit 0 diff --git a/.env.example b/.env.example new file mode 100644 index 0000000..9864bc3 --- /dev/null +++ b/.env.example @@ -0,0 +1,10 @@ +# Copy this file to .env and fill in real values. .env is gitignored — never commit it. +# Load in Python with: from dotenv import load_dotenv; load_dotenv() +# Read with: os.environ["NREL_API_KEY"] (raises if missing — fail loud, never default a secret) + +# NREL NSRDB / Developer Network — free key from https://developer.nrel.gov/signup/ +# Needed for solar irradiance lookups (NSRDB, PVWatts). +NREL_API_KEY= + +# Optional contact email some NREL endpoints request alongside the key. +NREL_API_EMAIL= diff --git a/.gitignore b/.gitignore index 805bf61..ce570aa 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,12 @@ tile_cache/ *.duckdb +*.duckdb.wal .venv/ __pycache__/ *.pyc .DS_Store + +# Secrets — never commit. Real values live only in .env (untracked). +.env +.env.* +!.env.example diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..d67809c --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,84 @@ +# SolarGrader — Engineering Guidelines for Claude Code + +Score every PA home (~5M) for solar sales potential from free public data, then surface +A/A+ homes to reps with optimized routes. **Hard constraint: $0 budget** — every tool and +data source must be free or open source. Background context lives in `memory/` (project +overview, tech stack, pipeline design); read it, don't duplicate it here. + +These principles are binding. When one conflicts with a quick fix, the principle wins — +or stop and flag the tradeoff. Don't silently take the easy path. + +## 1. Design with structure — classes where there's state, functions where there isn't + +Look for the pattern before writing code; build clean, named abstractions that scale to 5M +rows. But this is a numerical/geospatial pipeline, so **match structure to the work**: + +- **Pure functions** for stateless transforms (tile lookup, RANSAC plane fit, pitch/azimuth, + irradiance calc). Easy to test in isolation, no hidden state. This is most of the pipeline. +- **Classes** for genuine state or a real interface: a `TileProcessor`/`Pipeline` orchestrator, + a `Config` dataclass (replace the module-level constants in `test_pipeline.py`), a DB writer + that owns the DuckDB connection, a swappable irradiance source. +- Don't wrap stateless math in a class for ceremony's sake — it makes code harder to test, not + easier. "Simple over complex" (see `memory/tech_stack.md`). +- One responsibility per unit. If a function spans tile I/O **and** RANSAC **and** DB writes, + split it. `test_pipeline.py` is a script today; production code gets factored into a package. + +## 2. Never hardcode secrets — environment variables only + +- API keys / passwords / tokens are **never** inline literals. Read from `os.environ`. +- Real values live only in `.env` (gitignored). `.env.example` documents required vars with + empty values. Keep them in sync — add a var to one, add it to the other. +- **Fail loud:** `os.environ["NREL_API_KEY"]` (KeyError if missing), never + `os.environ.get("NREL_API_KEY", "")`. A defaulted secret is a leaked secret. +- Concrete case now: the NREL NSRDB irradiance source needs a free key (`NREL_API_KEY`). + Wire it through env from day one — don't paste it into a notebook "just to test." +- Never log, print, or commit a secret's value. Before any commit, scan the diff for keys. + +## 3. Minimize USD cost — it's the whole point of this project + +The reason this pipeline exists is to avoid Google Solar API's ~$25k. Protect that: + +- **Free data sources only** (USGS 3DEP, MS Building Footprints, NSRDB/PVGIS, OSM, PASDA). + Before adding any source/service, confirm it's free at PA scale and say so. +- **Stream, don't store.** Process LiDAR tile-by-tile and delete each after — the user does + not have 2TB. Keep working-set storage bounded (~25GB). See `memory/pipeline_design.md`. +- **Cache to avoid re-fetching** (tiles already cache in `tile_cache/`). Don't re-download or + re-compute what's on disk. +- **Think before brute force.** 5M homes × wasteful per-row work is real wall-clock and, on + paid compute, real money. Vectorize with numpy/pandas; batch DB writes; design for the + embarrassingly-parallel tile model (independent tiles, resumable county-by-county). +- If something is free locally but costs money at scale, say so up front. + +## 4. Evaluate security for the operating domain + +Threat model shifts by layer — check the one you're in: + +- **Pipeline / data ingest:** treat external inputs (USGS/OSM/NSRDB responses, LiDAR files) + as untrusted. Validate before use; never `eval`/`exec` fetched content; pin the schema you + expect. Don't pass unsanitized strings into shell/PDAL pipeline JSON. +- **Database (DuckDB → maybe Postgres):** **parameterized queries only**, never string-built + SQL with external values. This applies the moment any user/address input reaches a query. +- **Future web layer (FastAPI + Leaflet):** guard XSS (escape/encode all user-facing output), + SQL injection (parameterized), CSRF on state-changing routes, authn/authz on rep data, + rate-limit public endpoints. Address PII (5M home records) is sensitive — scope access. +- Validate at trust boundaries; fail closed. + +## 5. Additional standards + +- **Type hints + docstrings** on functions/classes; docstring says *what & why*, not line-by-line. +- **Testable & reproducible:** pure functions over hidden globals; keep determinism knobs + explicit (RANSAC iterations/threshold/seed live in config). Validate accuracy on a known tile + before scaling — never trust an unvalidated run across 5M homes. +- **Fail loud, recover gracefully:** raise on programmer error; for a bad tile, log and skip so + a county run isn't lost. Long runs must be interruptible and resumable. +- **Reuse before adding deps** — the user keeps the stack lean (`memory/tech_stack.md`). Justify + any new dependency; prefer what's already imported. +- **Leave the campsite cleaner:** small, focused changes; match surrounding style; no dead code. + +## Working with this user + +- **Explain before running.** Narrate what a command does and why before executing; no + rapid-fire reruns. (`memory/` → explain-before-running) +- **The user runs `test_pipeline.py` himself** to watch the logs/density output. Don't run it + for him — make the change, explain it, hand it back. (`memory/` → user-runs-the-pipeline) +- Prefer simple over complex. Don't add servers, Docker, or dependencies unless necessary. From 50d1e34737949b6fd569c6a261ec80ec8538ab2a Mon Sep 17 00:00:00 2001 From: WagnerJ-Dev Date: Sat, 27 Jun 2026 13:32:34 -0400 Subject: [PATCH 03/13] Dual-grade solar scoring: residential cap + max-roof potential, validated energy model --- test_pipeline.py | 143 ++++++++++++++++++++++++++++++----------------- 1 file changed, 92 insertions(+), 51 deletions(-) diff --git a/test_pipeline.py b/test_pipeline.py index 2190008..d1dd5ae 100644 --- a/test_pipeline.py +++ b/test_pipeline.py @@ -54,6 +54,7 @@ PANEL_PACKING = 0.90 # fraction of the setback area that packs into whole modules SYSTEM_LOSSES = 0.86 # inverter + wiring + soiling (~14% total loss) MAX_AZIMUTH_OFFSET = 120 # skip planes facing >this many degrees from due south +MAX_SYSTEM_KW = 10.0 # cap for a typical sellable residential install (DC kW) PANEL_AREA_M2 = PANEL_WIDTH_M * PANEL_HEIGHT_M # m² per module MODULE_EFF = PANEL_WATTS / (PANEL_AREA_M2 * 1000.0) # derived STC efficiency @@ -375,18 +376,17 @@ def panels_on_plane(plane_area_m2): def calculate_annual_kwh(lat, lon, tmy, planes): """ - Estimate realistic annual production by packing whole panels onto each viable - roof plane and running pvlib's irradiance model for that plane's orientation. - Returns (total_kwh, best_plane, system_kw_dc, total_panels). + Pack whole panels onto each viable roof plane (pvlib irradiance per plane), + then report TWO systems: + - max: every panel that fits the roof (upside potential) + - residential: the best-producing panels up to MAX_SYSTEM_KW (sellable now) + Returns a dict with kWh / kW / panel counts for both, plus the best plane, + or None if no viable (south-ish, panel-fitting) planes exist. """ location = pvlib.location.Location(latitude=lat, longitude=lon, tz="Etc/GMT+5") solar_pos = location.get_solarposition(tmy.index) - total_kwh = 0.0 - total_panels = 0 - best_plane = None - best_kwh = 0.0 - + viable = [] # (per_panel_kwh, n_panels, plane), best-producing first later for plane in planes: # Skip planes facing too far from due south to be worth paneling if abs(plane["azimuth_deg"] - 180) > MAX_AZIMUTH_OFFSET: @@ -396,8 +396,6 @@ def calculate_annual_kwh(lat, lon, tmy, planes): if n_panels <= 0: continue - panel_area = n_panels * PANEL_AREA_M2 - poa = pvlib.irradiance.get_total_irradiance( surface_tilt=plane["tilt_deg"], surface_azimuth=plane["azimuth_deg"], @@ -407,19 +405,37 @@ def calculate_annual_kwh(lat, lon, tmy, planes): dni=tmy["dni"], dhi=tmy["dhi"], ) + # Annual production of ONE panel on this plane (kWh/yr) + poa_annual = float(poa["poa_global"].fillna(0).sum() / 1000) # kWh/m²/yr + per_panel_kwh = poa_annual * PANEL_AREA_M2 * MODULE_EFF * SYSTEM_LOSSES + viable.append((per_panel_kwh, n_panels, plane)) - # Annual plane-of-array irradiance (kWh/m²/yr) × panel area × efficiency - poa_annual = float(poa["poa_global"].fillna(0).sum() / 1000) - plane_kwh = poa_annual * panel_area * MODULE_EFF * SYSTEM_LOSSES - - total_kwh += plane_kwh - total_panels += n_panels - if plane_kwh > best_kwh: - best_kwh = plane_kwh - best_plane = plane + if not viable: + return None - system_kw_dc = total_panels * PANEL_WATTS / 1000.0 - return total_kwh, best_plane, system_kw_dc, total_panels + # Max system: every panel that physically fits the roof + max_panels = sum(n for _, n, _ in viable) + max_kwh = sum(ppk * n for ppk, n, _ in viable) + max_kw = max_panels * PANEL_WATTS / 1000.0 + + # Residential system: fill the best-producing panels first, up to the kW cap + cap_panels = int(MAX_SYSTEM_KW * 1000 / PANEL_WATTS) + res_panels = 0 + res_kwh = 0.0 + for per_panel_kwh, n_panels, _ in sorted(viable, key=lambda v: v[0], reverse=True): + take = min(n_panels, cap_panels - res_panels) + if take <= 0: + break + res_panels += take + res_kwh += take * per_panel_kwh + res_kw = res_panels * PANEL_WATTS / 1000.0 + + best_plane = max(viable, key=lambda v: v[0])[2] + return { + "res_kwh": res_kwh, "res_kw": res_kw, "res_panels": res_panels, + "max_kwh": max_kwh, "max_kw": max_kw, "max_panels": max_panels, + "best_plane": best_plane, + } # ── Step 7: Grading ─────────────────────────────────────────────────────────── @@ -452,6 +468,19 @@ def grade_home(annual_kwh, total_area_m2, primary_azimuth_deg): else: return "D", score +def grade_potential(max_kw): + """ + Tier a roof by its MAXIMUM installable capacity (kW DC) — the upside beyond a + standard residential system (battery, EV, future expansion, or commercial). + Separate from the residential sales grade, which is about the sellable system. + """ + if max_kw >= 30: return "P1" # very high — large/commercial-scale roof + elif max_kw >= 20: return "P2" + elif max_kw >= 12: return "P3" + elif max_kw >= 7: return "P4" + else: return "P5" # limited — barely fits a full residential set + + # ── Step 8: Database ────────────────────────────────────────────────────────── def setup_db(path): @@ -469,11 +498,15 @@ def setup_db(path): primary_tilt_deg DOUBLE, primary_azimuth_deg DOUBLE, roof_plane_count INTEGER, - annual_kwh_potential DOUBLE, - system_kw_dc DOUBLE, - panel_count INTEGER, + res_annual_kwh DOUBLE, + res_system_kw DOUBLE, + res_panel_count INTEGER, solar_score INTEGER, solar_grade VARCHAR(2), + max_annual_kwh DOUBLE, + max_system_kw DOUBLE, + max_panel_count INTEGER, + potential_grade VARCHAR(2), processed_at TIMESTAMP ) """) @@ -483,7 +516,7 @@ def setup_db(path): def save_result(con, building, result): con.execute( """ - INSERT OR REPLACE INTO homes VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) + INSERT OR REPLACE INTO homes VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) """, [ building["osm_id"], @@ -494,11 +527,15 @@ def save_result(con, building, result): result["primary_tilt"], result["primary_azimuth"], result["plane_count"], - result["annual_kwh"], - result["system_kw_dc"], - result["panel_count"], + result["res_kwh"], + result["res_kw"], + result["res_panels"], result["score"], result["grade"], + result["max_kwh"], + result["max_kw"], + result["max_panels"], + result["potential_grade"], pd.Timestamp.now(), ], ) @@ -573,17 +610,17 @@ def main(): counts["skipped_no_planes"] += 1 continue - # Solar calc - annual_kwh, best_plane, system_kw, n_panels = calculate_annual_kwh( - building["lat"], building["lon"], tmy, planes - ) - if annual_kwh == 0 or best_plane is None or n_panels == 0: + # Solar calc — returns both a residential-capped and a max-roof system + sysres = calculate_annual_kwh(building["lat"], building["lon"], tmy, planes) + if sysres is None or sysres["res_panels"] == 0: counts["skipped_no_planes"] += 1 continue - # Grade + # Grade: residential = sellable system (lead quality); potential = max roof + best_plane = sysres["best_plane"] total_area = sum(p["area_m2"] for p in planes) - grade, score = grade_home(annual_kwh, total_area, best_plane["azimuth_deg"]) + grade, score = grade_home(sysres["res_kwh"], total_area, best_plane["azimuth_deg"]) + potential_grade = grade_potential(sysres["max_kw"]) save_result( con, @@ -593,9 +630,13 @@ def main(): "primary_tilt": best_plane["tilt_deg"], "primary_azimuth": best_plane["azimuth_deg"], "plane_count": len(planes), - "annual_kwh": annual_kwh, - "system_kw_dc": system_kw, - "panel_count": n_panels, + "res_kwh": sysres["res_kwh"], + "res_kw": sysres["res_kw"], + "res_panels": sysres["res_panels"], + "max_kwh": sysres["max_kwh"], + "max_kw": sysres["max_kw"], + "max_panels": sysres["max_panels"], + "potential_grade": potential_grade, "grade": grade, "score": score, }, @@ -625,20 +666,20 @@ def main(): if n: print(f" {g:2s} {'█' * n} ({n})") - print(f"\n Top 10 A/A+ leads:") + print(f"\n Top 10 A/A+ leads (res = sellable ~{MAX_SYSTEM_KW:.0f}kW system, max = full-roof upside):") df = con.execute(""" SELECT lat, lon, - ROUND(annual_kwh_potential) AS kwh_yr, - ROUND(system_kw_dc, 1) AS kw_dc, - panel_count AS panels, - ROUND(footprint_m2) AS osm_m2, - ROUND(usable_roof_area_m2) AS roof_m2, + solar_grade AS grade, + ROUND(res_system_kw, 1) AS res_kw, + ROUND(res_annual_kwh) AS res_kwh, + potential_grade AS pot, + ROUND(max_system_kw, 1) AS max_kw, + ROUND(max_annual_kwh) AS max_kwh, ROUND(primary_tilt_deg) AS tilt, - ROUND(primary_azimuth_deg) AS azimuth, - solar_grade AS grade + ROUND(primary_azimuth_deg) AS azimuth FROM homes WHERE solar_grade IN ('A+', 'A') - ORDER BY annual_kwh_potential DESC + ORDER BY res_annual_kwh DESC, max_annual_kwh DESC LIMIT 10 """).fetchdf() @@ -651,13 +692,13 @@ def main(): # Realism check: specific yield (kWh per kW). PA should land ~1,100–1,300. print(f"\n Realism check (kWh per kW installed — PA expected ~1,100–1,300):") yields = con.execute(""" - SELECT annual_kwh_potential / system_kw_dc AS specific_yield, system_kw_dc - FROM homes WHERE system_kw_dc > 0 + SELECT res_annual_kwh / res_system_kw AS specific_yield, res_system_kw + FROM homes WHERE res_system_kw > 0 """).fetchdf() if not yields.empty: med_yield = float(yields["specific_yield"].median()) - med_kw = float(yields["system_kw_dc"].median()) - print(f" Median specific yield: ~{med_yield:,.0f} kWh/kW (median system: {med_kw:.1f} kW DC)") + med_kw = float(yields["res_system_kw"].median()) + print(f" Median specific yield: ~{med_yield:,.0f} kWh/kW (median res system: {med_kw:.1f} kW DC)") if 950 <= med_yield <= 1450: print(" OK — production is in the physically realistic range.") else: From 4d1a9b47b251c06f154b339dcb3f931077ea5843 Mon Sep 17 00:00:00 2001 From: WagnerJ-Dev Date: Sat, 27 Jun 2026 13:53:20 -0400 Subject: [PATCH 04/13] LiDAR skyline shading: per-roof horizon, self-shading excluded, shading in grade --- test_pipeline.py | 137 +++++++++++++++++++++++++++++++++++++---------- 1 file changed, 110 insertions(+), 27 deletions(-) diff --git a/test_pipeline.py b/test_pipeline.py index d1dd5ae..aa2d30f 100644 --- a/test_pipeline.py +++ b/test_pipeline.py @@ -55,6 +55,8 @@ SYSTEM_LOSSES = 0.86 # inverter + wiring + soiling (~14% total loss) MAX_AZIMUTH_OFFSET = 120 # skip planes facing >this many degrees from due south MAX_SYSTEM_KW = 10.0 # cap for a typical sellable residential install (DC kW) +SHADE_RADIUS_M = 50.0 # search radius for shading obstructions (trees, neighbours) +HORIZON_BINS = 72 # azimuth bins for the roof skyline (5° each) PANEL_AREA_M2 = PANEL_WIDTH_M * PANEL_HEIGHT_M # m² per module MODULE_EFF = PANEL_WATTS / (PANEL_AREA_M2 * 1000.0) # derived STC efficiency @@ -374,19 +376,70 @@ def panels_on_plane(plane_area_m2): return int(usable_area // PANEL_AREA_M2) -def calculate_annual_kwh(lat, lon, tmy, planes): +def compute_horizon(tile_points, observer_xyz, exclude_poly=None, + radius_m=SHADE_RADIUS_M, n_bins=HORIZON_BINS): + """ + Build the roof's skyline: the maximum obstruction elevation (degrees above the + roof) in each compass direction, using surrounding LiDAR points — trees, + neighbouring buildings — within radius_m. Points inside exclude_poly (the + building's OWN footprint, in tile CRS) are dropped so a roof can't shade + itself. Returns an array indexed by azimuth bin (0=N, clockwise); all zeros = + wide-open sky. + """ + pts, X, Y = tile_points + ox, oy, oz = observer_xyz + horizon = np.zeros(n_bins) + + m = (X >= ox - radius_m) & (X <= ox + radius_m) & (Y >= oy - radius_m) & (Y <= oy + radius_m) + if not m.any(): + return horizon + + sx, sy = X[m], Y[m] + sz = np.asarray(pts["Z"], dtype=float)[m] + # Drop the building's own roof so it can't count as an obstruction to itself + if exclude_poly is not None: + outside = ~shapely.contains_xy(exclude_poly, sx, sy) + sx, sy, sz = sx[outside], sy[outside], sz[outside] + + dx = sx - ox + dy = sy - oy + dz = sz - oz + horiz = np.sqrt(dx * dx + dy * dy) + # Keep points within radius, not right on top of us, and ABOVE the roof + keep = (horiz <= radius_m) & (horiz > 1.0) & (dz > 0) + if not keep.any(): + return horizon + + dx, dy, dz, horiz = dx[keep], dy[keep], dz[keep], horiz[keep] + elev = np.degrees(np.arctan2(dz, horiz)) # elevation angle + az = (np.degrees(np.arctan2(dx, dy)) + 360) % 360 # 0=N, 90=E (compass) + bins = (az / (360.0 / n_bins)).astype(int) % n_bins + np.maximum.at(horizon, bins, elev) # tallest per direction + return horizon + + +def calculate_annual_kwh(lat, lon, tmy, planes, horizon=None): """ Pack whole panels onto each viable roof plane (pvlib irradiance per plane), then report TWO systems: - max: every panel that fits the roof (upside potential) - residential: the best-producing panels up to MAX_SYSTEM_KW (sellable now) - Returns a dict with kWh / kW / panel counts for both, plus the best plane, - or None if no viable (south-ish, panel-fitting) planes exist. + If a horizon (skyline) is given, the direct beam is removed during hours when + the sun sits below the skyline, and shade_loss_pct reports the production lost. + Returns a dict (incl. shade_loss_pct) or None if no viable planes exist. """ location = pvlib.location.Location(latitude=lat, longitude=lon, tz="Etc/GMT+5") solar_pos = location.get_solarposition(tmy.index) - viable = [] # (per_panel_kwh, n_panels, plane), best-producing first later + # Per-hour direct-beam shading mask: sun up but below the roof's skyline + sun_blocked = None + if horizon is not None: + sun_elev = 90.0 - solar_pos["apparent_zenith"].values + sun_az = solar_pos["azimuth"].values + sun_bins = (sun_az / (360.0 / len(horizon))).astype(int) % len(horizon) + sun_blocked = (sun_elev > 0) & (sun_elev <= horizon[sun_bins]) + + viable = [] # (per_panel_shaded, per_panel_unshaded, n_panels, plane) for plane in planes: # Skip planes facing too far from due south to be worth paneling if abs(plane["azimuth_deg"] - 180) > MAX_AZIMUTH_OFFSET: @@ -405,42 +458,52 @@ def calculate_annual_kwh(lat, lon, tmy, planes): dni=tmy["dni"], dhi=tmy["dhi"], ) - # Annual production of ONE panel on this plane (kWh/yr) - poa_annual = float(poa["poa_global"].fillna(0).sum() / 1000) # kWh/m²/yr - per_panel_kwh = poa_annual * PANEL_AREA_M2 * MODULE_EFF * SYSTEM_LOSSES - viable.append((per_panel_kwh, n_panels, plane)) + glob = poa["poa_global"].fillna(0).values + unshaded_annual = glob.sum() / 1000 # kWh/m²/yr with full sun + if sun_blocked is not None: + # Remove the beam component during blocked hours (diffuse still counts) + direct = poa["poa_direct"].fillna(0).values + shaded_annual = (glob - direct * sun_blocked).sum() / 1000 + else: + shaded_annual = unshaded_annual + + ppk_shaded = shaded_annual * PANEL_AREA_M2 * MODULE_EFF * SYSTEM_LOSSES + ppk_unshaded = unshaded_annual * PANEL_AREA_M2 * MODULE_EFF * SYSTEM_LOSSES + viable.append((ppk_shaded, ppk_unshaded, n_panels, plane)) if not viable: return None # Max system: every panel that physically fits the roof - max_panels = sum(n for _, n, _ in viable) - max_kwh = sum(ppk * n for ppk, n, _ in viable) + max_panels = sum(n for _, _, n, _ in viable) + max_kwh = sum(s * n for s, _, n, _ in viable) + max_unshaded = sum(u * n for _, u, n, _ in viable) max_kw = max_panels * PANEL_WATTS / 1000.0 + shade_loss_pct = 100.0 * (1 - max_kwh / max_unshaded) if max_unshaded > 0 else 0.0 # Residential system: fill the best-producing panels first, up to the kW cap cap_panels = int(MAX_SYSTEM_KW * 1000 / PANEL_WATTS) res_panels = 0 res_kwh = 0.0 - for per_panel_kwh, n_panels, _ in sorted(viable, key=lambda v: v[0], reverse=True): + for ppk_shaded, _, n_panels, _ in sorted(viable, key=lambda v: v[0], reverse=True): take = min(n_panels, cap_panels - res_panels) if take <= 0: break res_panels += take - res_kwh += take * per_panel_kwh + res_kwh += take * ppk_shaded res_kw = res_panels * PANEL_WATTS / 1000.0 - best_plane = max(viable, key=lambda v: v[0])[2] + best_plane = max(viable, key=lambda v: v[0])[3] return { "res_kwh": res_kwh, "res_kw": res_kw, "res_panels": res_panels, "max_kwh": max_kwh, "max_kw": max_kw, "max_panels": max_panels, - "best_plane": best_plane, + "shade_loss_pct": shade_loss_pct, "best_plane": best_plane, } # ── Step 7: Grading ─────────────────────────────────────────────────────────── -def grade_home(annual_kwh, total_area_m2, primary_azimuth_deg): +def grade_home(annual_kwh, total_area_m2, primary_azimuth_deg, shade_loss_pct): score = 0 if annual_kwh >= 12000: score += 50 @@ -457,8 +520,10 @@ def grade_home(annual_kwh, total_area_m2, primary_azimuth_deg): elif offset <= 45: score += 14 elif offset <= 75: score += 7 - # Shade factor placeholder (always 1.0 until shade analysis is implemented) - score += 10 + # Shade credit from the LiDAR skyline analysis (less shading = better lead) + if shade_loss_pct <= 5: score += 10 + elif shade_loss_pct <= 15: score += 7 + elif shade_loss_pct <= 30: score += 4 if score >= 85: return "A+", score elif score >= 70: return "A", score @@ -507,6 +572,7 @@ def setup_db(path): max_system_kw DOUBLE, max_panel_count INTEGER, potential_grade VARCHAR(2), + shade_loss_pct DOUBLE, processed_at TIMESTAMP ) """) @@ -516,7 +582,7 @@ def setup_db(path): def save_result(con, building, result): con.execute( """ - INSERT OR REPLACE INTO homes VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) + INSERT OR REPLACE INTO homes VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) """, [ building["osm_id"], @@ -536,6 +602,7 @@ def save_result(con, building, result): result["max_kw"], result["max_panels"], result["potential_grade"], + result["shade_loss_pct"], pd.Timestamp.now(), ], ) @@ -610,8 +677,16 @@ def main(): counts["skipped_no_planes"] += 1 continue + # Skyline shading from surrounding trees/buildings (observer = roof centroid), + # excluding the building's own footprint so it can't shade itself + observer = points.mean(axis=0) + own_footprint = shapely_transform(transformer.transform, building["geometry"]).buffer(1.0) + horizon = compute_horizon(tile_points, observer, own_footprint) + # Solar calc — returns both a residential-capped and a max-roof system - sysres = calculate_annual_kwh(building["lat"], building["lon"], tmy, planes) + sysres = calculate_annual_kwh( + building["lat"], building["lon"], tmy, planes, horizon + ) if sysres is None or sysres["res_panels"] == 0: counts["skipped_no_planes"] += 1 continue @@ -619,7 +694,9 @@ def main(): # Grade: residential = sellable system (lead quality); potential = max roof best_plane = sysres["best_plane"] total_area = sum(p["area_m2"] for p in planes) - grade, score = grade_home(sysres["res_kwh"], total_area, best_plane["azimuth_deg"]) + grade, score = grade_home( + sysres["res_kwh"], total_area, best_plane["azimuth_deg"], sysres["shade_loss_pct"] + ) potential_grade = grade_potential(sysres["max_kw"]) save_result( @@ -637,6 +714,7 @@ def main(): "max_kw": sysres["max_kw"], "max_panels": sysres["max_panels"], "potential_grade": potential_grade, + "shade_loss_pct": sysres["shade_loss_pct"], "grade": grade, "score": score, }, @@ -672,6 +750,7 @@ def main(): solar_grade AS grade, ROUND(res_system_kw, 1) AS res_kw, ROUND(res_annual_kwh) AS res_kwh, + ROUND(shade_loss_pct) AS shade_pct, potential_grade AS pot, ROUND(max_system_kw, 1) AS max_kw, ROUND(max_annual_kwh) AS max_kwh, @@ -689,20 +768,24 @@ def main(): else: print(df.to_string(index=False)) - # Realism check: specific yield (kWh per kW). PA should land ~1,100–1,300. - print(f"\n Realism check (kWh per kW installed — PA expected ~1,100–1,300):") + # Realism check: validate the PHYSICS on an UNSHADED basis (~1,100–1,300 in PA), + # and report real-world shading separately so it doesn't look like a model error. + print(f"\n Realism check (UNSHADED kWh per kW — PA expected ~1,100–1,300):") yields = con.execute(""" - SELECT res_annual_kwh / res_system_kw AS specific_yield, res_system_kw + SELECT (res_annual_kwh / NULLIF(1 - shade_loss_pct / 100.0, 0)) / res_system_kw + AS unshaded_yield, + shade_loss_pct FROM homes WHERE res_system_kw > 0 """).fetchdf() if not yields.empty: - med_yield = float(yields["specific_yield"].median()) - med_kw = float(yields["res_system_kw"].median()) - print(f" Median specific yield: ~{med_yield:,.0f} kWh/kW (median res system: {med_kw:.1f} kW DC)") + med_yield = float(yields["unshaded_yield"].median()) + med_shade = float(yields["shade_loss_pct"].median()) + print(f" Median unshaded yield: ~{med_yield:,.0f} kWh/kW") if 950 <= med_yield <= 1450: - print(" OK — production is in the physically realistic range.") + print(" OK — the panel/energy physics is in the realistic range.") else: print(" WARNING — outside the expected band; the panel/energy model needs review.") + print(f" Median shading loss: {med_shade:.0f}% (real reduction from trees/neighbours)") # Diagnostic: LiDAR point density check (averaged over buildings that had LiDAR) print(f"\n LiDAR density check:") From 31e42fa51417a26cf7beea6dd6efab5b2a6e12e3 Mon Sep 17 00:00:00 2001 From: WagnerJ-Dev Date: Sat, 27 Jun 2026 14:40:10 -0400 Subject: [PATCH 05/13] Streaming multi-tile pipeline: collection filter, resume, retry, accumulate --- pipeline.py | 312 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 312 insertions(+) create mode 100644 pipeline.py diff --git a/pipeline.py b/pipeline.py new file mode 100644 index 0000000..aa5b19b --- /dev/null +++ b/pipeline.py @@ -0,0 +1,312 @@ +""" +Solar Grader — Streaming Multi-Tile Pipeline +Processes every LiDAR tile covering a region, one at a time, accumulating scored +homes into DuckDB. Streaming + tile deletion keeps disk bounded; a tiles_done +table makes the run resumable (re-running skips finished tiles). + +Reuses the VALIDATED scoring functions from test_pipeline.py unchanged — only the +orchestration (tile discovery, accumulate, resume, cleanup) is new here. + +Run with: + source .venv/bin/activate + python pipeline.py +""" + +import os +import time +import warnings + +import duckdb +import pandas as pd +import requests +from pyproj import Transformer +from shapely.ops import transform as shapely_transform + +# The scoring algorithm — imported as-is from the validated test pipeline +from test_pipeline import ( + download_tile, + get_tile_srs, + load_tile_points, + get_buildings_from_osm, + get_tmy_data, + clip_lidar_to_building, + extract_roof_planes, + compute_horizon, + calculate_annual_kwh, + grade_home, + grade_potential, + save_result, +) + +warnings.filterwarnings("ignore") + +# ── Configuration ───────────────────────────────────────────────────────────── + +# First validation target: West Chester, PA. +# Widen this to a county once the streaming/resume mechanics are confirmed. +REGION_BBOX = (-75.620, 39.945, -75.595, 39.970) # (west, south, east, north) + +# CRITICAL: TNM serves overlapping LiDAR surveys of wildly different vintage/quality +# for the same area. Pin to ONE modern, high-density collection so results are +# consistent — mixing a 2024 survey with 2006 data corrupts the scores. Tiles whose +# title doesn't contain this string are dropped. (PA_17County_D24 = 2024, ~35 pts/m², +# covers the Chester County rollout. Full-PA would map the best collection per region.) +PREFERRED_COLLECTION = "PA_17County_D24" + +# Cap tiles for a fast first validation; set to None to process the whole region. +MAX_TILES = 3 + +DB_PATH = "solar_grader.duckdb" # accumulating production DB (separate from the test) +TILE_CACHE_DIR = "tile_stream" # tiles downloaded here, then deleted after processing +DELETE_TILES_AFTER = True # the storage-capping mechanic; resume covers re-runs + + +# ── Tile discovery ──────────────────────────────────────────────────────────── + +def find_all_tiles(bbox, page=50): + """Query USGS TNM for ALL LiDAR tiles covering bbox (paginated, de-duped).""" + west, south, east, north = bbox + url = "https://tnmaccess.nationalmap.gov/api/v1/products" + tiles, offset = {}, 0 + while True: + params = { + "datasets": "Lidar Point Cloud (LPC)", + "bbox": f"{west},{south},{east},{north}", + "prodFormats": "LAZ", + "max": page, + "offset": offset, + } + r = requests.get(url, params=params, timeout=60) + r.raise_for_status() + data = r.json() + items = data.get("items", []) + for it in items: + dl = it.get("downloadURL") + if dl: + tiles[dl] = it + offset += len(items) + if not items or offset >= data.get("total", 0): + break + return list(tiles.values()) + + +def tile_id_of(tile): + """Stable id for a tile = its LAZ filename.""" + return tile["downloadURL"].split("/")[-1].split("?")[0] + + +def with_retries(fn, *args, attempts=4, base_delay=5, **kwargs): + """Call fn with exponential backoff on transient failures (e.g. Overpass 504, + PVGIS hiccups). Re-raises only after the final attempt.""" + for i in range(attempts): + try: + return fn(*args, **kwargs) + except Exception as e: + if i == attempts - 1: + raise + wait = base_delay * (2 ** i) + print(f" {type(e).__name__}: {e}\n retry {i + 1}/{attempts - 1} in {wait}s...") + time.sleep(wait) + + +# ── Database (accumulating, resumable) ──────────────────────────────────────── + +def setup_db(path): + con = duckdb.connect(path) + # Same homes schema as the test pipeline, but NEVER dropped — we accumulate. + con.execute(""" + CREATE TABLE IF NOT EXISTS homes ( + osm_id BIGINT PRIMARY KEY, + lat DOUBLE, + lon DOUBLE, + footprint_m2 DOUBLE, + usable_roof_area_m2 DOUBLE, + primary_tilt_deg DOUBLE, + primary_azimuth_deg DOUBLE, + roof_plane_count INTEGER, + res_annual_kwh DOUBLE, + res_system_kw DOUBLE, + res_panel_count INTEGER, + solar_score INTEGER, + solar_grade VARCHAR(2), + max_annual_kwh DOUBLE, + max_system_kw DOUBLE, + max_panel_count INTEGER, + potential_grade VARCHAR(2), + shade_loss_pct DOUBLE, + processed_at TIMESTAMP + ) + """) + # Resume ledger: which tiles are already done. + con.execute(""" + CREATE TABLE IF NOT EXISTS tiles_done ( + tile_id VARCHAR PRIMARY KEY, + n_scored INTEGER, + processed_at TIMESTAMP + ) + """) + return con + + +# ── Per-building scoring (mirrors test_pipeline.main's inner loop) ───────────── + +def score_building(building, tile_points, transformer, tmy): + """Full scoring chain for one building. Returns a save_result payload or None.""" + points = clip_lidar_to_building(tile_points, building["geometry"], transformer) + if points is None: + return None + planes = extract_roof_planes(points) + if not planes: + return None + + observer = points.mean(axis=0) + own_footprint = shapely_transform(transformer.transform, building["geometry"]).buffer(1.0) + horizon = compute_horizon(tile_points, observer, own_footprint) + + sysres = calculate_annual_kwh(building["lat"], building["lon"], tmy, planes, horizon) + if sysres is None or sysres["res_panels"] == 0: + return None + + best_plane = sysres["best_plane"] + total_area = sum(p["area_m2"] for p in planes) + grade, score = grade_home( + sysres["res_kwh"], total_area, best_plane["azimuth_deg"], sysres["shade_loss_pct"] + ) + return { + "usable_area_m2": total_area, + "primary_tilt": best_plane["tilt_deg"], + "primary_azimuth": best_plane["azimuth_deg"], + "plane_count": len(planes), + "res_kwh": sysres["res_kwh"], + "res_kw": sysres["res_kw"], + "res_panels": sysres["res_panels"], + "max_kwh": sysres["max_kwh"], + "max_kw": sysres["max_kw"], + "max_panels": sysres["max_panels"], + "potential_grade": grade_potential(sysres["max_kw"]), + "shade_loss_pct": sysres["shade_loss_pct"], + "grade": grade, + "score": score, + } + + +def tile_latlon_bbox(tile_points, tile_srs): + """Derive the tile's lon/lat bbox from its actual point extent (UTM → 4326).""" + _, X, Y = tile_points + inv = Transformer.from_crs(tile_srs, "EPSG:4326", always_xy=True) + xs = [X.min(), X.max(), X.min(), X.max()] + ys = [Y.min(), Y.min(), Y.max(), Y.max()] + lons, lats = inv.transform(xs, ys) + return (min(lons), min(lats), max(lons), max(lats)) + + +# ── Main streaming loop ─────────────────────────────────────────────────────── + +def run(bbox, db_path=DB_PATH, tile_cache=TILE_CACHE_DIR): + print("=" * 60) + print("SOLAR GRADER — Streaming Multi-Tile Pipeline") + print(f"Region bbox: {bbox}") + print("=" * 60) + + print("\nDiscovering LiDAR tiles...") + all_tiles = find_all_tiles(bbox) + tiles = [t for t in all_tiles if PREFERRED_COLLECTION in t.get("title", "")] + dropped = len(all_tiles) - len(tiles) + print(f" {len(all_tiles)} tiles cover this region; {len(tiles)} match " + f"'{PREFERRED_COLLECTION}', {dropped} other-collection tiles dropped.") + if not tiles: + print(" ABORT: no tiles match the preferred collection — refusing to mix qualities.") + return + if MAX_TILES is not None: + tiles = tiles[:MAX_TILES] + print(f" MAX_TILES set — processing only the first {len(tiles)} this run.") + + con = setup_db(db_path) + done = {row[0] for row in con.execute("SELECT tile_id FROM tiles_done").fetchall()} + if done: + print(f" {len(done)} already processed — resuming, will skip those.") + + t_start = time.time() + for i, tile in enumerate(tiles, 1): + tid = tile_id_of(tile) + if tid in done: + print(f"\n[tile {i}/{len(tiles)}] {tid} (skip — already done)") + continue + + print(f"\n[tile {i}/{len(tiles)}] {tid}") + try: + path = download_tile(tile["downloadURL"], tile_cache) + tile_points = load_tile_points(path) + tile_srs = get_tile_srs(path) + transformer = Transformer.from_crs("EPSG:4326", tile_srs, always_xy=True) + + tbbox = tile_latlon_bbox(tile_points, tile_srs) + buildings = with_retries(get_buildings_from_osm, tbbox) + if not buildings: + _mark_done(con, tid, 0) + _cleanup(path) + continue + + # One irradiance fetch per tile (its center) covers all its buildings + cx = (tbbox[0] + tbbox[2]) / 2 + cy = (tbbox[1] + tbbox[3]) / 2 + tmy = with_retries(get_tmy_data, cy, cx) + + n_scored = 0 + for j, building in enumerate(buildings): + print(f"\r scoring {j+1}/{len(buildings)} (scored={n_scored}) ", end="") + result = score_building(building, tile_points, transformer, tmy) + if result: + save_result(con, building, result) + n_scored += 1 + print(f"\r scored {n_scored} of {len(buildings)} buildings in tile. ") + + _mark_done(con, tid, n_scored) + _cleanup(path) + + except Exception as e: + # Isolate failures — one bad tile must not kill the whole run + print(f"\n ERROR on {tid}: {type(e).__name__}: {e} — skipping tile.") + continue + + _summary(con, time.time() - t_start) + con.close() + + +def _mark_done(con, tid, n_scored): + con.execute( + "INSERT OR REPLACE INTO tiles_done VALUES (?, ?, ?)", + [tid, n_scored, pd.Timestamp.now()], + ) + + +def _cleanup(path): + if DELETE_TILES_AFTER and os.path.exists(path): + os.remove(path) + + +def _summary(con, elapsed): + print(f"\n\n{'=' * 60}") + print("PIPELINE SUMMARY") + print(f"{'=' * 60}") + n_homes = con.execute("SELECT COUNT(*) FROM homes").fetchone()[0] + n_tiles = con.execute("SELECT COUNT(*) FROM tiles_done").fetchone()[0] + print(f" Tiles processed: {n_tiles}") + print(f" Homes scored: {n_homes}") + print(f" Total time: {elapsed:.1f}s") + + print(f"\n Grade breakdown:") + rows = con.execute(""" + SELECT solar_grade, COUNT(*) FROM homes GROUP BY solar_grade + """).fetchall() + counts = dict(rows) + for g in ["A+", "A", "B+", "B", "C", "D"]: + n = counts.get(g, 0) + if n: + print(f" {g:2s} {'█' * min(n, 50)} ({n})") + + print(f"\n Results saved to: {DB_PATH}") + + +if __name__ == "__main__": + run(REGION_BBOX) From 6f245e8a3bf46e93105388446cdee57fb770a006 Mon Sep 17 00:00:00 2001 From: WagnerJ-Dev Date: Tue, 30 Jun 2026 21:35:20 -0400 Subject: [PATCH 06/13] working poc --- enrich_addresses.py | 187 ++++++++++++++++++++++++++++++++++++++++++++ make_map.py | 147 ++++++++++++++++++++++++++++++++++ map.html | 48 ++++++++++++ pipeline.py | 2 +- 4 files changed, 383 insertions(+), 1 deletion(-) create mode 100644 enrich_addresses.py create mode 100644 make_map.py create mode 100644 map.html diff --git a/enrich_addresses.py b/enrich_addresses.py new file mode 100644 index 0000000..b4c567a --- /dev/null +++ b/enrich_addresses.py @@ -0,0 +1,187 @@ +""" +Solar Grader — Address Enrichment (Chester County, PA) + +Post-processing step: attaches real street addresses to scored homes by matching +each building centroid to its nearest official address point from Chester County's +PASDA ArcGIS service. Decoupled from scoring — operates on the existing DuckDB, +re-runnable, no LiDAR reprocessing. + +Coverage note: the address-point service is Chester-County-specific. Other counties +have their own PASDA/GIS services; we'd map the right one per region when scaling +beyond Chester. + +Run with: + source .venv/bin/activate + python enrich_addresses.py +""" + +import warnings + +import duckdb +import numpy as np +import pandas as pd +import requests +from scipy.spatial import cKDTree + +warnings.filterwarnings("ignore") + +DB_PATH = "solar_grader.duckdb" +ADDR_LAYER = ( + "https://mapservices.pasda.psu.edu/server/rest/services/" + "pasda/ChesterCounty/MapServer/12/query" +) +MATCH_RADIUS_M = 50.0 # max centroid→address-point distance to accept a match +BBOX_MARGIN_DEG = 0.002 # pad the query bbox so edge homes still find a point + +# Residential filter: exclude institutional (tax-exempt), multi-unit/commercial +# (address has a unit marker), or too-large-to-be-a-house footprints. Heuristic — +# refine later with parcel land-use codes if needed. +RESIDENTIAL_MAX_FOOTPRINT_M2 = 500 +UNIT_MARKERS = (" APT", " STE", " UNIT", " FL ", " #", " BLDG", " RM ") + + +# ── Pull address points from the county service ─────────────────────────────── + +def fetch_address_points(bbox): + """All Chester County address points intersecting bbox, via ArcGIS REST. + Uses an id-only query then batched fetches (robust to the 1000-record cap).""" + west, south, east, north = bbox + geom = { + "geometry": f"{west},{south},{east},{north}", + "geometryType": "esriGeometryEnvelope", + "inSR": "4326", + "spatialRel": "esriSpatialRelIntersects", + } + + ids = requests.get( + ADDR_LAYER, + params={**geom, "where": "1=1", "returnIdsOnly": "true", "f": "json"}, + timeout=90, + ).json().get("objectIds") or [] + print(f" {len(ids)} address points in region; fetching...") + + points = [] + for i in range(0, len(ids), 200): + batch = ids[i:i + 200] + # POST (not GET) so a long objectIds list can't blow the URL length limit + resp = requests.post( + ADDR_LAYER, + data={ + "objectIds": ",".join(map(str, batch)), + "outFields": "ADDR_NUM,ROAD_NAME,FULL_ADDRE,MUNI_NAME,ZIP_CODE_S,TAX_EXEMPT", + "returnGeometry": "true", + "outSR": "4326", + "f": "json", + }, + timeout=90, + ) + resp.raise_for_status() + data = resp.json() + for f in data.get("features", []): + g, a = f.get("geometry"), f.get("attributes", {}) + if g: + points.append((g["x"], g["y"], a)) + return points + + +# ── Local equirectangular projection for fast nearest-neighbour matching ────── + +def _projector(lat0): + cos0 = np.cos(np.radians(lat0)) + return lambda lon, lat: (np.asarray(lon) * 111320.0 * cos0, + np.asarray(lat) * 110540.0) + + +def main(): + try: + con = duckdb.connect(DB_PATH) + except duckdb.IOException: + print("Database is locked — is pipeline.py still running?") + print("Let the scoring run finish, then re-run this (it's a separate step).") + return + homes = con.execute("SELECT osm_id, lat, lon, footprint_m2 FROM homes").fetchdf() + if homes.empty: + print("No homes in the database yet — run pipeline.py first.") + return + print(f"Enriching {len(homes)} scored homes with addresses...") + + bbox = ( + homes["lon"].min() - BBOX_MARGIN_DEG, homes["lat"].min() - BBOX_MARGIN_DEG, + homes["lon"].max() + BBOX_MARGIN_DEG, homes["lat"].max() + BBOX_MARGIN_DEG, + ) + pts = fetch_address_points(bbox) + if not pts: + print("No address points returned — check the region/service.") + return + + lat0 = float(homes["lat"].median()) + proj = _projector(lat0) + ax, ay = proj([p[0] for p in pts], [p[1] for p in pts]) + tree = cKDTree(np.column_stack([ax, ay])) + + hx, hy = proj(homes["lon"].values, homes["lat"].values) + dist, idx = tree.query(np.column_stack([hx, hy]), distance_upper_bound=MATCH_RADIUS_M) + + footprints = dict(zip(homes["osm_id"].values, homes["footprint_m2"].values)) + rows, matched = [], 0 + for osm_id, d, j in zip(homes["osm_id"].values, dist, idx): + osm_id = int(osm_id) + foot = footprints.get(osm_id) + small = foot is None or foot < RESIDENTIAL_MAX_FOOTPRINT_M2 + if np.isinf(d): + # No address match — classify residential on footprint alone + rows.append((osm_id, None, None, None, None, None, None, small)) + continue + matched += 1 + a = pts[j][2] + full = (a.get("FULL_ADDRE") or "") + tax = a.get("TAX_EXEMPT") + has_unit = any(mk in f" {full.upper()}" for mk in UNIT_MARKERS) + is_res = (tax != "Yes") and (not has_unit) and small + rows.append((osm_id, full or None, a.get("ADDR_NUM"), a.get("ROAD_NAME"), + a.get("MUNI_NAME"), a.get("ZIP_CODE_S"), tax, is_res)) + + # Add columns (idempotent) and bulk-update via a join + for col, typ in [("full_address", "VARCHAR"), ("house_number", "VARCHAR"), + ("street", "VARCHAR"), ("city", "VARCHAR"), ("zip_code", "VARCHAR"), + ("tax_exempt", "VARCHAR"), ("is_residential", "BOOLEAN")]: + try: + con.execute(f"ALTER TABLE homes ADD COLUMN {col} {typ}") + except duckdb.Error: + pass # already exists on a re-run + + addr_df = pd.DataFrame(rows, columns=[ + "osm_id", "full_address", "house_number", "street", "city", "zip_code", + "tax_exempt", "is_residential"]) + con.register("addr_df", addr_df) + con.execute(""" + UPDATE homes SET + full_address = a.full_address, house_number = a.house_number, + street = a.street, city = a.city, zip_code = a.zip_code, + tax_exempt = a.tax_exempt, is_residential = a.is_residential + FROM addr_df a + WHERE homes.osm_id = a.osm_id + """) + + # ── Report ──────────────────────────────────────────────────────────────── + total = len(homes) + n_res = con.execute("SELECT COUNT(*) FROM homes WHERE is_residential").fetchone()[0] + print(f"\n Matched {matched} / {total} homes ({100 * matched / total:.0f}%) to an address.") + print(f" Residential (single-family): {n_res} / {total} " + f"({100 * n_res / total:.0f}%) — the rest are institutional/multi-unit/large.") + + print("\n Top RESIDENTIAL A+/A leads:") + df = con.execute(""" + SELECT full_address, city, zip_code, solar_grade AS grade, + ROUND(res_annual_kwh) AS res_kwh, potential_grade AS pot + FROM homes + WHERE solar_grade IN ('A+', 'A') AND is_residential AND full_address IS NOT NULL + ORDER BY res_annual_kwh DESC + LIMIT 10 + """).fetchdf() + print(df.to_string(index=False) if not df.empty else " (none yet)") + con.close() + + +if __name__ == "__main__": + main() diff --git a/make_map.py b/make_map.py new file mode 100644 index 0000000..1531419 --- /dev/null +++ b/make_map.py @@ -0,0 +1,147 @@ +""" +Solar Grader — Leaflet Map Builder + +Reads scored + address-enriched homes from DuckDB and writes a self-contained +map.html (data embedded, no server needed — just open it in a browser). Pins are +colored by solar grade; click one for address, grade, production, and potential. + +Defaults to residential leads only; set INCLUDE_NON_RESIDENTIAL = True to show all. + +Run with: + source .venv/bin/activate + python make_map.py + open map.html +""" + +import json +import warnings + +import duckdb + +warnings.filterwarnings("ignore") + +DB_PATH = "solar_grader.duckdb" +OUT_HTML = "map.html" +INCLUDE_NON_RESIDENTIAL = False +GRADES_SHOWN = ("A+", "A", "B+", "B", "C", "D") # which grades to plot + +GRADE_COLORS = { + "A+": "#1a9850", "A": "#66bd63", "B+": "#a6d96a", + "B": "#fee08b", "C": "#fdae61", "D": "#d73027", +} + + +def load_features(): + con = duckdb.connect(DB_PATH, read_only=True) + cols = {c[0] for c in con.execute("DESCRIBE homes").fetchall()} + if "is_residential" not in cols: + con.close() + raise SystemExit("No address/residential columns — run enrich_addresses.py first.") + + where = ["solar_grade IN ({})".format(",".join(f"'{g}'" for g in GRADES_SHOWN))] + if not INCLUDE_NON_RESIDENTIAL: + where.append("is_residential") + rows = con.execute(f""" + SELECT lat, lon, solar_grade, solar_score, + ROUND(res_annual_kwh) AS res_kwh, ROUND(res_system_kw, 1) AS res_kw, + potential_grade, ROUND(max_system_kw, 1) AS max_kw, + ROUND(shade_loss_pct) AS shade, COALESCE(full_address, '(no address)') AS addr, + COALESCE(city, '') AS city + FROM homes + WHERE {' AND '.join(where)} + """).fetchdf() + con.close() + + features = [] + for r in rows.itertuples(index=False): + features.append({ + "type": "Feature", + "geometry": {"type": "Point", "coordinates": [r.lon, r.lat]}, + "properties": { + "grade": r.solar_grade, "score": int(r.solar_score), + "res_kwh": None if r.res_kwh != r.res_kwh else int(r.res_kwh), + "res_kw": r.res_kw, "pot": r.potential_grade, "max_kw": r.max_kw, + "shade": None if r.shade != r.shade else int(r.shade), + "addr": r.addr, "city": r.city, + }, + }) + return rows, features + + +def build_html(rows, features): + center_lat = float(rows["lat"].mean()) + center_lon = float(rows["lon"].mean()) + geojson = json.dumps({"type": "FeatureCollection", "features": features}) + colors = json.dumps(GRADE_COLORS) + + # Legend rows + legend = "".join( + f'
{g}
' + for g, c in GRADE_COLORS.items() + ) + + return f""" + +Solar Grader — Leads + + + + +
+""" + + +def main(): + rows, features = load_features() + if not features: + print("No homes to map — run pipeline.py and enrich_addresses.py first.") + return + html = build_html(rows, features) + with open(OUT_HTML, "w") as f: + f.write(html) + kind = "all" if INCLUDE_NON_RESIDENTIAL else "residential" + print(f"Wrote {OUT_HTML} with {len(features)} {kind} leads.") + print(f"Open it with: open {OUT_HTML}") + + +if __name__ == "__main__": + main() diff --git a/map.html b/map.html new file mode 100644 index 0000000..b7fac34 --- /dev/null +++ b/map.html @@ -0,0 +1,48 @@ + + +Solar Grader — Leads + + + + +
+ \ No newline at end of file diff --git a/pipeline.py b/pipeline.py index aa5b19b..a08b132 100644 --- a/pipeline.py +++ b/pipeline.py @@ -54,7 +54,7 @@ PREFERRED_COLLECTION = "PA_17County_D24" # Cap tiles for a fast first validation; set to None to process the whole region. -MAX_TILES = 3 +MAX_TILES = None DB_PATH = "solar_grader.duckdb" # accumulating production DB (separate from the test) TILE_CACHE_DIR = "tile_stream" # tiles downloaded here, then deleted after processing From 3c0d613b2c3cc62715729b1eb7e3276d40ed62c4 Mon Sep 17 00:00:00 2001 From: WagnerJ-Dev Date: Tue, 30 Jun 2026 23:30:44 -0400 Subject: [PATCH 07/13] County footprints, residential filter, map, rockyweb tile-discovery fallback --- pipeline.py | 129 ++++++++++++++++++++++++++++++++++++++++- solar-pipeline-plan.md | 16 ++--- 2 files changed, 135 insertions(+), 10 deletions(-) diff --git a/pipeline.py b/pipeline.py index a08b132..362c0f0 100644 --- a/pipeline.py +++ b/pipeline.py @@ -12,7 +12,9 @@ python pipeline.py """ +import math import os +import re import time import warnings @@ -20,6 +22,7 @@ import pandas as pd import requests from pyproj import Transformer +from shapely.geometry import Polygon from shapely.ops import transform as shapely_transform # The scoring algorithm — imported as-is from the validated test pipeline @@ -56,6 +59,23 @@ # Cap tiles for a fast first validation; set to None to process the whole region. MAX_TILES = None +# Building inventory source: "county" (Chester ArcGIS footprints — most complete +# here, no Overpass rate limits) or "osm" (Overpass — incomplete, prototyping only). +BUILDING_SOURCE = "county" +COUNTY_FOOTPRINTS = ( + "https://mapservices.pasda.psu.edu/server/rest/services/" + "pasda/ChesterCounty/MapServer/14/query" +) + +# Fallback tile discovery when TNM's spatial query is down: browse the project's +# staged LAZ directories on rockyweb and select tiles by the bbox→tile-name decode. +# (Tile name USGS_LPC_..._18SVK = UTM 18N easting EEE*1000, northing +# 4,000,000 + NNN*1000. Verified against known tiles.) Scoped to this project. +ROCKYWEB_PROJECT = ( + "https://rockyweb.usgs.gov/vdelivery/Datasets/Staged/Elevation/LPC/" + "Projects/PA_17County_D24" +) + DB_PATH = "solar_grader.duckdb" # accumulating production DB (separate from the test) TILE_CACHE_DIR = "tile_stream" # tiles downloaded here, then deleted after processing DELETE_TILES_AFTER = True # the storage-capping mechanic; resume covers re-runs @@ -95,6 +115,100 @@ def tile_id_of(tile): return tile["downloadURL"].split("/")[-1].split("?")[0] +def _covering_tile_codes(bbox): + """The set of 6-digit tile codes (EEENNN) whose 1 km UTM cells cover bbox.""" + fwd = Transformer.from_crs("EPSG:4326", "EPSG:26918", always_xy=True) # UTM 18N + w, s, e, n = bbox + xs, ys = fwd.transform([w, e, w, e], [s, s, n, n]) + e_range = range(int(min(xs) // 1000), int(max(xs) // 1000) + 1) + n_range = range(int((min(ys) - 4_000_000) // 1000), int((max(ys) - 4_000_000) // 1000) + 1) + return {f"{ee}{nn}" for ee in e_range for nn in n_range} + + +def find_tiles_rockyweb(bbox): + """Fallback for when TNM is down: list the project's staged LAZ dirs on + rockyweb and pick tiles covering bbox via the tile-name UTM decode.""" + codes = _covering_tile_codes(bbox) + index = requests.get(f"{ROCKYWEB_PROJECT}/", timeout=60).text + subprojects = sorted(set(re.findall(r'href="(PA_17Co_\d+_D24)/"', index))) + + tiles, seen = [], set() + for sub in subprojects: + laz_url = f"{ROCKYWEB_PROJECT}/{sub}/LAZ/" + listing = requests.get(laz_url, timeout=90).text + for fn in re.findall(r"USGS_LPC_PA_17County_D24_18[A-Z]{3}\d{6}\.laz", listing): + code = fn[-10:-4] # the 6 digits before ".laz" + if code in codes and fn not in seen: + seen.add(fn) + mgrs = fn[:-4].split("_")[-1] + tiles.append({ + "downloadURL": laz_url + fn, + "title": f"USGS Lidar Point Cloud PA_17County_D24 {mgrs}", + "sizeInBytes": 0, + }) + if len(seen) >= len(codes): + break # found every covering tile + return tiles + + +# ── Building inventory sources ──────────────────────────────────────────────── + +def get_buildings_from_county(bbox): + """Chester County building footprints (PASDA MapServer layer 14), bbox-queried. + Returns the same dict shape as get_buildings_from_osm so scoring is unchanged. + Uses an id-only query then batched geometry fetches (robust to the 1000 cap).""" + west, south, east, north = bbox + geom = { + "geometry": f"{west},{south},{east},{north}", + "geometryType": "esriGeometryEnvelope", + "inSR": "4326", + "spatialRel": "esriSpatialRelIntersects", + } + print(" Querying county building footprints...") + ids = requests.get( + COUNTY_FOOTPRINTS, + params={**geom, "where": "1=1", "returnIdsOnly": "true", "f": "json"}, + timeout=90, + ).json().get("objectIds") or [] + + buildings = [] + for i in range(0, len(ids), 200): + batch = ids[i:i + 200] + resp = requests.post( + COUNTY_FOOTPRINTS, + data={"objectIds": ",".join(map(str, batch)), "outFields": "OBJECTID", + "returnGeometry": "true", "outSR": "4326", "f": "json"}, + timeout=90, + ) + resp.raise_for_status() + for f in resp.json().get("features", []): + rings = (f.get("geometry") or {}).get("rings") + if not rings: + continue + try: + poly = Polygon(rings[0]) # exterior ring + if not poly.is_valid or poly.area == 0: + continue + except Exception: + continue + cx, cy = poly.centroid.x, poly.centroid.y + deg2m2 = (111320 ** 2) * math.cos(math.radians(cy)) + buildings.append({ + "osm_id": int(f["attributes"]["OBJECTID"]), # unique building id + "geometry": poly, "lat": cy, "lon": cx, + "footprint_m2": poly.area * deg2m2, + }) + print(f" Found {len(buildings)} buildings.") + return buildings + + +def get_buildings(bbox): + """Dispatch to the configured building source.""" + if BUILDING_SOURCE == "county": + return get_buildings_from_county(bbox) + return get_buildings_from_osm(bbox) + + def with_retries(fn, *args, attempts=4, base_delay=5, **kwargs): """Call fn with exponential backoff on transient failures (e.g. Overpass 504, PVGIS hiccups). Re-raises only after the final attempt.""" @@ -209,13 +323,22 @@ def run(bbox, db_path=DB_PATH, tile_cache=TILE_CACHE_DIR): print("=" * 60) print("\nDiscovering LiDAR tiles...") - all_tiles = find_all_tiles(bbox) + all_tiles = with_retries(find_all_tiles, bbox) + if not all_tiles: + print(" TNM returned 0 (spatial query down) — falling back to direct") + print(" rockyweb tile discovery...") + all_tiles = with_retries(find_tiles_rockyweb, bbox) + print(f" rockyweb found {len(all_tiles)} tiles.") + if not all_tiles: + print(" ABORT: no tiles from TNM or rockyweb. Both may be down — retry later.") + return tiles = [t for t in all_tiles if PREFERRED_COLLECTION in t.get("title", "")] dropped = len(all_tiles) - len(tiles) print(f" {len(all_tiles)} tiles cover this region; {len(tiles)} match " f"'{PREFERRED_COLLECTION}', {dropped} other-collection tiles dropped.") if not tiles: - print(" ABORT: no tiles match the preferred collection — refusing to mix qualities.") + print(f" ABORT: {len(all_tiles)} tiles found but none match " + f"'{PREFERRED_COLLECTION}' — refusing to mix collection qualities.") return if MAX_TILES is not None: tiles = tiles[:MAX_TILES] @@ -241,7 +364,7 @@ def run(bbox, db_path=DB_PATH, tile_cache=TILE_CACHE_DIR): transformer = Transformer.from_crs("EPSG:4326", tile_srs, always_xy=True) tbbox = tile_latlon_bbox(tile_points, tile_srs) - buildings = with_retries(get_buildings_from_osm, tbbox) + buildings = with_retries(get_buildings, tbbox) if not buildings: _mark_done(con, tid, 0) _cleanup(path) diff --git a/solar-pipeline-plan.md b/solar-pipeline-plan.md index a1490c5..b1e759c 100644 --- a/solar-pipeline-plan.md +++ b/solar-pipeline-plan.md @@ -48,15 +48,17 @@ A one-time data pipeline that produces a permanent, per-address solar score for # Use USGS National Map Downloader to enumerate available tiles for PA bounding box ``` -### 1B. Building Footprints (Microsoft) +### 1B. Building Footprints — pluggable source (the key to full coverage) -**What it gives you:** 2D polygon outlines of every building. Used to mask LiDAR points to only those belonging to a specific building. +**What it gives you:** 2D polygon outlines of every building. Used to mask LiDAR points to only those belonging to a specific building. Completeness here = how many houses end up scored, so this is the single biggest coverage lever. -**How to get it:** -- https://github.com/microsoft/buildings -- Pennsylvania file: direct GeoJSON download (~300–400 MB) -- Schema: geometry (polygon) + height estimate -- Coverage: ~5.1M buildings in PA +**Design:** a pluggable `get_buildings(bbox)` with interchangeable backends (scoring/address/grading downstream is identical regardless of source). Backends, from validation → scale: + +1. **OSM Overpass** — used in initial validation. Simple, per-bbox API, but **incomplete** (measured ~2,500 vs the county's 3,391 in West Chester = ~26% of buildings missing) and rate-limits/504s at scale. Fine for prototyping, not for coverage. +2. **County ArcGIS footprint layer** (CURRENT for Chester rollout) — e.g. PASDA Chester County `MapServer/14` Building Footprints, bbox-queryable like the address layer. More complete than OSM, and **no Overpass rate-limit/504 problem**. Best where a county publishes footprints; downside is per-county endpoints. +3. **Microsoft US Building Footprints** (for statewide → other states) — https://github.com/microsoft/USBuildingFootprints, ODbL, ~130M US buildings, one GeoJSON per state (~PA 5M). Quality ≥ OSM per Microsoft. Pattern: download state file once → load into a local spatial index (DuckDB spatial extension) → query per-tile locally. **Uniform national coverage = the "PA then state-by-state" unlock.** + +**Also:** record un-scoreable buildings (no roof / no viable plane) as a low/"not-viable" status instead of dropping them, so every structure appears in the output, not just solar-viable ones. ### 1C. Solar Irradiance (NREL NSRDB) From 288f23f3f4754cc73d45a3e48da0ff5263be8509 Mon Sep 17 00:00:00 2001 From: WagnerJ-Dev Date: Tue, 30 Jun 2026 23:31:33 -0400 Subject: [PATCH 08/13] map updates --- map.html | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/map.html b/map.html index b7fac34..82a0bb4 100644 --- a/map.html +++ b/map.html @@ -14,9 +14,9 @@
- -
-""" - - -def main(): - rows, features = load_features() - if not features: - print("No homes to map — run pipeline.py and enrich_addresses.py first.") - return - html = build_html(rows, features) - with open(OUT_HTML, "w") as f: - f.write(html) - kind = "all" if INCLUDE_NON_RESIDENTIAL else "residential" - print(f"Wrote {OUT_HTML} with {len(features)} {kind} leads.") - print(f"Open it with: open {OUT_HTML}") +import sys +from solargrader.cli import main if __name__ == "__main__": - main() + main(["map", *sys.argv[1:]]) diff --git a/pipeline.py b/pipeline.py index d52c5fb..cd7757a 100644 --- a/pipeline.py +++ b/pipeline.py @@ -1,935 +1,13 @@ -""" -Solar Grader — Streaming Multi-Tile Pipeline -Processes every LiDAR tile covering a region, one at a time, accumulating scored -homes into DuckDB. Streaming + tile deletion keeps disk bounded; a tiles_done -table makes the run resumable (re-running skips finished tiles). - -Reuses the VALIDATED scoring functions from test_pipeline.py unchanged — only the -orchestration (tile discovery, accumulate, resume, cleanup) is new here. +#!/usr/bin/env python +"""Backward-compatible shim → ``solargrader run``. -Run with (targets accumulate into one DB; every run is resumable): - source .venv/bin/activate - python pipeline.py --region harrisburg # a named area - python pipeline.py --county Dauphin,Cumberland # one or more counties (gridded) - python pipeline.py --region harrisburg-metro # the whole metro in one command - python pipeline.py --bbox -77.0 40.2 -76.8 40.4 # an explicit lon/lat box - python pipeline.py --list # show named regions - python pipeline.py # default REGION_BBOX (back-compat) +The pipeline now lives in the ``solargrader`` package. This keeps the old command +working: python pipeline.py --region harrisburg == solargrader run --region harrisburg """ -import gzip -import json -import math -import os -import re -import time -import warnings -from concurrent.futures import ProcessPoolExecutor, as_completed - -import duckdb -import numpy as np -import pandas as pd -import requests -from pyproj import Transformer -from shapely.geometry import Polygon -from shapely.ops import transform as shapely_transform - -# The scoring algorithm — imported as-is from the validated test pipeline -from test_pipeline import ( - download_tile, - get_tile_srs, - load_tile_points, - get_buildings_from_osm, - get_tmy_data, - clip_lidar_to_building, - extract_roof_planes, - compute_horizon, - calculate_annual_kwh, - grade_home, - grade_potential, - save_result, -) - -warnings.filterwarnings("ignore") - -# ── Configuration ───────────────────────────────────────────────────────────── - -# Current target: Harrisburg, PA + near suburbs (~56 tiles, ~1 hr; resumable). -# Outside Chester County, so BUILDING_SOURCE must be "ms" (national footprints). -REGION_BBOX = (-76.930, 40.240, -76.830, 40.310) # (west, south, east, north) -# Prior target — West Chester: (-75.620, 39.945, -75.595, 39.970) with source "county". - -# TNM serves overlapping LiDAR surveys of wildly different vintage/quality for the -# same area; mixing them corrupts scores. Leave this None to AUTO-SELECT the best -# collection per location (density-ranked greedy mosaic — see select_collections()). -# Set to a name substring (e.g. "PA_17County_D24") to force one collection manually. -PREFERRED_COLLECTION = None - -# Cap tiles for a fast first validation; set to None to process the whole region. -MAX_TILES = None - -# Whole county/state: set to a large bbox to process it as a grid of sub-regions -# (resumable per sub-region AND per tile). None = a single REGION_BBOX run. -WHOLE_REGION_BBOX = None -STEP_DEG = 0.05 # sub-region size (~5 km) when iterating a whole region - -# Building inventory source: "ms" (Microsoft, national — required outside Chester), -# "county" (Chester ArcGIS footprints), or "osm" (Overpass — prototyping only). -BUILDING_SOURCE = "ms" -COUNTY_FOOTPRINTS = ( - "https://mapservices.pasda.psu.edu/server/rest/services/" - "pasda/ChesterCounty/MapServer/14/query" -) - -# Microsoft Building Footprints (national, ODbL) — the scale-out source for areas -# without a county service and for other states. Quadkey-partitioned GeoJSONL: -# download the covering quadkey file once, then filter per tile. See DATA_SOURCES.md. -MS_DATASET_LINKS = "https://minedbuildings.z5.web.core.windows.net/global-buildings/dataset-links.csv" -MS_CACHE_DIR = "ms_cache" -MS_ZOOM = 9 - -# Fallback tile discovery for when TNM's product query is down. Region-agnostic: -# the USGS 3DEP index (below) reports which lidar PROJECT(s) cover a bbox anywhere -# in the US, each with its staged rockyweb LAZ directory (lpc_link). We browse that -# directory and select covering tiles by decoding the UTM grid cell embedded in -# tile filenames (USGS_LPC_..._.laz). Projects that name tiles -# some other way (state-plane, sequential ids) can't be spatially selected by name -# — those are logged and skipped; the primary TNM path (naming-agnostic) still gets -# them. No project name is hardcoded here anymore. -INDEX_LPC = ("https://index.nationalmap.gov/arcgis/rest/services/" - "3DEPElevationIndex/MapServer/8/query") - -DB_PATH = "solar_grader.duckdb" # accumulating production DB (separate from the test) -TILE_CACHE_DIR = "tile_stream" # tiles downloaded here, then deleted after processing -DELETE_TILES_AFTER = True # the storage-capping mechanic; resume covers re-runs - -# Process-level parallelism: tiles are independent, so score N of them at once in -# separate processes. Separate processes (not threads) sidestep the GIL for the -# RANSAC/pvlib CPU work, and each worker deletes its own tile after scoring, so peak -# RAM/disk stays ~N_WORKERS tiles — the streaming/bounded-storage design is intact. -# The single DuckDB writer stays in the main process (DuckDB is single-writer); -# workers return compact results and never touch the DB. Tune down if RAM is tight. -N_WORKERS = min(4, (os.cpu_count() or 2)) - -# Required data attribution — show wherever results are displayed. See DATA_SOURCES.md. -ATTRIBUTION = ("Data: USGS 3DEP · Building footprints © Microsoft (ODbL) · " - "© OpenStreetMap contributors · Chester County GIS/PASDA · EU PVGIS") - - -# ── Tile discovery ──────────────────────────────────────────────────────────── - -def find_all_tiles(bbox, page=50): - """Query USGS TNM for ALL LiDAR tiles covering bbox (paginated, de-duped).""" - west, south, east, north = bbox - url = "https://tnmaccess.nationalmap.gov/api/v1/products" - tiles, offset = {}, 0 - while True: - params = { - "datasets": "Lidar Point Cloud (LPC)", - "bbox": f"{west},{south},{east},{north}", - "prodFormats": "LAZ", - "max": page, - "offset": offset, - } - r = requests.get(url, params=params, timeout=60) - r.raise_for_status() - data = r.json() - items = data.get("items", []) - for it in items: - dl = it.get("downloadURL") - if dl: - tiles[dl] = it - offset += len(items) - if not items or offset >= data.get("total", 0): - break - return list(tiles.values()) - - -def tile_id_of(tile): - """Stable id for a tile = its LAZ filename.""" - return tile["downloadURL"].split("/")[-1].split("?")[0] - - -def _covering_cells(bbox): - """UTM 1 km cells covering bbox, in the bbox's own UTM zone. Returns - {code -> (easting_km, northing_km)} with the zone and its EPSG, where `code` is - the 6-digit 'EEENNN' (easting km, northing km mod 1000) that USGS embeds in tile - filenames. Zone is derived from the bbox, so this works in any UTM zone.""" - w, s, e, n = bbox - zone = int(((w + e) / 2 + 180) // 6) + 1 - epsg = 26900 + zone # NAD83 / UTM zone N - fwd = Transformer.from_crs("EPSG:4326", f"EPSG:{epsg}", always_xy=True) - xs, ys = fwd.transform([w, e, w, e], [s, s, n, n]) - cells = {} - for ee in range(int(min(xs) // 1000), int(max(xs) // 1000) + 1): - for nn in range(int(min(ys) // 1000), int(max(ys) // 1000) + 1): - cells[f"{ee:03d}{nn % 1000:03d}"] = (ee, nn) - return cells, zone, epsg - - -def _lidar_projects_covering(bbox): - """Region-agnostic: the USGS 3DEP index → [(project, staged LAZ dir URL)] for - every lidar project whose extent intersects bbox, anywhere in the US.""" - w, s, e, n = bbox - data = with_retries(lambda: requests.get( - INDEX_LPC, - params={"geometry": f"{w},{s},{e},{n}", "geometryType": "esriGeometryEnvelope", - "inSR": "4326", "spatialRel": "esriSpatialRelIntersects", - "outFields": "project,lpc_link", "returnGeometry": "false", "f": "json"}, - timeout=60).json()) - out = [] - for f in data.get("features", []): - a = f.get("attributes", {}) - link = (a.get("lpc_link") or "").rstrip("/") - if link: - out.append((a.get("project") or link.split("/")[-1], link)) - return out - - -# USGS LPC filename with a UTM/MGRS grid tail: ..._.laz -_TILE_RE = re.compile(r"(USGS_LPC_[\w.\-]+?_(\d{2})[A-Z]{3}(\d{6}))\.laz") - - -def find_tiles_rockyweb(bbox): - """Fallback tile discovery when TNM's product query is down. Discovers covering - projects via the 3DEP index (any US region — no hardcoded project), browses each - project's staged LAZ dir, and selects tiles whose UTM-named grid cell covers - bbox. Projects that name tiles non-spatially (state-plane, sequential ids) can't - be decoded here and are logged + skipped (fetch those via the primary TNM path).""" - cells, zone, epsg = _covering_cells(bbox) - inv = Transformer.from_crs(f"EPSG:{epsg}", "EPSG:4326", always_xy=True) - tiles, seen = [], set() - for project, laz_dir in _lidar_projects_covering(bbox): - laz_url = f"{laz_dir}/LAZ/" - try: - listing = with_retries(lambda: requests.get(laz_url, timeout=90).text) - except Exception as ex: - print(f" {project}: LAZ dir unreachable ({type(ex).__name__}) — skip.") - continue - matched = 0 - for fn_stem, ztile, code in {m for m in _TILE_RE.findall(listing)}: - if int(ztile) != zone or code not in cells or fn_stem in seen: - continue - seen.add(fn_stem) - matched += 1 - e_km, n_km = cells[code] - e0, n0 = e_km * 1000, n_km * 1000 - lons, lats = inv.transform([e0, e0 + 1000], [n0, n0 + 1000]) - tiles.append({ - "downloadURL": f"{laz_url}{fn_stem}.laz", - "title": f"USGS Lidar Point Cloud {project} {fn_stem.split('_')[-1]}", - "sizeInBytes": 0, - "boundingBox": {"minX": min(lons), "maxX": max(lons), - "minY": min(lats), "maxY": max(lats)}, - }) - if matched: - print(f" {project}: {matched} covering tiles (UTM zone {zone}).") - else: - why = ("non-UTM tile naming — skip (use TNM for this one)" - if ".laz" in listing else "no LAZ found") - print(f" {project}: 0 covering tiles ({why}).") - return tiles - - -# ── Building inventory sources ──────────────────────────────────────────────── - -def get_buildings_from_county(bbox): - """Chester County building footprints (PASDA MapServer layer 14), bbox-queried. - Returns the same dict shape as get_buildings_from_osm so scoring is unchanged. - Uses an id-only query then batched geometry fetches (robust to the 1000 cap).""" - west, south, east, north = bbox - geom = { - "geometry": f"{west},{south},{east},{north}", - "geometryType": "esriGeometryEnvelope", - "inSR": "4326", - "spatialRel": "esriSpatialRelIntersects", - } - print(" Querying county building footprints...") - ids = requests.get( - COUNTY_FOOTPRINTS, - params={**geom, "where": "1=1", "returnIdsOnly": "true", "f": "json"}, - timeout=90, - ).json().get("objectIds") or [] - - buildings = [] - for i in range(0, len(ids), 200): - batch = ids[i:i + 200] - resp = requests.post( - COUNTY_FOOTPRINTS, - data={"objectIds": ",".join(map(str, batch)), "outFields": "OBJECTID", - "returnGeometry": "true", "outSR": "4326", "f": "json"}, - timeout=90, - ) - resp.raise_for_status() - for f in resp.json().get("features", []): - rings = (f.get("geometry") or {}).get("rings") - if not rings: - continue - try: - poly = Polygon(rings[0]) # exterior ring - if not poly.is_valid or poly.area == 0: - continue - except Exception: - continue - cx, cy = poly.centroid.x, poly.centroid.y - deg2m2 = (111320 ** 2) * math.cos(math.radians(cy)) - buildings.append({ - "osm_id": int(f["attributes"]["OBJECTID"]), # unique building id - "geometry": poly, "lat": cy, "lon": cx, - "footprint_m2": poly.area * deg2m2, - }) - print(f" Found {len(buildings)} buildings.") - return buildings - - -def _quadkey(lon, lat, z=MS_ZOOM): - """Bing Maps quadkey for a lon/lat at zoom z (MS footprint partition scheme).""" - n = 2 ** z - x = int((lon + 180) / 360 * n) - y = int((1 - math.asinh(math.tan(math.radians(lat))) / math.pi) / 2 * n) - qk = "" - for i in range(z, 0, -1): - digit, mask = 0, 1 << (i - 1) - if x & mask: digit += 1 - if y & mask: digit += 2 - qk += str(digit) - return qk - - -def _ms_link_for(qk): - """Look up the download URL for a US quadkey in MS's dataset index (cached).""" - os.makedirs(MS_CACHE_DIR, exist_ok=True) - links = os.path.join(MS_CACHE_DIR, "dataset-links.csv") - if not os.path.exists(links): - print(" Fetching MS dataset index (one-time)...") - r = requests.get(MS_DATASET_LINKS, timeout=180) - r.raise_for_status() - with open(links, "w") as f: - f.write(r.text) - with open(links) as f: - for line in f: - p = line.split(",") - if len(p) >= 3 and p[0] == "UnitedStates" and p[1] == qk: - return p[2] - return None - - -_MS_CACHE = {} # quadkey -> (centroids ndarray, list-of-rings) - - -def _load_ms_quadkey(qk): - """Download (once) and parse a quadkey's GeoJSONL into centroids + rings.""" - if qk in _MS_CACHE: - return _MS_CACHE[qk] - url = _ms_link_for(qk) - if not url: - _MS_CACHE[qk] = (np.empty((0, 2)), []) - return _MS_CACHE[qk] - path = os.path.join(MS_CACHE_DIR, f"{qk}.csv.gz") - if not os.path.exists(path): - print(f" Downloading MS footprints quadkey {qk} (once)...") - r = requests.get(url, timeout=300) - r.raise_for_status() - with open(path, "wb") as f: - f.write(r.content) - centroids, rings = [], [] - with gzip.open(path, "rt") as f: - for line in f: - try: - coords = json.loads(line)["geometry"]["coordinates"][0] - except Exception: - continue - xs = [c[0] for c in coords] - ys = [c[1] for c in coords] - centroids.append((sum(xs) / len(xs), sum(ys) / len(ys))) - rings.append(coords) - _MS_CACHE[qk] = (np.array(centroids), rings) - print(f" Loaded {len(rings):,} MS buildings for quadkey {qk}.") - return _MS_CACHE[qk] - - -def get_buildings_from_ms(bbox): - """Microsoft Building Footprints within bbox (national/other-states source). - Returns the same dict shape as the other backends.""" - w, s, e, n = bbox - quadkeys = {_quadkey(lo, la) for lo in (w, e) for la in (s, n)} - buildings = [] - for qk in quadkeys: - centroids, rings = _load_ms_quadkey(qk) - if len(centroids) == 0: - continue - cx, cy = centroids[:, 0], centroids[:, 1] - mask = (cx >= w) & (cx <= e) & (cy >= s) & (cy <= n) - qk_int = int(qk) - for idx in np.nonzero(mask)[0]: - try: - poly = Polygon(rings[idx]) - if not poly.is_valid or poly.area == 0: - continue - except Exception: - continue - lon, lat = float(cx[idx]), float(cy[idx]) - deg2m2 = (111320 ** 2) * math.cos(math.radians(lat)) - buildings.append({ - "osm_id": qk_int * 10_000_000 + int(idx), # stable unique id - "geometry": poly, "lat": lat, "lon": lon, - "footprint_m2": poly.area * deg2m2, - }) - print(f" Found {len(buildings)} buildings (MS footprints).") - return buildings - - -def get_buildings(bbox): - """Dispatch to the configured building source.""" - if BUILDING_SOURCE == "county": - return get_buildings_from_county(bbox) - if BUILDING_SOURCE == "ms": - return get_buildings_from_ms(bbox) - return get_buildings_from_osm(bbox) - - -# ── LiDAR collection picker (best-quality mosaic per region) ─────────────────── - -def _collection_of(tile): - """Collection = the staged '/Projects//' segment of the download URL.""" - m = re.search(r"/Projects/([^/]+)/", tile.get("downloadURL", "")) - return m.group(1) if m else tile.get("title", "unknown") - - -def _tile_bbox(tile): - bb = tile.get("boundingBox") or {} - if all(k in bb for k in ("minX", "minY", "maxX", "maxY")): - return (bb["minX"], bb["minY"], bb["maxX"], bb["maxY"]) - return None - - -def _collection_density(ctiles): - """Density proxy = median file bytes per km² (what actually drives RANSAC).""" - vals = [] - for t in ctiles: - bb, sz = _tile_bbox(t), (t.get("sizeInBytes") or 0) - if not bb or sz <= 0: - continue - x0, y0, x1, y1 = bb - lat = (y0 + y1) / 2 - area = abs((x1 - x0) * 111.32 * math.cos(math.radians(lat)) * (y1 - y0) * 110.54) - if area > 0: - vals.append(sz / area) - return float(np.median(vals)) if vals else 0.0 - - -def _rank_key(name, ctiles): - """Rank collections by density, then recency, then QL hint (all descending).""" - years = [int(m.group(1)) for t in ctiles - for m in [re.search(r"(20\d\d)", t.get("publicationDate", ""))] if m] - ql = 2 if "QL1" in name else (1 if "QL2" in name else 0) - return (_collection_density(ctiles), max(years) if years else 0, ql) - - -def select_collections(tiles, bbox, cell_deg=0.01): - """Choose a best-quality LiDAR mosaic: rank collections by density/recency, then - greedily assign each ~1 km cell to the highest-ranked collection covering it, - filling gaps with lower-ranked ones. Returns the tiles needed for that mosaic.""" - colls = {} - for t in tiles: - colls.setdefault(_collection_of(t), []).append(t) - - if len(colls) == 1: - name = next(iter(colls)) - print(f" 1 collection covers this region: {name} ({len(tiles)} tiles).") - return tiles - - ranked = sorted(colls.items(), key=lambda kv: _rank_key(*kv), reverse=True) - print(f" {len(colls)} overlapping collections — ranking by density/recency:") - for name, ct in ranked: - print(f" {name}: ~{_collection_density(ct) / 1e3:.0f} KB/km², {len(ct)} tiles") - - w, s, e, n = bbox - xs, ys = np.arange(w, e, cell_deg), np.arange(s, n, cell_deg) - cells = [(x + cell_deg / 2, y + cell_deg / 2) for x in xs for y in ys] \ - or [((w + e) / 2, (s + n) / 2)] - covered = [False] * len(cells) - chosen, chosen_urls, usage = [], set(), {} - - for name, ct in ranked: - boxes = [(_tile_bbox(t), t) for t in ct] - for ci, (cx, cy) in enumerate(cells): - if covered[ci]: - continue - for bb, t in boxes: - if bb and bb[0] <= cx <= bb[2] and bb[1] <= cy <= bb[3]: - covered[ci] = True - usage[name] = usage.get(name, 0) + 1 - if t["downloadURL"] not in chosen_urls: - chosen_urls.add(t["downloadURL"]) - chosen.append(t) - break - - total = len(cells) - print(" Selected mosaic:") - for name, ct in ranked: - if usage.get(name): - print(f" {name}: {100 * usage[name] / total:.0f}% of area") - dropped = [name for name, ct in ranked if not usage.get(name)] - if dropped: - print(f" dropped (redundant / lower quality): {', '.join(dropped)}") - return chosen - - -def with_retries(fn, *args, attempts=4, base_delay=5, **kwargs): - """Call fn with exponential backoff on transient failures (e.g. Overpass 504, - PVGIS hiccups). Re-raises only after the final attempt.""" - for i in range(attempts): - try: - return fn(*args, **kwargs) - except Exception as e: - if i == attempts - 1: - raise - wait = base_delay * (2 ** i) - print(f" {type(e).__name__}: {e}\n retry {i + 1}/{attempts - 1} in {wait}s...") - time.sleep(wait) - - -# ── Database (accumulating, resumable) ──────────────────────────────────────── - -def setup_db(path): - con = duckdb.connect(path) - # Same homes schema as the test pipeline, but NEVER dropped — we accumulate. - con.execute(""" - CREATE TABLE IF NOT EXISTS homes ( - osm_id BIGINT PRIMARY KEY, - lat DOUBLE, - lon DOUBLE, - footprint_m2 DOUBLE, - usable_roof_area_m2 DOUBLE, - primary_tilt_deg DOUBLE, - primary_azimuth_deg DOUBLE, - roof_plane_count INTEGER, - res_annual_kwh DOUBLE, - res_system_kw DOUBLE, - res_panel_count INTEGER, - solar_score INTEGER, - solar_grade VARCHAR(2), - max_annual_kwh DOUBLE, - max_system_kw DOUBLE, - max_panel_count INTEGER, - potential_grade VARCHAR(2), - shade_loss_pct DOUBLE, - processed_at TIMESTAMP - ) - """) - # Resume ledger: which tiles are already done. - con.execute(""" - CREATE TABLE IF NOT EXISTS tiles_done ( - tile_id VARCHAR PRIMARY KEY, - n_scored INTEGER, - processed_at TIMESTAMP - ) - """) - return con - - -# ── Per-building scoring (mirrors test_pipeline.main's inner loop) ───────────── - -def score_building(building, tile_points, transformer, tmy): - """Full scoring chain for one building. Returns a save_result payload or None.""" - points = clip_lidar_to_building(tile_points, building["geometry"], transformer) - if points is None: - return None - planes = extract_roof_planes(points) - if not planes: - return None - - observer = points.mean(axis=0) - own_footprint = shapely_transform(transformer.transform, building["geometry"]).buffer(1.0) - horizon = compute_horizon(tile_points, observer, own_footprint) - - sysres = calculate_annual_kwh(building["lat"], building["lon"], tmy, planes, horizon) - if sysres is None or sysres["res_panels"] == 0: - return None - - best_plane = sysres["best_plane"] - total_area = sum(p["area_m2"] for p in planes) - grade, score = grade_home( - sysres["res_kwh"], sysres["res_kw"], best_plane["azimuth_deg"] - ) - return { - "usable_area_m2": total_area, - "primary_tilt": best_plane["tilt_deg"], - "primary_azimuth": best_plane["azimuth_deg"], - "plane_count": len(planes), - "res_kwh": sysres["res_kwh"], - "res_kw": sysres["res_kw"], - "res_panels": sysres["res_panels"], - "max_kwh": sysres["max_kwh"], - "max_kw": sysres["max_kw"], - "max_panels": sysres["max_panels"], - "potential_grade": grade_potential(sysres["max_kw"]), - "shade_loss_pct": sysres["shade_loss_pct"], - "grade": grade, - "score": score, - } - - -def tile_latlon_bbox(tile_points, tile_srs): - """Derive the tile's lon/lat bbox from its actual point extent (UTM → 4326).""" - _, X, Y = tile_points - inv = Transformer.from_crs(tile_srs, "EPSG:4326", always_xy=True) - xs = [X.min(), X.max(), X.min(), X.max()] - ys = [Y.min(), Y.min(), Y.max(), Y.max()] - lons, lats = inv.transform(xs, ys) - return (min(lons), min(lats), max(lons), max(lats)) - - -# ── Parallel worker ─────────────────────────────────────────────────────────── - -def _process_tile(tile, tile_cache=TILE_CACHE_DIR): - """Fully process ONE tile in a worker process: download, discover buildings, - score each, then delete the tile file. Returns a compact, picklable result — - no shapely geometry crosses the process boundary (save_result only needs - osm_id/lat/lon/footprint). Never touches the DB; the parent does all writes. - All failures are caught and returned so one bad tile can't kill the pool.""" - tid = tile_id_of(tile) - try: - path = download_tile(tile["downloadURL"], tile_cache) - tile_points = load_tile_points(path) - tile_srs = get_tile_srs(path) - transformer = Transformer.from_crs("EPSG:4326", tile_srs, always_xy=True) - - tbbox = tile_latlon_bbox(tile_points, tile_srs) - buildings = with_retries(get_buildings, tbbox) - if not buildings: - _cleanup(path) - return {"tid": tid, "n_buildings": 0, "scored": [], "error": None} - - # One irradiance fetch per tile (its center) covers all its buildings - cx = (tbbox[0] + tbbox[2]) / 2 - cy = (tbbox[1] + tbbox[3]) / 2 - tmy = with_retries(get_tmy_data, cy, cx) - - scored = [] - for building in buildings: - result = score_building(building, tile_points, transformer, tmy) - if result: - meta = {k: building[k] for k in ("osm_id", "lat", "lon", "footprint_m2")} - scored.append((meta, result)) - - _cleanup(path) - return {"tid": tid, "n_buildings": len(buildings), "scored": scored, "error": None} - - except Exception as e: - # Return the failure instead of raising — keeps the other workers running. - return {"tid": tid, "n_buildings": 0, "scored": [], - "error": f"{type(e).__name__}: {e}"} - - -# ── Main streaming loop ─────────────────────────────────────────────────────── - -def run(bbox, db_path=DB_PATH, tile_cache=TILE_CACHE_DIR): - print("=" * 60) - print("SOLAR GRADER — Streaming Multi-Tile Pipeline") - print(f"Region bbox: {bbox}") - print("=" * 60) - - print("\nDiscovering LiDAR tiles...") - all_tiles = with_retries(find_all_tiles, bbox) - if not all_tiles: - print(" TNM returned 0 (spatial query down) — falling back to direct") - print(" rockyweb tile discovery...") - all_tiles = with_retries(find_tiles_rockyweb, bbox) - print(f" rockyweb found {len(all_tiles)} tiles.") - if not all_tiles: - print(" ABORT: no tiles from TNM or rockyweb. Both may be down — retry later.") - return - print(f" {len(all_tiles)} tiles cover this region.") - if PREFERRED_COLLECTION: - tiles = [t for t in all_tiles if PREFERRED_COLLECTION in _collection_of(t)] - print(f" Forced collection '{PREFERRED_COLLECTION}': {len(tiles)} tiles.") - else: - tiles = select_collections(all_tiles, bbox) - if not tiles: - print(" ABORT: no tiles selected for this region.") - return - if MAX_TILES is not None: - tiles = tiles[:MAX_TILES] - print(f" MAX_TILES set — processing only the first {len(tiles)} this run.") - - con = setup_db(db_path) - done = {row[0] for row in con.execute("SELECT tile_id FROM tiles_done").fetchall()} - pending = [t for t in tiles if tile_id_of(t) not in done] - if done: - print(f" {len(done)} already processed — resuming, will skip those.") - print(f"\nScoring {len(pending)} tiles across {N_WORKERS} worker processes " - f"(single DB writer in this process)...") - - # Dispatch each tile to a worker; the parent stays the sole DB writer, saving - # results and marking tiles done as they stream back (order is non-deterministic). - t_start = time.time() - completed = 0 - with ProcessPoolExecutor(max_workers=N_WORKERS) as pool: - futures = {pool.submit(_process_tile, t, tile_cache): tile_id_of(t) - for t in pending} - for fut in as_completed(futures): - completed += 1 - res = fut.result() # worker catches its own errors; never raises here - tid = res["tid"] - if res["error"]: - print(f"[{completed}/{len(pending)}] {tid} — ERROR: {res['error']} (skipped)") - continue - for meta, result in res["scored"]: - save_result(con, meta, result) - _mark_done(con, tid, len(res["scored"])) - print(f"[{completed}/{len(pending)}] {tid} — " - f"scored {len(res['scored'])}/{res['n_buildings']}") - - _summary(con, time.time() - t_start) - con.close() - - -def _mark_done(con, tid, n_scored): - con.execute( - "INSERT OR REPLACE INTO tiles_done VALUES (?, ?, ?)", - [tid, n_scored, pd.Timestamp.now()], - ) - - -def _cleanup(path): - if DELETE_TILES_AFTER and os.path.exists(path): - os.remove(path) - - -def _summary(con, elapsed): - print(f"\n\n{'=' * 60}") - print("PIPELINE SUMMARY") - print(f"{'=' * 60}") - n_homes = con.execute("SELECT COUNT(*) FROM homes").fetchone()[0] - n_tiles = con.execute("SELECT COUNT(*) FROM tiles_done").fetchone()[0] - print(f" Tiles processed: {n_tiles}") - print(f" Homes scored: {n_homes}") - print(f" Total time: {elapsed:.1f}s") - - print(f"\n Grade breakdown:") - rows = con.execute(""" - SELECT solar_grade, COUNT(*) FROM homes GROUP BY solar_grade - """).fetchall() - counts = dict(rows) - for g in ["A+", "A", "B+", "B", "C", "D"]: - n = counts.get(g, 0) - if n: - print(f" {g:2s} {'█' * min(n, 50)} ({n})") - - print(f"\n Results saved to: {DB_PATH}") - - -def run_region(bbox, step_deg=STEP_DEG, db_path=DB_PATH, tile_cache=TILE_CACHE_DIR): - """Process a whole county/state by tiling it into sub-region bboxes and running - each through the normal pipeline. Resumable: a regions_done ledger skips finished - sub-regions, and tiles_done (inside run) dedups tiles shared across sub-regions.""" - w, s, e, n = bbox - nx = max(1, math.ceil((e - w) / step_deg)) # integer step counts avoid float drift - ny = max(1, math.ceil((n - s) / step_deg)) - subs = [] - for ix in range(nx): - x0, x1 = w + ix * step_deg, min(w + (ix + 1) * step_deg, e) - for iy in range(ny): - y0, y1 = s + iy * step_deg, min(s + (iy + 1) * step_deg, n) - if x1 - x0 > 1e-9 and y1 - y0 > 1e-9: # skip degenerate slivers - subs.append((round(x0, 4), round(y0, 4), round(x1, 4), round(y1, 4))) - - print("#" * 60) - print(f"REGION RUN — {len(subs)} sub-regions (~{step_deg}°) over {bbox}") - print("#" * 60) - - # regions_done ledger. Opened only BETWEEN run() calls so there's never a - # second connection open while run() holds the DuckDB writer. - con = duckdb.connect(db_path) - con.execute("CREATE TABLE IF NOT EXISTS regions_done " - "(region_key VARCHAR PRIMARY KEY, processed_at TIMESTAMP)") - done = {r[0] for r in con.execute("SELECT region_key FROM regions_done").fetchall()} - con.close() - if done: - print(f" Resuming — {len(done)} sub-regions already complete.") - - for i, sub in enumerate(subs, 1): - key = ",".join(f"{v:.4f}" for v in sub) - if key in done: - print(f"\n### sub-region {i}/{len(subs)} {sub} — already done, skip.") - continue - print(f"\n### sub-region {i}/{len(subs)} {sub}") - try: - run(sub, db_path, tile_cache) - except Exception as ex: - print(f"### sub-region {i} ERROR: {type(ex).__name__}: {ex} — continuing.") - continue - con = duckdb.connect(db_path) - con.execute("INSERT OR REPLACE INTO regions_done VALUES (?, ?)", - [key, pd.Timestamp.now()]) - con.close() - - con = duckdb.connect(db_path, read_only=True) - total = con.execute("SELECT COUNT(*) FROM homes").fetchone()[0] - con.close() - print("\n" + "#" * 60) - print(f"REGION RUN COMPLETE — {total} homes scored across {len(subs)} sub-regions.") - print("#" * 60) - - -# ── Area targeting: name/county → bbox (so new runs need no source edits) ───── - -# US Census TIGERweb county boundaries (free, US Gov). Resolves any county name to -# an exact lon/lat bbox — no hand-typed coordinates, scales to every US county. -TIGERWEB_COUNTIES = ( - "https://tigerweb.geo.census.gov/arcgis/rest/services/" - "TIGERweb/State_County/MapServer/1/query" -) -_STATE_FIPS = { - "AL": "01", "AK": "02", "AZ": "04", "AR": "05", "CA": "06", "CO": "08", - "CT": "09", "DE": "10", "DC": "11", "FL": "12", "GA": "13", "HI": "15", - "ID": "16", "IL": "17", "IN": "18", "IA": "19", "KS": "20", "KY": "21", - "LA": "22", "ME": "23", "MD": "24", "MA": "25", "MI": "26", "MN": "27", - "MS": "28", "MO": "29", "MT": "30", "NE": "31", "NV": "32", "NH": "33", - "NJ": "34", "NM": "35", "NY": "36", "NC": "37", "ND": "38", "OH": "39", - "OK": "40", "OR": "41", "PA": "42", "RI": "44", "SC": "45", "SD": "46", - "TN": "47", "TX": "48", "UT": "49", "VT": "50", "VA": "51", "WA": "53", - "WV": "54", "WI": "55", "WY": "56", -} - -# Convenience aliases for the Harrisburg region — type `--region `. Each is -# either an explicit bbox or one/more counties resolved live via TIGERweb. -REGIONS = { - "harrisburg": {"bbox": (-76.930, 40.240, -76.830, 40.310)}, # city + near suburbs - "dauphin": {"counties": ["Dauphin"]}, - "cumberland": {"counties": ["Cumberland"]}, - "perry": {"counties": ["Perry"]}, - "york": {"counties": ["York"]}, - "lancaster": {"counties": ["Lancaster"]}, - "lebanon": {"counties": ["Lebanon"]}, - # The whole Harrisburg metro in one command (each county is a resumable grid). - "harrisburg-metro": {"counties": ["Dauphin", "Cumberland", "Perry"]}, -} - - -def county_bbox(county, state="PA"): - """Exact lon/lat bbox for a US county via TIGERweb. `county` is the bare name - (e.g. 'Dauphin', not 'Dauphin County'). Validated against a strict whitelist - before it reaches the REST where-clause (trust-boundary hygiene).""" - if not re.fullmatch(r"[A-Za-z][A-Za-z .'\-]{0,48}", county or ""): - raise ValueError(f"Invalid county name: {county!r}") - fips = _STATE_FIPS.get(state.upper()) - if not fips: - raise ValueError(f"Unknown state: {state!r}") - r = with_retries(lambda: requests.get( - TIGERWEB_COUNTIES, - params={"where": f"BASENAME='{county}' AND STATE='{fips}'", - "outFields": "NAME", "returnGeometry": "true", - "outSR": "4326", "f": "json"}, - timeout=60, - ).json()) - feats = r.get("features") or [] - if not feats: - raise ValueError(f"No county '{county}' found in {state.upper()}.") - rings = feats[0]["geometry"]["rings"] - xs = [p[0] for ring in rings for p in ring] - ys = [p[1] for ring in rings for p in ring] - return (min(xs), min(ys), max(xs), max(ys)), feats[0]["attributes"]["NAME"] - - -def all_counties(state="PA"): - """Every county name in a state, alphabetized, via TIGERweb. Powers whole-state - runs — each county then runs as its own resumable gridded region.""" - fips = _STATE_FIPS.get(state.upper()) - if not fips: - raise ValueError(f"Unknown state: {state!r}") - data = with_retries(lambda: requests.get( - TIGERWEB_COUNTIES, - params={"where": f"STATE='{fips}'", "outFields": "BASENAME", - "returnGeometry": "false", "orderByFields": "BASENAME", "f": "json"}, - timeout=90).json()) - return [f["attributes"]["BASENAME"] for f in data.get("features", [])] - - -def _targets_from_args(args): - """Resolve CLI args into a list of (label, bbox, use_grid) run targets. All - targets accumulate into the same DB; the ledgers dedup any overlap.""" - if getattr(args, "all_counties", False): - names = all_counties(args.state) - print(f"Resolving {len(names)} counties in {args.state.upper()}...") - targets = [] - for name in names: - bbox, label = county_bbox(name, args.state) - targets.append((label, bbox, True)) # each county is a gridded run - return targets - if args.bbox: - return [("bbox", tuple(args.bbox), args.grid)] - if args.county: - names = [c.strip() for c in args.county.split(",") if c.strip()] - out = [] - for name in names: - bbox, label = county_bbox(name, args.state) - out.append((label, bbox, True)) # a county is always a gridded run - return out - if args.region: - spec = REGIONS.get(args.region.lower()) - if not spec: - raise SystemExit(f"Unknown region '{args.region}'. " - f"Options: {', '.join(sorted(REGIONS))}") - if "bbox" in spec: - return [(args.region, spec["bbox"], args.grid)] - out = [] - for name in spec["counties"]: - bbox, label = county_bbox(name, spec.get("state", "PA")) - out.append((label, bbox, True)) - return out - # No target given — preserve the original default behaviour. - if WHOLE_REGION_BBOX: - return [("default-region", WHOLE_REGION_BBOX, True)] - return [("default", REGION_BBOX, False)] - - -def main_cli(): - import argparse - p = argparse.ArgumentParser( - description="Solar Grader — score a region's rooftops. Targets accumulate " - "into one DB; every run is resumable (Ctrl-C and re-run).") - g = p.add_mutually_exclusive_group() - g.add_argument("--county", metavar="A,B,C", - help="one or more county names (comma-separated), e.g. Dauphin,Cumberland") - g.add_argument("--region", metavar="NAME", - help=f"a named area: {', '.join(sorted(REGIONS))}") - g.add_argument("--bbox", nargs=4, type=float, metavar=("W", "S", "E", "N"), - help="explicit lon/lat bounding box") - g.add_argument("--all-counties", action="store_true", - help="every county in --state — a whole-state run (long, but " - "fully resumable county-by-county and tile-by-tile)") - p.add_argument("--state", default="PA", help="state for --county / --all-counties (default PA)") - p.add_argument("--grid", action="store_true", - help="force gridded run_region for a --bbox/--region (auto-on for counties)") - p.add_argument("--step", type=float, default=STEP_DEG, - help=f"sub-region size in degrees when gridding (default {STEP_DEG})") - p.add_argument("--db", default=DB_PATH, help=f"DuckDB path (default {DB_PATH})") - p.add_argument("--list", action="store_true", help="list named regions and exit") - args = p.parse_args() - - if args.list: - print("Named regions (--region):") - for name, spec in sorted(REGIONS.items()): - what = spec.get("bbox") or "counties: " + ", ".join(spec["counties"]) - print(f" {name:18s} {what}") - return - - targets = _targets_from_args(args) - print(f"Targets to process ({len(targets)}), all into {args.db}:") - for label, bbox, grid in targets: - print(f" • {label:16s} {tuple(round(v, 3) for v in bbox)} " - f"{'[grid]' if grid else '[single]'}") - - for label, bbox, grid in targets: - print(f"\n{'=' * 60}\n▶ {label}\n{'=' * 60}") - if grid: - run_region(bbox, step_deg=args.step, db_path=args.db) - else: - run(bbox, db_path=args.db) +import sys +from solargrader.cli import main if __name__ == "__main__": - main_cli() + main(["run", *sys.argv[1:]]) diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..10a4d5a --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,61 @@ +[build-system] +requires = ["setuptools>=61"] +build-backend = "setuptools.build_meta" + +[project] +name = "solargrader" +version = "0.1.0" +description = "Score every rooftop in a region for solar sales potential from free public data" +readme = "README.md" +requires-python = ">=3.9" +license = { text = "MIT" } +authors = [{ name = "SolarGrader contributors" }] +keywords = ["solar", "lidar", "gis", "lead-scoring", "pvlib", "usgs-3dep"] +classifiers = [ + "Programming Language :: Python :: 3", + "License :: OSI Approved :: MIT License", + "Topic :: Scientific/Engineering :: GIS", +] + +# Core runtime deps that pip installs cleanly. The LiDAR reader (PDAL) is a system +# library — install it separately (see the "lidar" extra and the README). +dependencies = [ + "numpy", + "scipy", + "pandas", + "requests", + "shapely", + "pyproj", + "duckdb", + "pvlib", + "open3d", +] + +[project.optional-dependencies] +# PDAL's Python bindings need a system libpdal (best via conda-forge). Kept optional +# so `pip install solargrader` works for everything except reading .laz tiles. +lidar = ["pdal"] +dev = ["pytest", "ruff"] + +[project.scripts] +solargrader = "solargrader.cli:main" + +[project.urls] +Homepage = "https://github.com/" +Documentation = "https://github.com/" + +[tool.setuptools.packages.find] +include = ["solargrader*"] + +[tool.pytest.ini_options] +testpaths = ["tests"] + +[tool.ruff] +line-length = 100 +target-version = "py39" + +[tool.ruff.lint] +# Sensible defaults; the pipeline deliberately catches broad exceptions at trust +# boundaries (network, untrusted geometry), so BLE001 is silenced inline there. +select = ["E", "F", "W", "I", "UP", "B"] +ignore = ["E203"] diff --git a/regrade.py b/regrade.py index 2c06c29..44c60d7 100644 --- a/regrade.py +++ b/regrade.py @@ -1,66 +1,9 @@ -""" -Solar Grader — Re-grade in place (no LiDAR re-run) +#!/usr/bin/env python +"""Backward-compatible shim → ``solargrader regrade``.""" -Re-applies the current grade_home() / grade_potential() logic to every home already -scored in the DB, using the stored physical measurements (yield, kW, azimuth). Grading -is a cheap pure function of columns we already persisted, so tuning the grade never -requires reprocessing point clouds — run this instead. - -Run with: - source .venv/bin/activate - python regrade.py -""" - -import warnings - -import duckdb -import pandas as pd - -from test_pipeline import grade_home, grade_potential - -warnings.filterwarnings("ignore") - -DB_PATH = "solar_grader.duckdb" - - -def main(): - con = duckdb.connect(DB_PATH) - homes = con.execute( - "SELECT osm_id, res_annual_kwh, res_system_kw, primary_azimuth_deg, " - "max_system_kw, solar_grade FROM homes" - ).fetchdf() - if homes.empty: - print("No homes to re-grade — run pipeline.py first.") - return - - before = homes["solar_grade"].value_counts().to_dict() - - rows = [] - for r in homes.itertuples(index=False): - grade, score = grade_home(r.res_annual_kwh, r.res_system_kw, r.primary_azimuth_deg) - rows.append((int(r.osm_id), grade, int(score), grade_potential(r.max_system_kw))) - - new_df = pd.DataFrame(rows, columns=["osm_id", "grade", "score", "pot"]) - con.register("new_df", new_df) - con.execute(""" - UPDATE homes SET - solar_grade = n.grade, solar_score = n.score, potential_grade = n.pot - FROM new_df n WHERE homes.osm_id = n.osm_id - """) - - total = len(homes) - order = ["A+", "A", "B+", "B", "C", "D"] - after = con.execute( - "SELECT solar_grade, COUNT(*) FROM homes GROUP BY 1" - ).fetchdf().set_index("solar_grade")["count_star()"].to_dict() - - print(f"Re-graded {total} homes.\n") - print(f" {'grade':>5} {'before':>8} {'after':>8}") - for g in order: - b, a = before.get(g, 0), after.get(g, 0) - print(f" {g:>5} {b:>8} {a:>8} ({100*a/total:>4.1f}%)") - con.close() +import sys +from solargrader.cli import main if __name__ == "__main__": - main() + main(["regrade", *sys.argv[1:]]) diff --git a/requirements.txt b/requirements.txt index a7b6f69..6e50196 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,9 +1,18 @@ -pdal -pvlib -open3d -duckdb -geopandas +# Runtime dependencies. The package is the source of truth — prefer: +# pip install -e ".[dev]" +# This file mirrors the pip-installable deps for a plain `pip install -r`. +# +# The LiDAR reader (PDAL) is a SYSTEM library; install it via conda-forge: +# conda install -c conda-forge pdal python-pdal +# then it is picked up here. Everything else installs from PyPI. + numpy scipy -requests pandas +requests +shapely +pyproj +duckdb +pvlib +open3d +pdal # requires system libpdal (see note above) diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..02376b8 --- /dev/null +++ b/setup.py @@ -0,0 +1,7 @@ +"""Compatibility shim so `pip install -e .` works on older pip (< 21.3), which +can't do editable installs from pyproject.toml alone. All real metadata lives in +pyproject.toml — this just delegates to setuptools.""" + +from setuptools import setup + +setup() diff --git a/solargrader/__init__.py b/solargrader/__init__.py new file mode 100644 index 0000000..4dd850e --- /dev/null +++ b/solargrader/__init__.py @@ -0,0 +1,27 @@ +""" +SolarGrader — score every rooftop in a region for solar sales potential from free +public data (USGS 3DEP LiDAR, Microsoft/OSM/county building footprints, EU PVGIS +irradiance), then surface the best leads. + +The package is organized so the numerical core is pure functions (easy to test) and +only genuine state/interfaces are classes: + + config Config dataclass — every tunable in one place + models lightweight dataclasses: Building, RoofPlane, SystemResult, ScoredHome + geometry pure: CRS reprojection, LiDAR clipping, area conversion + roof pure: RANSAC roof-plane extraction + solar pure: panel packing, plane-of-array irradiance, annual energy, shading + grading pure: yield-based residential grade + potential tier + sources/ data acquisition: tiles, buildings, irradiance, regions + storage ResultStore — owns the DuckDB connection (schema, writes, resume ledgers) + pipeline streaming, parallel, resumable orchestration + cli command-line entry points + +Heavy optional dependencies (PDAL, Open3D, pvlib) are imported lazily inside the +functions that need them, so the pure modules import with only numpy/shapely/pyproj. +""" + +from .config import Config + +__version__ = "0.1.0" +__all__ = ["Config", "__version__"] diff --git a/solargrader/cli.py b/solargrader/cli.py new file mode 100644 index 0000000..9a187d8 --- /dev/null +++ b/solargrader/cli.py @@ -0,0 +1,142 @@ +""" +Command-line interface. + + solargrader run --region harrisburg # a named area + solargrader run --county Dauphin,Cumberland # one or more counties (gridded) + solargrader run --all-counties --state PA # a whole state (long, resumable) + solargrader run --bbox -77.0 40.2 -76.8 40.4 # an explicit lon/lat box + solargrader enrich # attach addresses + residential flag + solargrader map # write map.html + solargrader regrade # re-grade in place, no reprocessing + +Targets accumulate into one DB and every run is resumable (Ctrl-C and re-run). +""" + +from __future__ import annotations + +import argparse + +from .config import Config +from .sources.regions import REGIONS, resolve_targets + + +def _config_from_args(args) -> Config: + """Build a Config, overriding only the fields the user set.""" + overrides = {} + if getattr(args, "db", None): + overrides["db_path"] = args.db + if getattr(args, "source", None): + overrides["building_source"] = args.source + if getattr(args, "workers", None): + overrides["n_workers"] = args.workers + if getattr(args, "max_tiles", None) is not None: + overrides["max_tiles"] = args.max_tiles + if getattr(args, "preferred_collection", None): + overrides["preferred_collection"] = args.preferred_collection + if getattr(args, "step", None): + overrides["step_deg"] = args.step + return Config(**overrides) + + +def _cmd_run(args) -> None: + from .pipeline import run, run_region + + cfg = _config_from_args(args) + if args.list: + print("Named regions (--region):") + for name, spec in sorted(REGIONS.items()): + what = spec.get("bbox") or "counties: " + ", ".join(spec["counties"]) + print(f" {name:18s} {what}") + return + + # Default to the Harrisburg region when no target is given (backward compatible). + have_target = any([args.bbox, args.county, args.region, args.all_counties]) + targets = resolve_targets( + cfg, bbox=args.bbox, county=args.county, + region=args.region if have_target else "harrisburg", + all_counties_flag=args.all_counties, state=args.state, grid=args.grid, + ) + + print(f"Targets to process ({len(targets)}), all into {cfg.db_path}:") + for label, bbox, grid in targets: + print(f" • {label:16s} {tuple(round(v, 3) for v in bbox)} " + f"{'[grid]' if grid else '[single]'}") + for label, bbox, grid in targets: + print(f"\n{'=' * 60}\n▶ {label}\n{'=' * 60}") + (run_region if grid else run)(bbox, cfg) + + +def _cmd_enrich(args) -> None: + from .enrich import enrich + enrich(_config_from_args(args)) + + +def _cmd_map(args) -> None: + from .webmap import build_map + build_map(_config_from_args(args), out_html=args.out, + include_non_residential=args.all) + + +def _cmd_regrade(args) -> None: + from .regrade import regrade + regrade(_config_from_args(args)) + + +def _add_db_arg(p) -> None: + p.add_argument("--db", default=None, help="DuckDB path (default from Config)") + + +def build_parser() -> argparse.ArgumentParser: + p = argparse.ArgumentParser(prog="solargrader", description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + sub = p.add_subparsers(dest="command", required=True) + + # run + r = sub.add_parser("run", help="score rooftops in a region") + g = r.add_mutually_exclusive_group() + g.add_argument("--county", metavar="A,B,C", help="county names (comma-separated)") + g.add_argument("--region", metavar="NAME", help=f"named area: {', '.join(sorted(REGIONS))}") + g.add_argument("--bbox", nargs=4, type=float, metavar=("W", "S", "E", "N")) + g.add_argument("--all-counties", action="store_true", + help="every county in --state (whole-state, resumable)") + r.add_argument("--state", default="PA", help="state for --county/--all-counties (default PA)") + r.add_argument("--grid", action="store_true", help="force gridded run for a --bbox/--region") + r.add_argument("--step", type=float, default=None, help="sub-region degrees when gridding") + r.add_argument("--source", choices=["ms", "county", "osm"], default=None, + help="building footprint source (default ms)") + r.add_argument("--workers", type=int, default=None, help="parallel worker processes") + r.add_argument("--max-tiles", type=int, default=None, dest="max_tiles", + help="cap tiles processed this run") + r.add_argument("--preferred-collection", default=None, dest="preferred_collection", + help="force one LiDAR collection (name substring)") + r.add_argument("--list", action="store_true", help="list named regions and exit") + _add_db_arg(r) + r.set_defaults(func=_cmd_run) + + # enrich + e = sub.add_parser("enrich", help="attach addresses + residential flag") + _add_db_arg(e) + e.set_defaults(func=_cmd_enrich) + + # map + m = sub.add_parser("map", help="write a Leaflet map.html") + _add_db_arg(m) + m.add_argument("--out", default="map.html", help="output HTML path") + m.add_argument("--all", action="store_true", help="include non-residential leads") + m.set_defaults(func=_cmd_map) + + # regrade + rg = sub.add_parser("regrade", help="re-grade in place (no reprocessing)") + _add_db_arg(rg) + rg.set_defaults(func=_cmd_regrade) + + return p + + +def main(argv: list[str] | None = None) -> None: + args = build_parser().parse_args(argv) + args.func(args) + + +if __name__ == "__main__": + main() diff --git a/solargrader/config.py b/solargrader/config.py new file mode 100644 index 0000000..625b33c --- /dev/null +++ b/solargrader/config.py @@ -0,0 +1,105 @@ +""" +Central configuration for the whole pipeline. + +Everything tunable lives on one ``Config`` dataclass instead of scattered module +globals, so a run is fully described by a single object (easy to override per run, +easy to test, easy to serialize into worker processes). Secrets are never stored +here as literals — they are read from the environment on demand (see +``require_nrel_api_key``). +""" + +from __future__ import annotations + +import os +from dataclasses import dataclass, field + + +@dataclass +class Config: + """All pipeline settings. Construct with defaults and override fields as needed: + + >>> cfg = Config(building_source="ms", n_workers=8) + """ + + # ── Panel / system model (all physical — change to model different hardware) ── + panel_width_m: float = 1.05 # standard residential module footprint + panel_height_m: float = 1.74 + panel_watts: int = 400 # module rated DC power at STC (1000 W/m²) + roof_setback_m: float = 0.46 # fire-code clear pathway from plane edges (~18") + panel_packing: float = 0.90 # fraction of setback area that packs into modules + system_losses: float = 0.86 # inverter + wiring + soiling (~14% loss) + max_azimuth_offset: float = 120.0 # skip planes facing >this many ° from due south + max_system_kw: float = 10.0 # cap for a typical sellable residential install + shade_radius_m: float = 50.0 # search radius for shading obstructions + horizon_bins: int = 72 # azimuth bins for the roof skyline (5° each) + + # ── RANSAC roof-plane fitting ──────────────────────────────────────────────── + min_roof_points: int = 20 # skip buildings with fewer LiDAR points + ransac_iterations: int = 300 # lower = faster, higher = more accurate + ransac_distance_threshold: float = 0.20 # meters a point may sit off a plane + + # ── Pipeline / storage ─────────────────────────────────────────────────────── + db_path: str = "solar_grader.duckdb" # accumulating, resumable results DB + tile_cache_dir: str = "tile_stream" # tiles land here, deleted after scoring + delete_tiles_after: bool = True # storage cap; resume covers re-runs + n_workers: int = field(default_factory=lambda: min(4, os.cpu_count() or 2)) + step_deg: float = 0.05 # sub-region size when gridding a region + + # ── Building inventory source: "ms" | "county" | "osm" ─────────────────────── + building_source: str = "ms" + + # ── LiDAR collection selection ─────────────────────────────────────────────── + # None → auto-select the best-quality mosaic per region (density-ranked). Set to + # a project-name substring to force one collection. ``max_tiles`` caps a run for + # a quick first validation (None = process the whole region). + preferred_collection: str | None = None + max_tiles: int | None = None + + # ── Service endpoints (free data sources) ──────────────────────────────────── + tnm_url: str = "https://tnmaccess.nationalmap.gov/api/v1/products" + index_lpc_url: str = ("https://index.nationalmap.gov/arcgis/rest/services/" + "3DEPElevationIndex/MapServer/8/query") + county_footprints_url: str = ("https://mapservices.pasda.psu.edu/server/rest/services/" + "pasda/ChesterCounty/MapServer/14/query") + ms_dataset_links_url: str = ("https://minedbuildings.z5.web.core.windows.net/" + "global-buildings/dataset-links.csv") + ms_cache_dir: str = "ms_cache" + ms_zoom: int = 9 + overpass_url: str = "https://overpass-api.de/api/interpreter" + tigerweb_counties_url: str = ("https://tigerweb.geo.census.gov/arcgis/rest/services/" + "TIGERweb/State_County/MapServer/1/query") + + # ── Networking ─────────────────────────────────────────────────────────────── + http_user_agent: str = "SolarGrader/1.0 (https://github.com/)" + retry_attempts: int = 4 + retry_base_delay: float = 5.0 + + # ── Required attribution — display wherever results are shown (DATA_SOURCES.md) + attribution: str = ("Data: USGS 3DEP · Building footprints © Microsoft (ODbL) · " + "© OpenStreetMap contributors · Chester County GIS/PASDA · EU PVGIS") + + # ── Derived quantities (kept as properties so they track the fields above) ──── + @property + def panel_area_m2(self) -> float: + """Module area in m² (width × height).""" + return self.panel_width_m * self.panel_height_m + + @property + def module_efficiency(self) -> float: + """STC efficiency implied by the rated watts and module area.""" + return self.panel_watts / (self.panel_area_m2 * 1000.0) + + @property + def residential_panel_cap(self) -> int: + """Whole panels that fit under the residential kW cap.""" + return int(self.max_system_kw * 1000 / self.panel_watts) + + # ── Secrets: read from the environment, never defaulted (a defaulted secret is + # a leaked secret). Only needed if you wire the optional NREL NSRDB source. ─ + def require_nrel_api_key(self) -> str: + """Return NREL_API_KEY from the environment, raising if it is unset.""" + return os.environ["NREL_API_KEY"] + + +# A ready-to-use default instance for simple scripts and the CLI's defaults. +DEFAULT = Config() diff --git a/solargrader/enrich.py b/solargrader/enrich.py new file mode 100644 index 0000000..47d2055 --- /dev/null +++ b/solargrader/enrich.py @@ -0,0 +1,135 @@ +""" +Address enrichment — attach real street addresses and a residential flag to scored +homes by matching each building centroid to the nearest official address point. + +Decoupled from scoring: it operates on the existing DB, is re-runnable, and needs no +LiDAR reprocessing. The address-point service is county-specific (Chester County / +PASDA by default); where no service covers a region, homes are still classified +residential by footprint size, just without street addresses. +""" + +from __future__ import annotations + +import duckdb +import numpy as np +import pandas as pd +import requests +from scipy.spatial import cKDTree + +from .config import Config + +# Chester County address points (PASDA). Map the right service per county to extend. +ADDR_LAYER = ("https://mapservices.pasda.psu.edu/server/rest/services/" + "pasda/ChesterCounty/MapServer/12/query") +MATCH_RADIUS_M = 50.0 # max centroid→address-point distance to accept +BBOX_MARGIN_DEG = 0.002 # pad the query bbox so edge homes still match +RESIDENTIAL_MAX_FOOTPRINT_M2 = 500 +UNIT_MARKERS = (" APT", " STE", " UNIT", " FL ", " #", " BLDG", " RM ") + + +def _fetch_address_points(bbox): + """All address points intersecting bbox: an id-only query then batched POSTs + (robust to the service's 1000-record cap and URL-length limits).""" + west, south, east, north = bbox + geom = {"geometry": f"{west},{south},{east},{north}", + "geometryType": "esriGeometryEnvelope", "inSR": "4326", + "spatialRel": "esriSpatialRelIntersects"} + ids = requests.get(ADDR_LAYER, params={**geom, "where": "1=1", + "returnIdsOnly": "true", "f": "json"}, timeout=90) \ + .json().get("objectIds") or [] + print(f" {len(ids)} address points in region; fetching...") + + points = [] + for i in range(0, len(ids), 200): + resp = requests.post(ADDR_LAYER, data={ + "objectIds": ",".join(map(str, ids[i:i + 200])), + "outFields": "ADDR_NUM,ROAD_NAME,FULL_ADDRE,MUNI_NAME,ZIP_CODE_S,TAX_EXEMPT", + "returnGeometry": "true", "outSR": "4326", "f": "json"}, timeout=90) + resp.raise_for_status() + for f in resp.json().get("features", []): + g, a = f.get("geometry"), f.get("attributes", {}) + if g: + points.append((g["x"], g["y"], a)) + return points + + +def _projector(lat0): + """Local equirectangular lon/lat → meters, for fast nearest-neighbour matching.""" + cos0 = np.cos(np.radians(lat0)) + return lambda lon, lat: (np.asarray(lon) * 111320.0 * cos0, + np.asarray(lat) * 110540.0) + + +def enrich(cfg: Config) -> None: + try: + con = duckdb.connect(cfg.db_path) + except duckdb.IOException: + print("Database is locked — is the pipeline still running? Let it finish, " + "then re-run this (it's a separate step).") + return + homes = con.execute("SELECT osm_id, lat, lon, footprint_m2 FROM homes").fetchdf() + if homes.empty: + print("No homes in the database yet — run the pipeline first.") + con.close() + return + print(f"Enriching {len(homes)} scored homes with addresses...") + + bbox = (homes["lon"].min() - BBOX_MARGIN_DEG, homes["lat"].min() - BBOX_MARGIN_DEG, + homes["lon"].max() + BBOX_MARGIN_DEG, homes["lat"].max() + BBOX_MARGIN_DEG) + pts = _fetch_address_points(bbox) + footprints = dict(zip(homes["osm_id"].values, homes["footprint_m2"].values)) + rows, matched = [], 0 + + if not pts: + print(" No county address points for this region — classifying residential " + "by footprint only (street addresses need that county's service).") + for osm_id in homes["osm_id"].values: + osm_id = int(osm_id) + foot = footprints.get(osm_id) + small = foot is None or foot < RESIDENTIAL_MAX_FOOTPRINT_M2 + rows.append((osm_id, None, None, None, None, None, None, small)) + else: + proj = _projector(float(homes["lat"].median())) + ax, ay = proj([p[0] for p in pts], [p[1] for p in pts]) + tree = cKDTree(np.column_stack([ax, ay])) + hx, hy = proj(homes["lon"].values, homes["lat"].values) + dist, idx = tree.query(np.column_stack([hx, hy]), distance_upper_bound=MATCH_RADIUS_M) + for osm_id, d, j in zip(homes["osm_id"].values, dist, idx): + osm_id = int(osm_id) + foot = footprints.get(osm_id) + small = foot is None or foot < RESIDENTIAL_MAX_FOOTPRINT_M2 + if np.isinf(d): + rows.append((osm_id, None, None, None, None, None, None, small)) + continue + matched += 1 + a = pts[j][2] + full = a.get("FULL_ADDRE") or "" + tax = a.get("TAX_EXEMPT") + has_unit = any(mk in f" {full.upper()}" for mk in UNIT_MARKERS) + is_res = (tax != "Yes") and (not has_unit) and small + rows.append((osm_id, full or None, a.get("ADDR_NUM"), a.get("ROAD_NAME"), + a.get("MUNI_NAME"), a.get("ZIP_CODE_S"), tax, is_res)) + + for col, typ in [("full_address", "VARCHAR"), ("house_number", "VARCHAR"), + ("street", "VARCHAR"), ("city", "VARCHAR"), ("zip_code", "VARCHAR"), + ("tax_exempt", "VARCHAR"), ("is_residential", "BOOLEAN")]: + try: + con.execute(f"ALTER TABLE homes ADD COLUMN {col} {typ}") + except duckdb.Error: + pass # already exists on a re-run + + con.register("addr_df", pd.DataFrame(rows, columns=[ + "osm_id", "full_address", "house_number", "street", "city", "zip_code", + "tax_exempt", "is_residential"])) + con.execute(""" + UPDATE homes SET full_address = a.full_address, house_number = a.house_number, + street = a.street, city = a.city, zip_code = a.zip_code, + tax_exempt = a.tax_exempt, is_residential = a.is_residential + FROM addr_df a WHERE homes.osm_id = a.osm_id + """) + + total = len(homes) + n_res = con.execute("SELECT COUNT(*) FROM homes WHERE is_residential").fetchone()[0] + print(f"\n Matched {matched} / {total} homes ({100 * matched / total:.0f}%) to an address.") + print(f" Residential (single-family): {n_res} / {total} ({100 * n_res / total:.0f}%).") + con.close() diff --git a/solargrader/geometry.py b/solargrader/geometry.py new file mode 100644 index 0000000..d747f73 --- /dev/null +++ b/solargrader/geometry.py @@ -0,0 +1,88 @@ +""" +Pure geometry / coordinate helpers — no I/O, no heavy dependencies. + +The pipeline mixes two coordinate systems: building footprints arrive in lon/lat +(EPSG:4326) while LiDAR tiles are in projected meters (a UTM zone). These functions +handle the reprojection and the clip of a loaded point cloud to one footprint. +""" + +from __future__ import annotations + +import math + +import numpy as np +import shapely +from pyproj import Transformer +from shapely.geometry import Polygon +from shapely.ops import transform as shapely_transform + +from .config import Config + +# A degree of latitude ≈ 111,320 m; a degree of longitude ≈ that × cos(lat). +_M_PER_DEG = 111320.0 + +# LiDAR point tuple: (structured array, X float view, Y float view). Loading a tile +# once and masking in RAM is far faster than re-reading the file per building. +TilePoints = tuple[np.ndarray, np.ndarray, np.ndarray] + + +def polygon_area_m2(poly: Polygon) -> float: + """Approximate a lon/lat polygon's area in m² using a local equal-area scale.""" + lat = poly.centroid.y + deg2_to_m2 = (_M_PER_DEG ** 2) * math.cos(math.radians(lat)) + return poly.area * deg2_to_m2 + + +def transformer_to(tile_srs: str) -> Transformer: + """Transformer from lon/lat (EPSG:4326) into a tile's projected CRS.""" + return Transformer.from_crs("EPSG:4326", tile_srs, always_xy=True) + + +def clip_lidar_to_building(tile_points: TilePoints, footprint_poly: Polygon, + transformer: Transformer, cfg: Config, + buffer_m: float = 0.5) -> np.ndarray | None: + """Extract roof-ish LiDAR points inside one building footprint. + + The footprint (lon/lat) is reprojected into the tile CRS, a fast bounding-box + mask narrows the candidates, then an exact point-in-polygon test selects the + building's points. Building-classified returns (class 6) are preferred; failing + that, points more than 2 m above the local minimum are used. Returns an (N, 3) + array of [X, Y, Z] in tile meters, or ``None`` if too few points are found. + """ + pts, X, Y = tile_points + proj_poly = shapely_transform(transformer.transform, footprint_poly).buffer(buffer_m) + + minx, miny, maxx, maxy = proj_poly.bounds + bbox_mask = (X >= minx) & (X <= maxx) & (Y >= miny) & (Y <= maxy) + if not bbox_mask.any(): + return None + + sub = pts[bbox_mask] + inside = shapely.contains_xy(proj_poly, sub["X"], sub["Y"]) + sub = sub[inside] + if len(sub) == 0: + return None + + if "Classification" in sub.dtype.names: + roof_pts = sub[sub["Classification"] == 6] + if len(roof_pts) < cfg.min_roof_points: + z_min = sub["Z"].min() + roof_pts = sub[sub["Z"] > z_min + 2.0] + else: + z_min = sub["Z"].min() + roof_pts = sub[sub["Z"] > z_min + 2.0] + + if len(roof_pts) < cfg.min_roof_points: + return None + + return np.column_stack([roof_pts["X"], roof_pts["Y"], roof_pts["Z"]]) + + +def tile_latlon_bbox(tile_points: TilePoints, tile_srs: str) -> tuple[float, float, float, float]: + """Derive a tile's lon/lat bbox (w, s, e, n) from its actual point extent.""" + _, X, Y = tile_points + inv = Transformer.from_crs(tile_srs, "EPSG:4326", always_xy=True) + xs = [X.min(), X.max(), X.min(), X.max()] + ys = [Y.min(), Y.min(), Y.max(), Y.max()] + lons, lats = inv.transform(xs, ys) + return (min(lons), min(lats), max(lons), max(lats)) diff --git a/solargrader/grading.py b/solargrader/grading.py new file mode 100644 index 0000000..6899af2 --- /dev/null +++ b/solargrader/grading.py @@ -0,0 +1,77 @@ +""" +Grading policy — pure functions turning physical results into sales grades. + +Two independent axes: + * ``grade_home`` residential lead quality (A+…D), driven by specific yield. + * ``grade_potential`` maximum-roof capacity tier (P1…P5), the upside beyond a + standard residential system. + +Thresholds are module-level constants so the scoring policy is easy to read and +tune in one place. Re-grading a stored DB never needs LiDAR reprocessing — see +``solargrader.regrade``. +""" + +from __future__ import annotations + +# Residential grade: specific-yield tiers (kWh per kW installed) and their points. +# Yield already integrates orientation, tilt, and shading into one physical number, +# so it is the dominant axis. PA rooftops span ~900–1500 kWh/kW. +_YIELD_POINTS = [(1400, 65), (1300, 52), (1200, 38), (1100, 24), (1000, 12)] +_YIELD_FLOOR = 4 + +# Orientation reinforces prime south-facing among similar-yield roofs. +_ORIENT_POINTS = [(15, 25), (45, 17), (75, 7)] + +# Letter-grade cutoffs on the 0–100 score. +_GRADE_CUTOFFS = [(85, "A+"), (70, "A"), (55, "B+"), (40, "B"), (25, "C")] + +# Potential tier cutoffs on max installable kW (DC). +_POTENTIAL_CUTOFFS = [(30, "P1"), (20, "P2"), (12, "P3"), (7, "P4")] + + +def grade_home(res_annual_kwh: float, res_system_kw: float, + primary_azimuth_deg: float) -> tuple[str, int]: + """Residential SALES grade (letter, 0–100 score). + + Driven by SPECIFIC YIELD (kWh/kW) rather than total energy: the residential + system is capped, so total kWh saturates and stops discriminating (it collapsed + every home to A/A+). Yield restores the spread — a south-facing unshaded PA roof + lands ~1400–1500, a poorly-oriented or shaded one ~900–1000. Roof SIZE is the + other grade's job (``grade_potential``), so it enters here only as a small-system + penalty, never a reward the capped majority would all collect. + """ + yield_kwh_per_kw = res_annual_kwh / res_system_kw if res_system_kw else 0.0 + score = _YIELD_FLOOR + for threshold, points in _YIELD_POINTS: + if yield_kwh_per_kw >= threshold: + score = points + break + + offset = abs(primary_azimuth_deg - 180) + for max_offset, points in _ORIENT_POINTS: + if offset <= max_offset: + score += points + break + + # Deal size: only penalize sub-residential roofs; don't reward the capped majority. + if res_system_kw >= 7: + score += 10 + elif res_system_kw >= 5: + score += 5 + elif res_system_kw < 3: + score -= 10 + + for cutoff, letter in _GRADE_CUTOFFS: + if score >= cutoff: + return letter, score + return "D", score + + +def grade_potential(max_kw: float) -> str: + """Tier a roof by its MAXIMUM installable capacity (kW DC) — the upside beyond a + standard residential system (battery, EV, expansion, or commercial). P1 = very + high (large/commercial roof) … P5 = barely fits a full residential set.""" + for cutoff, tier in _POTENTIAL_CUTOFFS: + if max_kw >= cutoff: + return tier + return "P5" diff --git a/solargrader/models.py b/solargrader/models.py new file mode 100644 index 0000000..3ee8ef0 --- /dev/null +++ b/solargrader/models.py @@ -0,0 +1,68 @@ +""" +Domain data models — lightweight ``@dataclass`` records passed between pipeline +stages. These replace opaque dicts so signatures are self-documenting and typed, +without adding behaviour-heavy classes. +""" + +from __future__ import annotations + +from dataclasses import dataclass + +from shapely.geometry import Polygon + + +@dataclass +class Building: + """A building footprint from any source (OSM / county / Microsoft). + + ``id`` is a stable, source-unique integer (kept in the DB as ``osm_id`` for + backward compatibility). ``geometry`` is in lon/lat (EPSG:4326). + """ + id: int + geometry: Polygon + lat: float + lon: float + footprint_m2: float + + +@dataclass +class RoofPlane: + """One planar roof face fitted from LiDAR.""" + tilt_deg: float # 0 = flat, 90 = vertical + azimuth_deg: float # compass bearing the plane faces (0=N, 90=E, 180=S) + area_m2: float # true tilted surface area + point_count: int + + +@dataclass +class SystemResult: + """Energy model output for one roof: a residential (capped) system and the + maximum-roof system, plus the shading loss backed out from the skyline.""" + res_kwh: float + res_kw: float + res_panels: int + max_kwh: float + max_kw: float + max_panels: int + shade_loss_pct: float + best_plane: RoofPlane + + +@dataclass +class ScoredHome: + """A fully scored building — the row persisted to the ``homes`` table.""" + building: Building + usable_area_m2: float + primary_tilt_deg: float + primary_azimuth_deg: float + plane_count: int + res_kwh: float + res_kw: float + res_panels: int + max_kwh: float + max_kw: float + max_panels: int + potential_grade: str + shade_loss_pct: float + grade: str + score: int diff --git a/solargrader/pipeline.py b/solargrader/pipeline.py new file mode 100644 index 0000000..1113204 --- /dev/null +++ b/solargrader/pipeline.py @@ -0,0 +1,197 @@ +""" +Pipeline orchestration — streaming, parallel, and resumable. + +Tiles are independent, so each is scored end-to-end in its own worker process; the +main process holds the single DuckDB writer and persists results as workers finish. +Each worker deletes its own tile after scoring, so peak disk stays ~``n_workers`` +tiles. ``Config`` is passed explicitly into every worker (it is picklable), so a run +is fully described by that object even under the spawn start method. +""" + +from __future__ import annotations + +import math +import os +import time +from concurrent.futures import ProcessPoolExecutor, as_completed + +from shapely.ops import transform as shapely_transform + +from .config import Config +from .geometry import clip_lidar_to_building, tile_latlon_bbox, transformer_to +from .grading import grade_home, grade_potential +from .models import Building, ScoredHome +from .roof import extract_roof_planes +from .solar import annual_energy, compute_horizon +from .sources.buildings import get_buildings +from .sources.irradiance import get_tmy +from .sources.tiles import discover_tiles, download_tile, get_tile_srs, load_tile_points, tile_id_of +from .storage import ResultStore +from .util import with_retries + + +def score_building(building: Building, tile_points, transformer, tmy, + cfg: Config) -> ScoredHome | None: + """Full scoring chain for one building: clip LiDAR → fit planes → shade → energy + → grade. Returns a ``ScoredHome`` or ``None`` if the building can't be scored.""" + points = clip_lidar_to_building(tile_points, building.geometry, transformer, cfg) + if points is None: + return None + planes = extract_roof_planes(points, cfg) + if not planes: + return None + + observer = points.mean(axis=0) + own_footprint = shapely_transform(transformer.transform, building.geometry).buffer(1.0) + horizon = compute_horizon(tile_points, observer, cfg, exclude_poly=own_footprint) + + sysres = annual_energy(building.lat, building.lon, tmy, planes, cfg, horizon) + if sysres is None or sysres.res_panels == 0: + return None + + best = sysres.best_plane + grade, score = grade_home(sysres.res_kwh, sysres.res_kw, best.azimuth_deg) + return ScoredHome( + building=building, + usable_area_m2=sum(p.area_m2 for p in planes), + primary_tilt_deg=best.tilt_deg, primary_azimuth_deg=best.azimuth_deg, + plane_count=len(planes), + res_kwh=sysres.res_kwh, res_kw=sysres.res_kw, res_panels=sysres.res_panels, + max_kwh=sysres.max_kwh, max_kw=sysres.max_kw, max_panels=sysres.max_panels, + potential_grade=grade_potential(sysres.max_kw), + shade_loss_pct=sysres.shade_loss_pct, grade=grade, score=score, + ) + + +def _delete_tile(path: str, cfg: Config) -> None: + if cfg.delete_tiles_after and os.path.exists(path): + os.remove(path) + + +def _process_tile(tile: dict, cfg: Config) -> dict: + """Worker body: fully process ONE tile in its own process, then delete the tile + file. Returns a compact, picklable result and never touches the DB. All failures + are caught and returned so one bad tile can't kill the pool.""" + tid = tile_id_of(tile) + try: + path = download_tile(tile["downloadURL"], cfg.tile_cache_dir) + tile_points = load_tile_points(path) + tile_srs = get_tile_srs(path) + transformer = transformer_to(tile_srs) + tbbox = tile_latlon_bbox(tile_points, tile_srs) + + buildings = with_retries(get_buildings, tbbox, cfg) + if not buildings: + _delete_tile(path, cfg) + return {"tid": tid, "n_buildings": 0, "scored": [], "error": None} + + # One irradiance fetch per tile (its center) covers all its buildings. + cx, cy = (tbbox[0] + tbbox[2]) / 2, (tbbox[1] + tbbox[3]) / 2 + tmy = with_retries(get_tmy, cy, cx) + + scored = [] + for building in buildings: + home = score_building(building, tile_points, transformer, tmy, cfg) + if home is not None: + # Geometry isn't needed downstream (save() uses id/lat/lon/footprint); + # drop it so the inter-process payload stays small. + home.building.geometry = None + scored.append(home) + + _delete_tile(path, cfg) + return {"tid": tid, "n_buildings": len(buildings), "scored": scored, "error": None} + except Exception as e: # noqa: BLE001 + return {"tid": tid, "n_buildings": 0, "scored": [], + "error": f"{type(e).__name__}: {e}"} + + +def run(bbox, cfg: Config, store: ResultStore | None = None) -> None: + """Score every tile covering ``bbox``, accumulating into the results DB.""" + print("=" * 60) + print(f"SOLAR GRADER — region {tuple(round(v, 3) for v in bbox)}") + print("=" * 60) + print("\nDiscovering LiDAR tiles...") + tiles = discover_tiles(bbox, cfg) + if not tiles: + print(" ABORT: no tiles for this region (TNM + fallback both empty).") + return + + own_store = store is None + store = store or ResultStore(cfg.db_path) + try: + done = store.done_tiles() + pending = [t for t in tiles if tile_id_of(t) not in done] + if done: + print(f" {len(done)} already processed — resuming.") + print(f"\nScoring {len(pending)} tiles across {cfg.n_workers} worker " + f"processes (single DB writer in this process)...") + + t_start = time.time() + completed = 0 + with ProcessPoolExecutor(max_workers=cfg.n_workers) as pool: + futures = {pool.submit(_process_tile, t, cfg): tile_id_of(t) for t in pending} + for fut in as_completed(futures): + completed += 1 + res = fut.result() # workers catch their own errors; never raises here + tid = res["tid"] + if res["error"]: + print(f"[{completed}/{len(pending)}] {tid} — ERROR: {res['error']} (skipped)") + continue + for home in res["scored"]: + store.save(home) + store.mark_tile_done(tid, len(res["scored"])) + print(f"[{completed}/{len(pending)}] {tid} — " + f"scored {len(res['scored'])}/{res['n_buildings']}") + + if own_store: + store.print_summary(time.time() - t_start) + finally: + if own_store: + store.close() + + +def run_region(bbox, cfg: Config, store: ResultStore | None = None) -> None: + """Process a whole county/state by tiling ``bbox`` into sub-regions and running + each. Resumable: a ``regions_done`` ledger skips finished sub-regions, and + ``tiles_done`` (inside ``run``) dedups tiles shared across sub-regions.""" + w, s, e, n = bbox + step = cfg.step_deg + nx = max(1, math.ceil((e - w) / step)) # integer step counts avoid float drift + ny = max(1, math.ceil((n - s) / step)) + subs = [] + for ix in range(nx): + x0, x1 = w + ix * step, min(w + (ix + 1) * step, e) + for iy in range(ny): + y0, y1 = s + iy * step, min(s + (iy + 1) * step, n) + if x1 - x0 > 1e-9 and y1 - y0 > 1e-9: + subs.append((round(x0, 4), round(y0, 4), round(x1, 4), round(y1, 4))) + + print("#" * 60) + print(f"REGION RUN — {len(subs)} sub-regions (~{step}°) over " + f"{tuple(round(v, 3) for v in bbox)}") + print("#" * 60) + + own_store = store is None + store = store or ResultStore(cfg.db_path) + try: + done = store.done_regions() + if done: + print(f" Resuming — {len(done)} sub-regions already complete.") + for i, sub in enumerate(subs, 1): + key = ",".join(f"{v:.4f}" for v in sub) + if key in done: + print(f"\n### sub-region {i}/{len(subs)} {sub} — already done, skip.") + continue + print(f"\n### sub-region {i}/{len(subs)} {sub}") + try: + run(sub, cfg, store=store) + except Exception as ex: # noqa: BLE001 + print(f"### sub-region {i} ERROR: {type(ex).__name__}: {ex} — continuing.") + continue + store.mark_region_done(key) + print("\n" + "#" * 60) + print(f"REGION RUN COMPLETE — {store.count_homes()} homes scored.") + print("#" * 60) + finally: + if own_store: + store.close() diff --git a/solargrader/regrade.py b/solargrader/regrade.py new file mode 100644 index 0000000..d3ea2d7 --- /dev/null +++ b/solargrader/regrade.py @@ -0,0 +1,53 @@ +""" +Re-grade stored homes in place — no LiDAR reprocessing. + +Grading is a pure function of columns already persisted (yield, kW, azimuth, max +kW), so tuning the grading policy never requires re-running the pipeline. This reads +every row, re-applies ``grade_home`` / ``grade_potential``, and reports the shift. +""" + +from __future__ import annotations + +import duckdb +import pandas as pd + +from .config import Config +from .grading import grade_home, grade_potential + +_GRADE_ORDER = ["A+", "A", "B+", "B", "C", "D"] + + +def regrade(cfg: Config) -> None: + con = duckdb.connect(cfg.db_path) + homes = con.execute( + "SELECT osm_id, res_annual_kwh, res_system_kw, primary_azimuth_deg, " + "max_system_kw, solar_grade FROM homes" + ).fetchdf() + if homes.empty: + print("No homes to re-grade — run the pipeline first.") + con.close() + return + + before = homes["solar_grade"].value_counts().to_dict() + rows = [] + for r in homes.itertuples(index=False): + grade, score = grade_home(r.res_annual_kwh, r.res_system_kw, r.primary_azimuth_deg) + rows.append((int(r.osm_id), grade, int(score), grade_potential(r.max_system_kw))) + + con.register("new_df", pd.DataFrame(rows, columns=["osm_id", "grade", "score", "pot"])) + con.execute(""" + UPDATE homes SET solar_grade = n.grade, solar_score = n.score, + potential_grade = n.pot + FROM new_df n WHERE homes.osm_id = n.osm_id + """) + + total = len(homes) + after = con.execute("SELECT solar_grade, COUNT(*) FROM homes GROUP BY 1").fetchdf() \ + .set_index("solar_grade")["count_star()"].to_dict() + con.close() + + print(f"Re-graded {total} homes.\n") + print(f" {'grade':>5} {'before':>8} {'after':>8}") + for g in _GRADE_ORDER: + b, a = before.get(g, 0), after.get(g, 0) + print(f" {g:>5} {b:>8} {a:>8} ({100 * a / total:>4.1f}%)") diff --git a/solargrader/roof.py b/solargrader/roof.py new file mode 100644 index 0000000..6baf69d --- /dev/null +++ b/solargrader/roof.py @@ -0,0 +1,96 @@ +""" +Roof-plane extraction from a building's LiDAR points via RANSAC (Open3D). + +Pure with respect to the pipeline — takes points + config, returns ``RoofPlane`` +objects. Open3D and SciPy are imported lazily so the rest of the package imports +without them. +""" + +from __future__ import annotations + +import numpy as np + +from .config import Config +from .models import RoofPlane + +# Surfaces steeper than this are walls/dormers, not paneling candidates. +_MAX_TILT_DEG = 72.0 +# Roughly the maximum number of distinct faces on a residential roof. +_MAX_PLANES = 6 + + +def _ground_density(points_xyz: np.ndarray) -> float: + """Points per m² of horizontal area, from ONE convex hull over all points. + + Sizing each plane by ``point_count / density`` avoids the area inflation you get + from summing overlapping per-plane hulls: every point belongs to exactly one + plane, so the summed areas can't exceed the building's real footprint. + """ + from scipy.spatial import ConvexHull # lazy: SciPy is heavy-ish + + all_xy = np.asarray(points_xyz)[:, :2] + try: + covered_m2 = float(ConvexHull(all_xy).volume) if len(all_xy) >= 3 else 0.0 + except Exception: + covered_m2 = 0.0 + return (len(points_xyz) / covered_m2) if covered_m2 > 0 else 0.0 + + +def extract_roof_planes(points_xyz: np.ndarray, cfg: Config) -> list[RoofPlane]: + """Fit up to a handful of roof planes to an (N, 3) point cloud (tile meters). + + Iteratively segments the dominant plane, records its tilt/azimuth/area, removes + its inliers, and repeats. Near-vertical and tiny planes are discarded. + """ + import open3d as o3d # lazy: Open3D is a large optional dependency + + pcd = o3d.geometry.PointCloud() + pcd.points = o3d.utility.Vector3dVector(points_xyz) + + ground_density = _ground_density(points_xyz) + + planes: list[RoofPlane] = [] + remaining = pcd + + for _ in range(_MAX_PLANES): + if len(remaining.points) < cfg.min_roof_points: + break + + plane_model, inliers = remaining.segment_plane( + distance_threshold=cfg.ransac_distance_threshold, + ransac_n=3, + num_iterations=cfg.ransac_iterations, + ) + if len(inliers) < cfg.min_roof_points: + break + + # Plane equation ax + by + cz + d = 0 → normal = [a, b, c]; force it upward. + normal = np.array(plane_model[:3]) + if normal[2] < 0: + normal = -normal + norm_len = np.linalg.norm(normal) + if norm_len == 0: + remaining = remaining.select_by_index(inliers, invert=True) + continue + nx, ny, nz = normal / norm_len + + tilt_deg = float(np.degrees(np.arccos(min(nz, 1.0)))) + if tilt_deg > _MAX_TILT_DEG: + remaining = remaining.select_by_index(inliers, invert=True) + continue + + # Compass azimuth: +X=east, +Y=north → atan2(nx, ny), 0=N, 90=E, 180=S. + azimuth_deg = float((np.degrees(np.arctan2(nx, ny)) + 360) % 360) + + # Horizontal area from point count ÷ density, lifted to the tilted surface. + if ground_density > 0: + ground_area_m2 = len(inliers) / ground_density + area_m2 = ground_area_m2 / max(np.cos(np.radians(tilt_deg)), 0.1) + else: + area_m2 = 0.0 + + planes.append(RoofPlane(tilt_deg=tilt_deg, azimuth_deg=azimuth_deg, + area_m2=area_m2, point_count=len(inliers))) + remaining = remaining.select_by_index(inliers, invert=True) + + return planes diff --git a/solargrader/solar.py b/solargrader/solar.py new file mode 100644 index 0000000..e85357b --- /dev/null +++ b/solargrader/solar.py @@ -0,0 +1,145 @@ +""" +Solar energy model — pure functions over roof planes and irradiance. + +Two systems are reported per roof: the *maximum* (every panel that fits) and the +*residential* (best-producing panels up to the kW cap — the sellable system). If a +LiDAR-derived skyline is supplied, the direct beam is removed during hours the sun +sits below an obstruction, and the lost fraction is reported as ``shade_loss_pct``. +""" + +from __future__ import annotations + +import numpy as np + +from .config import Config +from .geometry import TilePoints +from .models import RoofPlane, SystemResult + + +def panels_on_plane(plane_area_m2: float, cfg: Config) -> int: + """Whole panels that physically fit on a roof plane of this area. + + The plane is approximated as an equal-area square, eroded inward by the fire-code + setback on all sides, packed at ``panel_packing``, then divided by module area. + Small/cut-up roofs lose proportionally more — behaviour a flat usable-% can't + capture. + """ + side = np.sqrt(plane_area_m2) + usable_side = max(0.0, side - 2 * cfg.roof_setback_m) + usable_area = (usable_side ** 2) * cfg.panel_packing + return int(usable_area // cfg.panel_area_m2) + + +def compute_horizon(tile_points: TilePoints, observer_xyz, cfg: Config, + exclude_poly=None) -> np.ndarray: + """Build the roof's skyline: the maximum obstruction elevation (degrees above the + roof) per compass direction, from surrounding LiDAR points (trees, neighbours) + within ``shade_radius_m``. Points inside ``exclude_poly`` (the building's OWN + footprint, in tile CRS) are dropped so a roof can't shade itself. Returns an + array indexed by azimuth bin (0=N, clockwise); all zeros = wide-open sky. + """ + import shapely # lazy: only needed when excluding the own footprint + + pts, X, Y = tile_points + ox, oy, oz = observer_xyz + n_bins = cfg.horizon_bins + radius_m = cfg.shade_radius_m + horizon = np.zeros(n_bins) + + m = ((X >= ox - radius_m) & (X <= ox + radius_m) + & (Y >= oy - radius_m) & (Y <= oy + radius_m)) + if not m.any(): + return horizon + + sx, sy = X[m], Y[m] + sz = np.asarray(pts["Z"], dtype=float)[m] + if exclude_poly is not None: + outside = ~shapely.contains_xy(exclude_poly, sx, sy) + sx, sy, sz = sx[outside], sy[outside], sz[outside] + + dx, dy, dz = sx - ox, sy - oy, sz - oz + horiz = np.sqrt(dx * dx + dy * dy) + keep = (horiz <= radius_m) & (horiz > 1.0) & (dz > 0) # within radius, above roof + if not keep.any(): + return horizon + + dx, dy, dz, horiz = dx[keep], dy[keep], dz[keep], horiz[keep] + elev = np.degrees(np.arctan2(dz, horiz)) # elevation angle + az = (np.degrees(np.arctan2(dx, dy)) + 360) % 360 # 0=N, 90=E (compass) + bins = (az / (360.0 / n_bins)).astype(int) % n_bins + np.maximum.at(horizon, bins, elev) # tallest per direction + return horizon + + +def annual_energy(lat: float, lon: float, tmy, planes: list[RoofPlane], cfg: Config, + horizon: np.ndarray | None = None) -> SystemResult | None: + """Pack panels onto every viable plane, compute plane-of-array irradiance, and + return the residential + maximum systems (or ``None`` if no plane is viable). + """ + import pvlib # lazy: pvlib pulls pandas machinery + + location = pvlib.location.Location(latitude=lat, longitude=lon, tz="Etc/GMT+5") + solar_pos = location.get_solarposition(tmy.index) + + # Per-hour direct-beam shading mask: sun up but below the roof's skyline. + sun_blocked = None + if horizon is not None: + sun_elev = 90.0 - solar_pos["apparent_zenith"].values + sun_az = solar_pos["azimuth"].values + sun_bins = (sun_az / (360.0 / len(horizon))).astype(int) % len(horizon) + sun_blocked = (sun_elev > 0) & (sun_elev <= horizon[sun_bins]) + + viable = [] # (per_panel_shaded, per_panel_unshaded, n_panels, plane) + for plane in planes: + if abs(plane.azimuth_deg - 180) > cfg.max_azimuth_offset: + continue + n_panels = panels_on_plane(plane.area_m2, cfg) + if n_panels <= 0: + continue + + poa = pvlib.irradiance.get_total_irradiance( + surface_tilt=plane.tilt_deg, + surface_azimuth=plane.azimuth_deg, + solar_zenith=solar_pos["apparent_zenith"], + solar_azimuth=solar_pos["azimuth"], + ghi=tmy["ghi"], dni=tmy["dni"], dhi=tmy["dhi"], + ) + glob = poa["poa_global"].fillna(0).values + unshaded_annual = glob.sum() / 1000 # kWh/m²/yr, full sun + if sun_blocked is not None: + direct = poa["poa_direct"].fillna(0).values + shaded_annual = (glob - direct * sun_blocked).sum() / 1000 + else: + shaded_annual = unshaded_annual + + per_kwh = cfg.panel_area_m2 * cfg.module_efficiency * cfg.system_losses + viable.append((shaded_annual * per_kwh, unshaded_annual * per_kwh, n_panels, plane)) + + if not viable: + return None + + # Max system: every panel that physically fits the roof. + max_panels = sum(n for _, _, n, _ in viable) + max_kwh = sum(s * n for s, _, n, _ in viable) + max_unshaded = sum(u * n for _, u, n, _ in viable) + max_kw = max_panels * cfg.panel_watts / 1000.0 + shade_loss_pct = 100.0 * (1 - max_kwh / max_unshaded) if max_unshaded > 0 else 0.0 + + # Residential system: fill the best-producing panels first, up to the kW cap. + cap_panels = cfg.residential_panel_cap + res_panels = 0 + res_kwh = 0.0 + for ppk_shaded, _, n_panels, _ in sorted(viable, key=lambda v: v[0], reverse=True): + take = min(n_panels, cap_panels - res_panels) + if take <= 0: + break + res_panels += take + res_kwh += take * ppk_shaded + res_kw = res_panels * cfg.panel_watts / 1000.0 + + best_plane = max(viable, key=lambda v: v[0])[3] + return SystemResult( + res_kwh=res_kwh, res_kw=res_kw, res_panels=res_panels, + max_kwh=max_kwh, max_kw=max_kw, max_panels=max_panels, + shade_loss_pct=shade_loss_pct, best_plane=best_plane, + ) diff --git a/solargrader/sources/__init__.py b/solargrader/sources/__init__.py new file mode 100644 index 0000000..964174c --- /dev/null +++ b/solargrader/sources/__init__.py @@ -0,0 +1,7 @@ +"""Data acquisition — all free public sources, each in its own module. + + tiles USGS 3DEP LiDAR tile discovery + download + collection selection + buildings building footprints (Microsoft / county / OpenStreetMap) + irradiance EU PVGIS typical-meteorological-year data + regions county/state → bounding box (US Census TIGERweb) +""" diff --git a/solargrader/sources/buildings.py b/solargrader/sources/buildings.py new file mode 100644 index 0000000..cb33060 --- /dev/null +++ b/solargrader/sources/buildings.py @@ -0,0 +1,208 @@ +""" +Building footprint sources, all returning ``list[Building]`` so scoring is identical +regardless of provider: + + "ms" Microsoft Building Footprints — national, ODbL (the scale-out source) + "county" a county ArcGIS service (Chester County / PASDA by default) + "osm" OpenStreetMap via Overpass (handy for quick prototyping) +""" + +from __future__ import annotations + +import gzip +import json +import math +import os + +import numpy as np +import requests +from shapely.geometry import Polygon + +from ..config import Config +from ..geometry import polygon_area_m2 +from ..models import Building + +# ── OpenStreetMap (Overpass) ────────────────────────────────────────────────── + +def get_buildings_from_osm(bbox, cfg: Config) -> list[Building]: + """Building polygons from the Overpass API (prototyping source).""" + west, south, east, north = bbox + query = (f"[out:json][timeout:60];(way[\"building\"]" + f"({south},{west},{north},{east}););out geom;") + r = requests.post(cfg.overpass_url, data={"data": query}, + headers={"User-Agent": cfg.http_user_agent}, timeout=90) + r.raise_for_status() + + buildings = [] + for el in r.json().get("elements", []): + if el["type"] != "way" or "geometry" not in el: + continue + coords = [(n["lon"], n["lat"]) for n in el["geometry"]] + if len(coords) < 4: + continue + try: + poly = Polygon(coords) + if not poly.is_valid or poly.area == 0: + continue + except Exception: # noqa: BLE001 + continue + buildings.append(Building(id=el["id"], geometry=poly, + lat=poly.centroid.y, lon=poly.centroid.x, + footprint_m2=polygon_area_m2(poly))) + print(f" Found {len(buildings)} buildings (OSM).") + return buildings + + +# ── County ArcGIS service (PASDA / Chester County by default) ───────────────── + +def get_buildings_from_county(bbox, cfg: Config) -> list[Building]: + """County building footprints via ArcGIS REST. An id-only query then batched + geometry fetches keeps us robust to the service's 1000-record cap.""" + west, south, east, north = bbox + geom = {"geometry": f"{west},{south},{east},{north}", + "geometryType": "esriGeometryEnvelope", "inSR": "4326", + "spatialRel": "esriSpatialRelIntersects"} + print(" Querying county building footprints...") + ids = requests.get(cfg.county_footprints_url, + params={**geom, "where": "1=1", "returnIdsOnly": "true", + "f": "json"}, timeout=90).json().get("objectIds") or [] + + buildings = [] + for i in range(0, len(ids), 200): + batch = ids[i:i + 200] + resp = requests.post(cfg.county_footprints_url, data={ + "objectIds": ",".join(map(str, batch)), "outFields": "OBJECTID", + "returnGeometry": "true", "outSR": "4326", "f": "json"}, timeout=90) + resp.raise_for_status() + for f in resp.json().get("features", []): + rings = (f.get("geometry") or {}).get("rings") + if not rings: + continue + try: + poly = Polygon(rings[0]) + if not poly.is_valid or poly.area == 0: + continue + except Exception: # noqa: BLE001 + continue + buildings.append(Building(id=int(f["attributes"]["OBJECTID"]), geometry=poly, + lat=poly.centroid.y, lon=poly.centroid.x, + footprint_m2=polygon_area_m2(poly))) + print(f" Found {len(buildings)} buildings (county).") + return buildings + + +# ── Microsoft Building Footprints (national, quadkey-partitioned) ───────────── + +def _quadkey(lon: float, lat: float, z: int) -> str: + """Bing Maps quadkey for a lon/lat at zoom z (MS footprint partition scheme).""" + n = 2 ** z + x = int((lon + 180) / 360 * n) + y = int((1 - math.asinh(math.tan(math.radians(lat))) / math.pi) / 2 * n) + qk = "" + for i in range(z, 0, -1): + digit, mask = 0, 1 << (i - 1) + if x & mask: + digit += 1 + if y & mask: + digit += 2 + qk += str(digit) + return qk + + +def _ms_link_for(qk: str, cfg: Config): + """Download URL for a US quadkey in MS's dataset index (index cached on disk).""" + os.makedirs(cfg.ms_cache_dir, exist_ok=True) + links = os.path.join(cfg.ms_cache_dir, "dataset-links.csv") + if not os.path.exists(links): + print(" Fetching MS dataset index (one-time)...") + r = requests.get(cfg.ms_dataset_links_url, timeout=180) + r.raise_for_status() + with open(links, "w") as f: + f.write(r.text) + with open(links) as f: + for line in f: + p = line.split(",") + if len(p) >= 3 and p[0] == "UnitedStates" and p[1] == qk: + return p[2] + return None + + +# quadkey -> (centroids ndarray, list-of-rings); a quadkey covers a wide area and is +# reused across many tiles, so parsing it once per process is worth caching. +_MS_CACHE: dict = {} + + +def _load_ms_quadkey(qk: str, cfg: Config): + """Download (once) and parse a quadkey's GeoJSONL into centroids + rings.""" + if qk in _MS_CACHE: + return _MS_CACHE[qk] + url = _ms_link_for(qk, cfg) + if not url: + _MS_CACHE[qk] = (np.empty((0, 2)), []) + return _MS_CACHE[qk] + path = os.path.join(cfg.ms_cache_dir, f"{qk}.csv.gz") + if not os.path.exists(path): + print(f" Downloading MS footprints quadkey {qk} (once)...") + r = requests.get(url, timeout=300) + r.raise_for_status() + with open(path, "wb") as f: + f.write(r.content) + centroids, rings = [], [] + with gzip.open(path, "rt") as f: + for line in f: + try: + coords = json.loads(line)["geometry"]["coordinates"][0] + except Exception: # noqa: BLE001 + continue + xs = [c[0] for c in coords] + ys = [c[1] for c in coords] + centroids.append((sum(xs) / len(xs), sum(ys) / len(ys))) + rings.append(coords) + _MS_CACHE[qk] = (np.array(centroids), rings) + print(f" Loaded {len(rings):,} MS buildings for quadkey {qk}.") + return _MS_CACHE[qk] + + +def get_buildings_from_ms(bbox, cfg: Config) -> list[Building]: + """Microsoft Building Footprints within bbox (national / other-states source).""" + w, s, e, n = bbox + quadkeys = {_quadkey(lo, la, cfg.ms_zoom) for lo in (w, e) for la in (s, n)} + buildings = [] + for qk in quadkeys: + centroids, rings = _load_ms_quadkey(qk, cfg) + if len(centroids) == 0: + continue + cx, cy = centroids[:, 0], centroids[:, 1] + mask = (cx >= w) & (cx <= e) & (cy >= s) & (cy <= n) + qk_int = int(qk) + for idx in np.nonzero(mask)[0]: + try: + poly = Polygon(rings[idx]) + if not poly.is_valid or poly.area == 0: + continue + except Exception: # noqa: BLE001 + continue + lon, lat = float(cx[idx]), float(cy[idx]) + buildings.append(Building(id=qk_int * 10_000_000 + int(idx), geometry=poly, + lat=lat, lon=lon, footprint_m2=polygon_area_m2(poly))) + print(f" Found {len(buildings)} buildings (MS footprints).") + return buildings + + +# ── Dispatcher ──────────────────────────────────────────────────────────────── + +_SOURCES = { + "osm": get_buildings_from_osm, + "county": get_buildings_from_county, + "ms": get_buildings_from_ms, +} + + +def get_buildings(bbox, cfg: Config) -> list[Building]: + """Fetch buildings from the source named by ``cfg.building_source``.""" + try: + source = _SOURCES[cfg.building_source] + except KeyError: + raise ValueError(f"Unknown building_source {cfg.building_source!r}; " + f"choose one of {sorted(_SOURCES)}") from None + return source(bbox, cfg) diff --git a/solargrader/sources/irradiance.py b/solargrader/sources/irradiance.py new file mode 100644 index 0000000..8e693fb --- /dev/null +++ b/solargrader/sources/irradiance.py @@ -0,0 +1,21 @@ +""" +Solar irradiance source — EU PVGIS typical-meteorological-year (TMY) data. + +PVGIS is free and needs no API key, covers the US, and returns hourly GHI/DNI/DHI. +Isolated behind one function so an alternate source (e.g. NREL NSRDB, which would +use ``Config.require_nrel_api_key``) can be dropped in without touching the model. +""" + +from __future__ import annotations + + +def get_tmy(lat: float, lon: float): + """Fetch a Typical Meteorological Year for a location as a pandas DataFrame with + pvlib-standard columns (ghi, dni, dhi). One fetch covers all buildings nearby.""" + import pvlib # lazy + + tmy, _ = pvlib.iotools.get_pvgis_tmy( # pvlib >=0.11 returns (data, metadata) + latitude=lat, longitude=lon, outputformat="json", + usehorizon=True, map_variables=True, + ) + return tmy diff --git a/solargrader/sources/regions.py b/solargrader/sources/regions.py new file mode 100644 index 0000000..b06e0df --- /dev/null +++ b/solargrader/sources/regions.py @@ -0,0 +1,103 @@ +""" +Area targeting — turn a place (county, named region, or explicit bbox) into the +bounding box(es) a run should cover. County boundaries come from US Census TIGERweb, +so any US county resolves by name with no hand-typed coordinates. +""" + +from __future__ import annotations + +import re + +import requests + +from ..config import Config +from ..util import with_retries + +# A run target: (human label, bbox, whether to grid it into sub-regions). +Target = tuple[str, tuple[float, float, float, float], bool] + +_STATE_FIPS = { + "AL": "01", "AK": "02", "AZ": "04", "AR": "05", "CA": "06", "CO": "08", + "CT": "09", "DE": "10", "DC": "11", "FL": "12", "GA": "13", "HI": "15", + "ID": "16", "IL": "17", "IN": "18", "IA": "19", "KS": "20", "KY": "21", + "LA": "22", "ME": "23", "MD": "24", "MA": "25", "MI": "26", "MN": "27", + "MS": "28", "MO": "29", "MT": "30", "NE": "31", "NV": "32", "NH": "33", + "NJ": "34", "NM": "35", "NY": "36", "NC": "37", "ND": "38", "OH": "39", + "OK": "40", "OR": "41", "PA": "42", "RI": "44", "SC": "45", "SD": "46", + "TN": "47", "TX": "48", "UT": "49", "VT": "50", "VA": "51", "WA": "53", + "WV": "54", "WI": "55", "WY": "56", +} + +# Convenience aliases — type ``--region ``. Each is an explicit bbox or one or +# more counties resolved live via TIGERweb. +REGIONS = { + "harrisburg": {"bbox": (-76.930, 40.240, -76.830, 40.310)}, # city + suburbs + "dauphin": {"counties": ["Dauphin"]}, + "cumberland": {"counties": ["Cumberland"]}, + "perry": {"counties": ["Perry"]}, + "york": {"counties": ["York"]}, + "lancaster": {"counties": ["Lancaster"]}, + "lebanon": {"counties": ["Lebanon"]}, + "harrisburg-metro": {"counties": ["Dauphin", "Cumberland", "Perry"]}, +} + + +def county_bbox(county: str, cfg: Config, state: str = "PA") -> tuple[tuple, str]: + """Exact lon/lat bbox for a US county via TIGERweb. ``county`` is the bare name + (e.g. 'Dauphin', not 'Dauphin County'), validated against a strict whitelist + before it reaches the REST where-clause (trust-boundary hygiene).""" + if not re.fullmatch(r"[A-Za-z][A-Za-z .'\-]{0,48}", county or ""): + raise ValueError(f"Invalid county name: {county!r}") + fips = _STATE_FIPS.get(state.upper()) + if not fips: + raise ValueError(f"Unknown state: {state!r}") + data = with_retries(lambda: requests.get(cfg.tigerweb_counties_url, params={ + "where": f"BASENAME='{county}' AND STATE='{fips}'", "outFields": "NAME", + "returnGeometry": "true", "outSR": "4326", "f": "json"}, timeout=60).json()) + feats = data.get("features") or [] + if not feats: + raise ValueError(f"No county '{county}' found in {state.upper()}.") + rings = feats[0]["geometry"]["rings"] + xs = [p[0] for ring in rings for p in ring] + ys = [p[1] for ring in rings for p in ring] + return (min(xs), min(ys), max(xs), max(ys)), feats[0]["attributes"]["NAME"] + + +def all_counties(cfg: Config, state: str = "PA") -> list[str]: + """Every county name in a state, alphabetized, via TIGERweb (whole-state runs).""" + fips = _STATE_FIPS.get(state.upper()) + if not fips: + raise ValueError(f"Unknown state: {state!r}") + data = with_retries(lambda: requests.get(cfg.tigerweb_counties_url, params={ + "where": f"STATE='{fips}'", "outFields": "BASENAME", "returnGeometry": "false", + "orderByFields": "BASENAME", "f": "json"}, timeout=90).json()) + return [f["attributes"]["BASENAME"] for f in data.get("features", [])] + + +def resolve_targets(cfg: Config, *, bbox: tuple | None = None, county: str | None = None, + region: str | None = None, all_counties_flag: bool = False, + state: str = "PA", grid: bool = False) -> list[Target]: + """Resolve a user's area choice into run targets (label, bbox, use_grid). All + targets accumulate into the same DB; the resume ledgers dedup any overlap.""" + if all_counties_flag: + names = all_counties(cfg, state) + print(f"Resolving {len(names)} counties in {state.upper()}...") + return [(*_county_target(cfg, n, state),) for n in names] + if bbox: + return [("bbox", tuple(bbox), grid)] + if county: + names = [c.strip() for c in county.split(",") if c.strip()] + return [_county_target(cfg, n, state) for n in names] + if region: + spec = REGIONS.get(region.lower()) + if not spec: + raise ValueError(f"Unknown region '{region}'. Options: {', '.join(sorted(REGIONS))}") + if "bbox" in spec: + return [(region, spec["bbox"], grid)] + return [_county_target(cfg, n, spec.get("state", "PA")) for n in spec["counties"]] + raise ValueError("No target given — pass bbox, county, region, or all_counties_flag.") + + +def _county_target(cfg: Config, name: str, state: str) -> Target: + bbox, label = county_bbox(name, cfg, state) + return (label, bbox, True) # a county is always a gridded run diff --git a/solargrader/sources/tiles.py b/solargrader/sources/tiles.py new file mode 100644 index 0000000..aafdc50 --- /dev/null +++ b/solargrader/sources/tiles.py @@ -0,0 +1,287 @@ +""" +LiDAR tile discovery, download, and best-quality collection selection (USGS 3DEP). + +Primary discovery is the TNM products API (naming-agnostic, returns per-tile +download URLs). If that spatial query is down, a region-agnostic fallback uses the +3DEP index to find covering projects and decodes UTM-named tiles from their staged +directories. When multiple overlapping surveys cover a region, ``select_collections`` +builds a best-quality mosaic so scores aren't corrupted by mixing vintages. +""" + +from __future__ import annotations + +import json +import math +import os +import re + +import numpy as np +import requests +from pyproj import Transformer + +from ..config import Config +from ..geometry import TilePoints +from ..util import with_retries + +# USGS LPC filename with a UTM/MGRS grid tail: ..._.laz +_TILE_RE = re.compile(r"(USGS_LPC_[\w.\-]+?_(\d{2})[A-Z]{3}(\d{6}))\.laz") + + +# ── Primary discovery: TNM products API ─────────────────────────────────────── + +def find_all_tiles(bbox, cfg: Config, page: int = 50) -> list[dict]: + """All LiDAR tiles covering bbox from USGS TNM (paginated, de-duped by URL).""" + west, south, east, north = bbox + tiles, offset = {}, 0 + while True: + r = requests.get(cfg.tnm_url, params={ + "datasets": "Lidar Point Cloud (LPC)", + "bbox": f"{west},{south},{east},{north}", + "prodFormats": "LAZ", "max": page, "offset": offset, + }, timeout=60) + r.raise_for_status() + data = r.json() + items = data.get("items", []) + for it in items: + if it.get("downloadURL"): + tiles[it["downloadURL"]] = it + offset += len(items) + if not items or offset >= data.get("total", 0): + break + return list(tiles.values()) + + +def tile_id_of(tile: dict) -> str: + """Stable id for a tile = its LAZ filename.""" + return tile["downloadURL"].split("/")[-1].split("?")[0] + + +# ── Fallback discovery: 3DEP index → rockyweb staged directories ────────────── + +def _covering_cells(bbox) -> tuple[dict, int, int]: + """UTM 1 km cells covering bbox, in the bbox's own UTM zone. Returns + {code -> (easting_km, northing_km)} with the zone and EPSG, where ``code`` is + the 6-digit 'EEENNN' USGS embeds in tile filenames. Works in any UTM zone.""" + w, s, e, n = bbox + zone = int(((w + e) / 2 + 180) // 6) + 1 + epsg = 26900 + zone # NAD83 / UTM zone N + fwd = Transformer.from_crs("EPSG:4326", f"EPSG:{epsg}", always_xy=True) + xs, ys = fwd.transform([w, e, w, e], [s, s, n, n]) + cells = {} + for ee in range(int(min(xs) // 1000), int(max(xs) // 1000) + 1): + for nn in range(int(min(ys) // 1000), int(max(ys) // 1000) + 1): + cells[f"{ee:03d}{nn % 1000:03d}"] = (ee, nn) + return cells, zone, epsg + + +def _lidar_projects_covering(bbox, cfg: Config) -> list[tuple[str, str]]: + """The USGS 3DEP index → [(project, staged LAZ dir URL)] for every lidar project + whose extent intersects bbox, anywhere in the US.""" + w, s, e, n = bbox + data = with_retries(lambda: requests.get(cfg.index_lpc_url, params={ + "geometry": f"{w},{s},{e},{n}", "geometryType": "esriGeometryEnvelope", + "inSR": "4326", "spatialRel": "esriSpatialRelIntersects", + "outFields": "project,lpc_link", "returnGeometry": "false", "f": "json", + }, timeout=60).json()) + out = [] + for f in data.get("features", []): + a = f.get("attributes", {}) + link = (a.get("lpc_link") or "").rstrip("/") + if link: + out.append((a.get("project") or link.split("/")[-1], link)) + return out + + +def find_tiles_rockyweb(bbox, cfg: Config) -> list[dict]: + """Region-agnostic fallback when TNM is down. Discovers covering projects via the + 3DEP index, browses each staged LAZ dir, and selects tiles whose UTM-named grid + cell covers bbox. Non-UTM-named projects (state-plane, sequential ids) can't be + decoded by name and are logged + skipped (fetch those via TNM).""" + cells, zone, epsg = _covering_cells(bbox) + inv = Transformer.from_crs(f"EPSG:{epsg}", "EPSG:4326", always_xy=True) + tiles, seen = [], set() + for project, laz_dir in _lidar_projects_covering(bbox, cfg): + laz_url = f"{laz_dir}/LAZ/" + try: + listing = with_retries(lambda u=laz_url: requests.get(u, timeout=90).text) + except Exception as ex: # noqa: BLE001 + print(f" {project}: LAZ dir unreachable ({type(ex).__name__}) — skip.") + continue + matched = 0 + for fn_stem, ztile, code in {m for m in _TILE_RE.findall(listing)}: + if int(ztile) != zone or code not in cells or fn_stem in seen: + continue + seen.add(fn_stem) + matched += 1 + e_km, n_km = cells[code] + e0, n0 = e_km * 1000, n_km * 1000 + lons, lats = inv.transform([e0, e0 + 1000], [n0, n0 + 1000]) + tiles.append({ + "downloadURL": f"{laz_url}{fn_stem}.laz", + "title": f"USGS Lidar Point Cloud {project} {fn_stem.split('_')[-1]}", + "sizeInBytes": 0, + "boundingBox": {"minX": min(lons), "maxX": max(lons), + "minY": min(lats), "maxY": max(lats)}, + }) + if matched: + print(f" {project}: {matched} covering tiles (UTM zone {zone}).") + else: + why = ("non-UTM tile naming — skip (use TNM for this one)" + if ".laz" in listing else "no LAZ found") + print(f" {project}: 0 covering tiles ({why}).") + return tiles + + +# ── Best-quality collection mosaic ──────────────────────────────────────────── + +def _collection_of(tile: dict) -> str: + """Collection = the staged '/Projects//' segment of the download URL.""" + m = re.search(r"/Projects/([^/]+)/", tile.get("downloadURL", "")) + return m.group(1) if m else tile.get("title", "unknown") + + +def _tile_bbox(tile: dict): + bb = tile.get("boundingBox") or {} + if all(k in bb for k in ("minX", "minY", "maxX", "maxY")): + return (bb["minX"], bb["minY"], bb["maxX"], bb["maxY"]) + return None + + +def _collection_density(ctiles: list[dict]) -> float: + """Density proxy = median file bytes per km² (what actually drives RANSAC).""" + vals = [] + for t in ctiles: + bb, sz = _tile_bbox(t), (t.get("sizeInBytes") or 0) + if not bb or sz <= 0: + continue + x0, y0, x1, y1 = bb + lat = (y0 + y1) / 2 + area = abs((x1 - x0) * 111.32 * math.cos(math.radians(lat)) * (y1 - y0) * 110.54) + if area > 0: + vals.append(sz / area) + return float(np.median(vals)) if vals else 0.0 + + +def _rank_key(name: str, ctiles: list[dict]): + """Rank collections by density, then recency, then QL hint (all descending).""" + years = [int(m.group(1)) for t in ctiles + for m in [re.search(r"(20\d\d)", t.get("publicationDate", ""))] if m] + ql = 2 if "QL1" in name else (1 if "QL2" in name else 0) + return (_collection_density(ctiles), max(years) if years else 0, ql) + + +def select_collections(tiles: list[dict], bbox, cell_deg: float = 0.01) -> list[dict]: + """Choose a best-quality LiDAR mosaic: rank collections by density/recency, then + greedily assign each ~1 km cell to the highest-ranked collection covering it, + filling gaps with lower-ranked ones. Returns the tiles needed for that mosaic.""" + colls: dict = {} + for t in tiles: + colls.setdefault(_collection_of(t), []).append(t) + + if len(colls) == 1: + name = next(iter(colls)) + print(f" 1 collection covers this region: {name} ({len(tiles)} tiles).") + return tiles + + ranked = sorted(colls.items(), key=lambda kv: _rank_key(*kv), reverse=True) + print(f" {len(colls)} overlapping collections — ranking by density/recency:") + for name, ct in ranked: + print(f" {name}: ~{_collection_density(ct) / 1e3:.0f} KB/km², {len(ct)} tiles") + + w, s, e, n = bbox + xs, ys = np.arange(w, e, cell_deg), np.arange(s, n, cell_deg) + cells = [(x + cell_deg / 2, y + cell_deg / 2) for x in xs for y in ys] \ + or [((w + e) / 2, (s + n) / 2)] + covered = [False] * len(cells) + chosen, chosen_urls, usage = [], set(), {} + + for name, ct in ranked: + boxes = [(_tile_bbox(t), t) for t in ct] + for ci, (cx, cy) in enumerate(cells): + if covered[ci]: + continue + for bb, t in boxes: + if bb and bb[0] <= cx <= bb[2] and bb[1] <= cy <= bb[3]: + covered[ci] = True + usage[name] = usage.get(name, 0) + 1 + if t["downloadURL"] not in chosen_urls: + chosen_urls.add(t["downloadURL"]) + chosen.append(t) + break + + total = len(cells) + print(" Selected mosaic:") + for name, _ct in ranked: + if usage.get(name): + print(f" {name}: {100 * usage[name] / total:.0f}% of area") + dropped = [name for name, _ct in ranked if not usage.get(name)] + if dropped: + print(f" dropped (redundant / lower quality): {', '.join(dropped)}") + return chosen + + +def discover_tiles(bbox, cfg: Config) -> list[dict]: + """Full discovery for a region: TNM (with rockyweb fallback), then collection + selection (or a forced ``preferred_collection``), honoring ``max_tiles``.""" + tiles = with_retries(find_all_tiles, bbox, cfg) + if not tiles: + print(" TNM returned 0 (spatial query down) — falling back to rockyweb...") + tiles = with_retries(find_tiles_rockyweb, bbox, cfg) + print(f" rockyweb found {len(tiles)} tiles.") + if not tiles: + return [] + print(f" {len(tiles)} tiles cover this region.") + + if cfg.preferred_collection: + tiles = [t for t in tiles if cfg.preferred_collection in _collection_of(t)] + print(f" Forced collection '{cfg.preferred_collection}': {len(tiles)} tiles.") + else: + tiles = select_collections(tiles, bbox) + + if cfg.max_tiles is not None: + tiles = tiles[:cfg.max_tiles] + print(f" max_tiles set — processing only the first {len(tiles)} this run.") + return tiles + + +# ── Download + read ─────────────────────────────────────────────────────────── + +def download_tile(download_url: str, cache_dir: str) -> str: + """Download a tile to ``cache_dir`` (streamed), skipping if already present.""" + os.makedirs(cache_dir, exist_ok=True) + filename = download_url.split("/")[-1].split("?")[0] + local_path = os.path.join(cache_dir, filename) + if os.path.exists(local_path): + return local_path + + r = requests.get(download_url, stream=True, timeout=600) + r.raise_for_status() + with open(local_path, "wb") as f: + for chunk in r.iter_content(chunk_size=1024 * 1024): + f.write(chunk) + return local_path + + +def get_tile_srs(tile_path: str) -> str: + """Read the tile's spatial reference (WKT) from its LAS header.""" + import pdal # lazy: PDAL is a system/conda dependency + + p = pdal.Pipeline(json.dumps({"pipeline": [ + {"type": "readers.las", "filename": tile_path.replace("\\", "/"), "count": 1} + ]})) + p.execute() + md = p.metadata["metadata"]["readers.las"] + return md.get("comp_spatialreference") or md.get("spatialreference") + + +def load_tile_points(tile_path: str) -> TilePoints: + """Read the entire tile into memory once as (points, X, Y) float views. Reading + once and masking per-building in RAM is ~100× faster than re-opening the file.""" + import pdal # lazy + + p = pdal.Pipeline(json.dumps({"pipeline": [ + {"type": "readers.las", "filename": tile_path.replace("\\", "/")} + ]})) + p.execute() + pts = p.arrays[0] + return pts, np.asarray(pts["X"], dtype=float), np.asarray(pts["Y"], dtype=float) diff --git a/solargrader/storage.py b/solargrader/storage.py new file mode 100644 index 0000000..301e955 --- /dev/null +++ b/solargrader/storage.py @@ -0,0 +1,124 @@ +""" +Result storage — a thin class that OWNS the DuckDB connection. + +DuckDB is single-writer, so exactly one ``ResultStore`` (in the main process) holds +the connection; worker processes never touch the DB. The store manages the schema, +scored-home writes, and the resume ledgers (which tiles / sub-regions are done). +""" + +from __future__ import annotations + +import duckdb +import pandas as pd + +from .models import ScoredHome + +_HOMES_SCHEMA = """ + CREATE TABLE IF NOT EXISTS homes ( + osm_id BIGINT PRIMARY KEY, + lat DOUBLE, + lon DOUBLE, + footprint_m2 DOUBLE, + usable_roof_area_m2 DOUBLE, + primary_tilt_deg DOUBLE, + primary_azimuth_deg DOUBLE, + roof_plane_count INTEGER, + res_annual_kwh DOUBLE, + res_system_kw DOUBLE, + res_panel_count INTEGER, + solar_score INTEGER, + solar_grade VARCHAR(2), + max_annual_kwh DOUBLE, + max_system_kw DOUBLE, + max_panel_count INTEGER, + potential_grade VARCHAR(2), + shade_loss_pct DOUBLE, + processed_at TIMESTAMP + ) +""" + +_GRADE_ORDER = ["A+", "A", "B+", "B", "C", "D"] + + +class ResultStore: + """Owns the accumulating results DB. Use as a context manager: + + >>> with ResultStore("solar_grader.duckdb") as store: + ... store.save(scored_home) + """ + + def __init__(self, db_path: str, read_only: bool = False): + self.db_path = db_path + self.con = duckdb.connect(db_path, read_only=read_only) + if not read_only: + self._ensure_schema() + + def _ensure_schema(self) -> None: + self.con.execute(_HOMES_SCHEMA) + self.con.execute("CREATE TABLE IF NOT EXISTS tiles_done " + "(tile_id VARCHAR PRIMARY KEY, n_scored INTEGER, processed_at TIMESTAMP)") + self.con.execute("CREATE TABLE IF NOT EXISTS regions_done " + "(region_key VARCHAR PRIMARY KEY, processed_at TIMESTAMP)") + + # ── Writes ──────────────────────────────────────────────────────────────── + def save(self, home: ScoredHome) -> None: + """Insert (or replace) one scored home.""" + b = home.building + self.con.execute( + "INSERT OR REPLACE INTO homes VALUES " + "(?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)", + [b.id, b.lat, b.lon, b.footprint_m2, home.usable_area_m2, + home.primary_tilt_deg, home.primary_azimuth_deg, home.plane_count, + home.res_kwh, home.res_kw, home.res_panels, home.score, home.grade, + home.max_kwh, home.max_kw, home.max_panels, home.potential_grade, + home.shade_loss_pct, pd.Timestamp.now()], + ) + + # ── Resume ledgers ──────────────────────────────────────────────────────── + def mark_tile_done(self, tile_id: str, n_scored: int) -> None: + self.con.execute("INSERT OR REPLACE INTO tiles_done VALUES (?, ?, ?)", + [tile_id, n_scored, pd.Timestamp.now()]) + + def done_tiles(self) -> set[str]: + return {r[0] for r in self.con.execute("SELECT tile_id FROM tiles_done").fetchall()} + + def mark_region_done(self, region_key: str) -> None: + self.con.execute("INSERT OR REPLACE INTO regions_done VALUES (?, ?)", + [region_key, pd.Timestamp.now()]) + + def done_regions(self) -> set[str]: + return {r[0] for r in self.con.execute("SELECT region_key FROM regions_done").fetchall()} + + # ── Reporting ───────────────────────────────────────────────────────────── + def count_homes(self) -> int: + return self.con.execute("SELECT COUNT(*) FROM homes").fetchone()[0] + + def grade_breakdown(self) -> dict: + rows = self.con.execute("SELECT solar_grade, COUNT(*) FROM homes GROUP BY 1").fetchall() + return dict(rows) + + def print_summary(self, elapsed: float | None = None) -> None: + n_homes = self.count_homes() + n_tiles = self.con.execute("SELECT COUNT(*) FROM tiles_done").fetchone()[0] + print(f"\n{'=' * 60}\nPIPELINE SUMMARY\n{'=' * 60}") + print(f" Tiles processed: {n_tiles}") + print(f" Homes scored: {n_homes}") + if elapsed is not None: + print(f" Total time: {elapsed:.1f}s") + counts = self.grade_breakdown() + print("\n Grade breakdown:") + for g in _GRADE_ORDER: + n = counts.get(g, 0) + if n: + print(f" {g:2s} {'█' * min(n, 50)} ({n})") + print(f"\n Results saved to: {self.db_path}") + + # ── Lifecycle ───────────────────────────────────────────────────────────── + def close(self) -> None: + self.con.close() + + def __enter__(self) -> ResultStore: + return self + + def __exit__(self, *exc) -> None: + self.close() diff --git a/solargrader/util.py b/solargrader/util.py new file mode 100644 index 0000000..91440f3 --- /dev/null +++ b/solargrader/util.py @@ -0,0 +1,29 @@ +"""Small shared utilities with no domain logic.""" + +from __future__ import annotations + +import time +from typing import Callable, TypeVar + +T = TypeVar("T") + + +def with_retries(fn: Callable[..., T], *args, + attempts: int = 4, base_delay: float = 5.0, **kwargs) -> T: + """Call ``fn`` with exponential backoff on transient failures (e.g. an Overpass + 504 or a PVGIS hiccup). Re-raises only after the final attempt. + + >>> with_retries(lambda: 1 + 1) + 2 + """ + for i in range(attempts): + try: + return fn(*args, **kwargs) + except Exception as e: # noqa: BLE001 — deliberate: retry any transient error + if i == attempts - 1: + raise + wait = base_delay * (2 ** i) + print(f" {type(e).__name__}: {e}\n" + f" retry {i + 1}/{attempts - 1} in {wait:.0f}s...") + time.sleep(wait) + raise AssertionError("unreachable") # pragma: no cover diff --git a/solargrader/webmap.py b/solargrader/webmap.py new file mode 100644 index 0000000..33c3672 --- /dev/null +++ b/solargrader/webmap.py @@ -0,0 +1,125 @@ +""" +Leaflet map builder — render scored + enriched homes to a self-contained ``map.html`` +(data embedded, no server needed). Pins are colored by solar grade; click for +address, grade, production, and potential. Defaults to residential leads only. +""" + +from __future__ import annotations + +import json +import warnings + +import duckdb + +from .config import Config + +warnings.filterwarnings("ignore") + +GRADES_SHOWN = ("A+", "A", "B+", "B", "C", "D") +GRADE_COLORS = { + "A+": "#1a9850", "A": "#66bd63", "B+": "#a6d96a", + "B": "#fee08b", "C": "#fdae61", "D": "#d73027", +} + + +def _load_features(cfg: Config, include_non_residential: bool): + con = duckdb.connect(cfg.db_path, read_only=True) + cols = {c[0] for c in con.execute("DESCRIBE homes").fetchall()} + if "is_residential" not in cols: + con.close() + raise SystemExit("No address/residential columns — run enrichment first.") + + where = ["solar_grade IN ({})".format(",".join(f"'{g}'" for g in GRADES_SHOWN))] + if not include_non_residential: + where.append("is_residential") + rows = con.execute(f""" + SELECT lat, lon, solar_grade, solar_score, + ROUND(res_annual_kwh) AS res_kwh, ROUND(res_system_kw, 1) AS res_kw, + potential_grade, ROUND(max_system_kw, 1) AS max_kw, + ROUND(shade_loss_pct) AS shade, + COALESCE(full_address, '(no address)') AS addr, COALESCE(city, '') AS city + FROM homes WHERE {' AND '.join(where)} + """).fetchdf() + con.close() + + features = [] + for r in rows.itertuples(index=False): + features.append({ + "type": "Feature", + "geometry": {"type": "Point", "coordinates": [r.lon, r.lat]}, + "properties": { + "grade": r.solar_grade, "score": int(r.solar_score), + "res_kwh": None if r.res_kwh != r.res_kwh else int(r.res_kwh), + "res_kw": r.res_kw, "pot": r.potential_grade, "max_kw": r.max_kw, + "shade": None if r.shade != r.shade else int(r.shade), + "addr": r.addr, "city": r.city, + }, + }) + return rows, features + + +def _build_html(cfg: Config, rows, features) -> str: + center_lat, center_lon = float(rows["lat"].mean()), float(rows["lon"].mean()) + geojson = json.dumps({"type": "FeatureCollection", "features": features}) + colors = json.dumps(GRADE_COLORS) + legend = "".join(f'
{g}
' + for g, c in GRADE_COLORS.items()) + + return f""" + +Solar Grader — Leads + + + + +
+""" + + +def build_map(cfg: Config, out_html: str = "map.html", + include_non_residential: bool = False) -> None: + rows, features = _load_features(cfg, include_non_residential) + if not features: + print("No homes to map — run the pipeline and enrichment first.") + return + with open(out_html, "w") as f: + f.write(_build_html(cfg, rows, features)) + kind = "all" if include_non_residential else "residential" + print(f"Wrote {out_html} with {len(features)} {kind} leads. Open it with: open {out_html}") diff --git a/test_pipeline.py b/test_pipeline.py deleted file mode 100644 index 7f0a94b..0000000 --- a/test_pipeline.py +++ /dev/null @@ -1,827 +0,0 @@ -""" -Solar Grader — Single Tile Pipeline Test -Tests the full pipeline on one LiDAR tile from Chester County, PA. - -Run with: - conda activate solar-grader - python test_pipeline.py -""" - -import json -import os -import time -import warnings - -import duckdb -import numpy as np -import open3d as o3d -import pandas as pd -import pdal -import pvlib -import requests -import shapely -from pyproj import Transformer -from scipy.spatial import ConvexHull -from shapely.geometry import Polygon -from shapely.ops import transform as shapely_transform - -warnings.filterwarnings("ignore") - -# ── Configuration ───────────────────────────────────────────────────────────── - -# Small test area: residential single-family streets, West Chester borough, PA. -# Picked because the goal is scoring HOMES — the old Exton bbox (below) was all -# commercial/retail. Adjust these 4 numbers to target a different area; keep it -# small (0.02° × 0.02° max). (west, south, east, north) -TEST_BBOX = (-75.615, 39.955, -75.600, 39.970) -TEST_LAT = 39.962 -TEST_LON = -75.607 - -# Previous test area (Exton commercial strip — big-box roofs, no houses): -# TEST_BBOX = (-75.640, 40.015, -75.620, 40.030); TEST_LAT/LON = 40.022, -75.630 - -TILE_CACHE_DIR = "tile_cache" # LiDAR tiles cached here (re-used on re-runs) -DB_PATH = "solar_results.duckdb" - -# ── Panel / system model ────────────────────────────────────────────────────── -# All physical: change these to model different panels/configurations. The energy -# model packs WHOLE panels onto each roof plane after a fire-code edge setback, -# rather than assuming a fixed % of the roof is usable. -PANEL_WIDTH_M = 1.05 # standard residential module footprint (~1.05 × 1.74 m) -PANEL_HEIGHT_M = 1.74 -PANEL_WATTS = 400 # module rated DC power at STC (1000 W/m²) -ROOF_SETBACK_M = 0.46 # fire-code clear pathway from plane edges (~18 inches) -PANEL_PACKING = 0.90 # fraction of the setback area that packs into whole modules -SYSTEM_LOSSES = 0.86 # inverter + wiring + soiling (~14% total loss) -MAX_AZIMUTH_OFFSET = 120 # skip planes facing >this many degrees from due south -MAX_SYSTEM_KW = 10.0 # cap for a typical sellable residential install (DC kW) -SHADE_RADIUS_M = 50.0 # search radius for shading obstructions (trees, neighbours) -HORIZON_BINS = 72 # azimuth bins for the roof skyline (5° each) - -PANEL_AREA_M2 = PANEL_WIDTH_M * PANEL_HEIGHT_M # m² per module -MODULE_EFF = PANEL_WATTS / (PANEL_AREA_M2 * 1000.0) # derived STC efficiency - -MIN_ROOF_POINTS = 20 # skip buildings with fewer LiDAR points than this -RANSAC_ITERATIONS = 300 # lower = faster, higher = more accurate -RANSAC_DISTANCE_THRESHOLD = 0.20 # meters; how far a point can be from a plane - - -# ── Step 1: Find & download a LiDAR tile ───────────────────────────────────── - -def find_lidar_tile(bbox): - """Query USGS TNM API for a LiDAR tile covering the bounding box.""" - west, south, east, north = bbox - url = "https://tnmaccess.nationalmap.gov/api/v1/products" - params = { - "datasets": "Lidar Point Cloud (LPC)", - "bbox": f"{west},{south},{east},{north}", - "prodFormats": "LAZ", - "max": 5, - } - print(f" Querying USGS TNM API...") - r = requests.get(url, params=params, timeout=30) - r.raise_for_status() - - items = r.json().get("items", []) - if not items: - print(" No tiles found — widening search area...") - expand = 0.05 - params["bbox"] = f"{west-expand},{south-expand},{east+expand},{north+expand}" - r = requests.get(url, params=params, timeout=30) - items = r.json().get("items", []) - - if not items: - raise RuntimeError( - "No LiDAR tiles found. Check https://apps.nationalmap.gov/downloader/ " - "manually for PA coverage." - ) - - tile = items[0] - size_mb = tile.get("sizeInBytes", 0) / 1e6 - print(f" Found: {tile['title']}") - print(f" Size: {size_mb:.0f} MB") - if size_mb > 300: - print(" WARNING: large tile — download may take a few minutes.") - return tile - - -def download_tile(download_url, cache_dir): - """Download a tile to cache_dir, skipping if already present.""" - os.makedirs(cache_dir, exist_ok=True) - filename = download_url.split("/")[-1].split("?")[0] - local_path = os.path.join(cache_dir, filename) - - if os.path.exists(local_path): - print(f" Using cached tile: {local_path}") - return local_path - - print(f" Downloading to {local_path} ...") - r = requests.get(download_url, stream=True, timeout=600) - r.raise_for_status() - - total = int(r.headers.get("content-length", 0)) - downloaded = 0 - with open(local_path, "wb") as f: - for chunk in r.iter_content(chunk_size=1024 * 1024): - f.write(chunk) - downloaded += len(chunk) - if total: - print( - f"\r {downloaded / 1e6:.1f} / {total / 1e6:.1f} MB " - f"({downloaded / total * 100:.0f}%)", - end="", - ) - print(f"\n Done.") - return local_path - - -# ── Step 2: Get building footprints from OpenStreetMap ──────────────────────── - -def get_buildings_from_osm(bbox): - """ - Fetch building polygons from Overpass API. - Returns a list of dicts with geometry and metadata. - """ - west, south, east, north = bbox - # Overpass bbox order is: south, west, north, east - query = f""" -[out:json][timeout:60]; -( - way["building"]({south},{west},{north},{east}); -); -out geom; -""" - print(" Querying Overpass API for buildings...") - r = requests.post( - "https://overpass-api.de/api/interpreter", - data={"data": query}, - headers={"User-Agent": "SolarGrader/1.0 (joshuawagnersd@gmail.com)"}, - timeout=90, - ) - r.raise_for_status() - - buildings = [] - for el in r.json().get("elements", []): - if el["type"] != "way" or "geometry" not in el: - continue - coords = [(n["lon"], n["lat"]) for n in el["geometry"]] - if len(coords) < 4: - continue - try: - poly = Polygon(coords) - if not poly.is_valid or poly.area == 0: - continue - cx, cy = poly.centroid.x, poly.centroid.y - # Rough m² from degrees: 1° lat ≈ 111,320m; 1° lon ≈ 111,320 × cos(lat) - deg2m2 = (111320 ** 2) * np.cos(np.radians(cy)) - buildings.append( - { - "osm_id": el["id"], - "geometry": poly, - "lat": cy, - "lon": cx, - "footprint_m2": poly.area * deg2m2, - } - ) - except Exception: - pass - - print(f" Found {len(buildings)} buildings.") - return buildings - - -# ── Step 3: Get solar irradiance data (PVGIS, free, no API key) ─────────────── - -def get_tmy_data(lat, lon): - """ - Fetch a Typical Meteorological Year from the EU PVGIS service. - Free, no API key required, covers the US. - Returns a DataFrame with hourly irradiance columns: ghi, dni, dhi. - """ - print(f" Fetching TMY data from PVGIS for ({lat:.3f}, {lon:.3f})...") - tmy, _ = pvlib.iotools.get_pvgis_tmy( # pvlib >=0.11 returns (data, metadata) - latitude=lat, - longitude=lon, - outputformat="json", - usehorizon=True, - map_variables=True, # renames columns to pvlib standard: ghi, dni, dhi - ) - print(f" TMY data loaded: {len(tmy)} hourly records.") - return tmy - - -# ── Step 4: LiDAR clipping ──────────────────────────────────────────────────── - -def get_tile_srs(tile_path): - """Read the LiDAR tile's spatial reference (WKT) from its LAS header.""" - p = pdal.Pipeline(json.dumps( - {"pipeline": [{"type": "readers.las", "filename": tile_path.replace("\\", "/"), "count": 1}]} - )) - p.execute() - md = p.metadata["metadata"]["readers.las"] - return md.get("comp_spatialreference") or md.get("spatialreference") - - -def load_tile_points(tile_path): - """ - Read the entire LiDAR tile into memory ONCE as a structured numpy array. - Returns (points, X, Y) where X/Y are float views for fast bbox masking. - Reading once and filtering per-building in RAM is ~100x faster than - re-opening the 145 MB file for every building. - """ - p = pdal.Pipeline(json.dumps( - {"pipeline": [{"type": "readers.las", "filename": tile_path.replace("\\", "/")}]} - )) - p.execute() - pts = p.arrays[0] - return pts, np.asarray(pts["X"], dtype=float), np.asarray(pts["Y"], dtype=float) - - -def clip_lidar_to_building(tile_points, footprint_poly, transformer, buffer_m=0.5): - """ - Extract LiDAR points inside a building footprint from already-loaded tile - points. The footprint comes in lon/lat (EPSG:4326); the tile is in projected - meters, so we reproject the polygon into the tile's CRS before clipping. - Returns an (N, 3) numpy array of [X, Y, Z], or None if too few points. - """ - pts, X, Y = tile_points - proj_poly = shapely_transform(transformer.transform, footprint_poly).buffer(buffer_m) - - # Fast bounding-box prefilter, then exact point-in-polygon on the small subset - minx, miny, maxx, maxy = proj_poly.bounds - bbox_mask = (X >= minx) & (X <= maxx) & (Y >= miny) & (Y <= maxy) - if not bbox_mask.any(): - return None - - sub = pts[bbox_mask] - inside = shapely.contains_xy(proj_poly, sub["X"], sub["Y"]) - sub = sub[inside] - if len(sub) == 0: - return None - - # Prefer building-classified points (class 6); fall back to height filter - if "Classification" in sub.dtype.names: - roof_pts = sub[sub["Classification"] == 6] - if len(roof_pts) < MIN_ROOF_POINTS: - z_min = sub["Z"].min() - roof_pts = sub[sub["Z"] > z_min + 2.0] - else: - z_min = sub["Z"].min() - roof_pts = sub[sub["Z"] > z_min + 2.0] - - if len(roof_pts) < MIN_ROOF_POINTS: - return None - - return np.column_stack([roof_pts["X"], roof_pts["Y"], roof_pts["Z"]]) - - -# ── Step 5: Roof plane extraction (RANSAC) ──────────────────────────────────── - -def extract_roof_planes(points_xyz): - """ - Fit roof planes to LiDAR points using Open3D RANSAC. - Returns a list of plane dicts with tilt, azimuth, and area. - """ - pcd = o3d.geometry.PointCloud() - pcd.points = o3d.utility.Vector3dVector(points_xyz) - - # Ground point density (points per m² of horizontal area) from ONE convex - # hull over all points. Used to size each plane by its point count — this - # avoids the inflation you get from summing overlapping per-plane hulls. - # Coordinates are in projected meters (tile CRS), so areas are already m². - all_xy = np.asarray(points_xyz)[:, :2] - try: - covered_m2 = float(ConvexHull(all_xy).volume) if len(all_xy) >= 3 else 0.0 - except Exception: - covered_m2 = 0.0 - ground_density = (len(points_xyz) / covered_m2) if covered_m2 > 0 else 0.0 - - planes = [] - remaining = pcd - - for _ in range(6): # up to 6 roof faces - if len(remaining.points) < MIN_ROOF_POINTS: - break - - plane_model, inliers = remaining.segment_plane( - distance_threshold=RANSAC_DISTANCE_THRESHOLD, - ransac_n=3, - num_iterations=RANSAC_ITERATIONS, - ) - - if len(inliers) < MIN_ROOF_POINTS: - break - - # Plane equation: ax + by + cz + d = 0 → normal = [a, b, c] - normal = np.array(plane_model[:3]) - if normal[2] < 0: - normal = -normal # ensure normal points upward - - norm_len = np.linalg.norm(normal) - if norm_len == 0: - remaining = remaining.select_by_index(inliers, invert=True) - continue - - nx, ny, nz = normal / norm_len - - # Tilt from horizontal: 0° = flat, 90° = vertical wall - tilt_deg = float(np.degrees(np.arccos(min(nz, 1.0)))) - - # Skip near-vertical surfaces (walls, dormers) - if tilt_deg > 72: - remaining = remaining.select_by_index(inliers, invert=True) - continue - - # Azimuth: 0°=N, 90°=E, 180°=S, 270°=W (compass convention) - # In projected coords, +X=east, +Y=north → atan2(nx, ny) gives bearing - azimuth_deg = float((np.degrees(np.arctan2(nx, ny)) + 360) % 360) - - # Plane area = (inlier point count ÷ ground density) gives the plane's - # horizontal (ground-projected) area; dividing by cos(tilt) converts it - # to true tilted surface area. Each point belongs to exactly one plane, - # so summing these areas can never exceed the building's real footprint. - if ground_density > 0: - ground_area_m2 = len(inliers) / ground_density - area_m2 = ground_area_m2 / max(np.cos(np.radians(tilt_deg)), 0.1) - else: - area_m2 = 0.0 - - planes.append( - { - "tilt_deg": tilt_deg, - "azimuth_deg": azimuth_deg, - "area_m2": area_m2, - "point_count": len(inliers), - } - ) - - remaining = remaining.select_by_index(inliers, invert=True) - - return planes - - -# ── Step 6: Solar potential calculation (pvlib) ─────────────────────────────── - -def panels_on_plane(plane_area_m2): - """ - How many whole panels physically fit on a roof plane of this area. - Approximate the plane as an equal-area square, erode it inward by the - fire-code setback on all sides, apply a packing factor, then divide by the - module footprint. Small/cut-up roofs lose proportionally more area — the - realistic behaviour a flat usable-% can't capture. - """ - side = np.sqrt(plane_area_m2) - usable_side = max(0.0, side - 2 * ROOF_SETBACK_M) - usable_area = (usable_side ** 2) * PANEL_PACKING - return int(usable_area // PANEL_AREA_M2) - - -def compute_horizon(tile_points, observer_xyz, exclude_poly=None, - radius_m=SHADE_RADIUS_M, n_bins=HORIZON_BINS): - """ - Build the roof's skyline: the maximum obstruction elevation (degrees above the - roof) in each compass direction, using surrounding LiDAR points — trees, - neighbouring buildings — within radius_m. Points inside exclude_poly (the - building's OWN footprint, in tile CRS) are dropped so a roof can't shade - itself. Returns an array indexed by azimuth bin (0=N, clockwise); all zeros = - wide-open sky. - """ - pts, X, Y = tile_points - ox, oy, oz = observer_xyz - horizon = np.zeros(n_bins) - - m = (X >= ox - radius_m) & (X <= ox + radius_m) & (Y >= oy - radius_m) & (Y <= oy + radius_m) - if not m.any(): - return horizon - - sx, sy = X[m], Y[m] - sz = np.asarray(pts["Z"], dtype=float)[m] - # Drop the building's own roof so it can't count as an obstruction to itself - if exclude_poly is not None: - outside = ~shapely.contains_xy(exclude_poly, sx, sy) - sx, sy, sz = sx[outside], sy[outside], sz[outside] - - dx = sx - ox - dy = sy - oy - dz = sz - oz - horiz = np.sqrt(dx * dx + dy * dy) - # Keep points within radius, not right on top of us, and ABOVE the roof - keep = (horiz <= radius_m) & (horiz > 1.0) & (dz > 0) - if not keep.any(): - return horizon - - dx, dy, dz, horiz = dx[keep], dy[keep], dz[keep], horiz[keep] - elev = np.degrees(np.arctan2(dz, horiz)) # elevation angle - az = (np.degrees(np.arctan2(dx, dy)) + 360) % 360 # 0=N, 90=E (compass) - bins = (az / (360.0 / n_bins)).astype(int) % n_bins - np.maximum.at(horizon, bins, elev) # tallest per direction - return horizon - - -def calculate_annual_kwh(lat, lon, tmy, planes, horizon=None): - """ - Pack whole panels onto each viable roof plane (pvlib irradiance per plane), - then report TWO systems: - - max: every panel that fits the roof (upside potential) - - residential: the best-producing panels up to MAX_SYSTEM_KW (sellable now) - If a horizon (skyline) is given, the direct beam is removed during hours when - the sun sits below the skyline, and shade_loss_pct reports the production lost. - Returns a dict (incl. shade_loss_pct) or None if no viable planes exist. - """ - location = pvlib.location.Location(latitude=lat, longitude=lon, tz="Etc/GMT+5") - solar_pos = location.get_solarposition(tmy.index) - - # Per-hour direct-beam shading mask: sun up but below the roof's skyline - sun_blocked = None - if horizon is not None: - sun_elev = 90.0 - solar_pos["apparent_zenith"].values - sun_az = solar_pos["azimuth"].values - sun_bins = (sun_az / (360.0 / len(horizon))).astype(int) % len(horizon) - sun_blocked = (sun_elev > 0) & (sun_elev <= horizon[sun_bins]) - - viable = [] # (per_panel_shaded, per_panel_unshaded, n_panels, plane) - for plane in planes: - # Skip planes facing too far from due south to be worth paneling - if abs(plane["azimuth_deg"] - 180) > MAX_AZIMUTH_OFFSET: - continue - - n_panels = panels_on_plane(plane["area_m2"]) - if n_panels <= 0: - continue - - poa = pvlib.irradiance.get_total_irradiance( - surface_tilt=plane["tilt_deg"], - surface_azimuth=plane["azimuth_deg"], - solar_zenith=solar_pos["apparent_zenith"], - solar_azimuth=solar_pos["azimuth"], - ghi=tmy["ghi"], - dni=tmy["dni"], - dhi=tmy["dhi"], - ) - glob = poa["poa_global"].fillna(0).values - unshaded_annual = glob.sum() / 1000 # kWh/m²/yr with full sun - if sun_blocked is not None: - # Remove the beam component during blocked hours (diffuse still counts) - direct = poa["poa_direct"].fillna(0).values - shaded_annual = (glob - direct * sun_blocked).sum() / 1000 - else: - shaded_annual = unshaded_annual - - ppk_shaded = shaded_annual * PANEL_AREA_M2 * MODULE_EFF * SYSTEM_LOSSES - ppk_unshaded = unshaded_annual * PANEL_AREA_M2 * MODULE_EFF * SYSTEM_LOSSES - viable.append((ppk_shaded, ppk_unshaded, n_panels, plane)) - - if not viable: - return None - - # Max system: every panel that physically fits the roof - max_panels = sum(n for _, _, n, _ in viable) - max_kwh = sum(s * n for s, _, n, _ in viable) - max_unshaded = sum(u * n for _, u, n, _ in viable) - max_kw = max_panels * PANEL_WATTS / 1000.0 - shade_loss_pct = 100.0 * (1 - max_kwh / max_unshaded) if max_unshaded > 0 else 0.0 - - # Residential system: fill the best-producing panels first, up to the kW cap - cap_panels = int(MAX_SYSTEM_KW * 1000 / PANEL_WATTS) - res_panels = 0 - res_kwh = 0.0 - for ppk_shaded, _, n_panels, _ in sorted(viable, key=lambda v: v[0], reverse=True): - take = min(n_panels, cap_panels - res_panels) - if take <= 0: - break - res_panels += take - res_kwh += take * ppk_shaded - res_kw = res_panels * PANEL_WATTS / 1000.0 - - best_plane = max(viable, key=lambda v: v[0])[3] - return { - "res_kwh": res_kwh, "res_kw": res_kw, "res_panels": res_panels, - "max_kwh": max_kwh, "max_kw": max_kw, "max_panels": max_panels, - "shade_loss_pct": shade_loss_pct, "best_plane": best_plane, - } - - -# ── Step 7: Grading ─────────────────────────────────────────────────────────── - -def grade_home(res_annual_kwh, res_system_kw, primary_azimuth_deg): - """ - Residential SALES grade (0-100) — how good a lead this roof is for a standard - home system. Driven by SPECIFIC YIELD (kWh per kW installed): a single physical - number that already integrates orientation, tilt, and shading. Total energy is - deliberately NOT the main axis — the residential system is capped (MAX_SYSTEM_KW), - so ~all adequately-sized roofs saturate that cap and total kWh stopped - discriminating (it collapsed every home to A/A+). Yield restores the spread: - a south-facing unshaded PA roof lands ~1400-1500 kWh/kW; a poorly-oriented or - shaded one ~900-1000. Roof SIZE is intentionally the other grade's job - (grade_potential / P1-P5), so it enters here only as a small-system penalty, - not a reward the capped majority would all collect. - """ - yield_kwh_per_kw = res_annual_kwh / res_system_kw if res_system_kw else 0.0 - score = 0 - - # Quality (dominant): specific yield. PA rooftop range ~900-1500 kWh/kW. - if yield_kwh_per_kw >= 1400: score += 65 - elif yield_kwh_per_kw >= 1300: score += 52 - elif yield_kwh_per_kw >= 1200: score += 38 - elif yield_kwh_per_kw >= 1100: score += 24 - elif yield_kwh_per_kw >= 1000: score += 12 - else: score += 4 - - # Orientation: reinforces prime south-facing among similar-yield roofs. - offset = abs(primary_azimuth_deg - 180) - if offset <= 15: score += 25 - elif offset <= 45: score += 17 - elif offset <= 75: score += 7 - - # Deal size: don't reward the capped majority; only penalize sub-residential roofs. - if res_system_kw >= 7: score += 10 - elif res_system_kw >= 5: score += 5 - elif res_system_kw < 3: score -= 10 - - if score >= 85: return "A+", score - elif score >= 70: return "A", score - elif score >= 55: return "B+", score - elif score >= 40: return "B", score - elif score >= 25: return "C", score - else: return "D", score - - -def grade_potential(max_kw): - """ - Tier a roof by its MAXIMUM installable capacity (kW DC) — the upside beyond a - standard residential system (battery, EV, future expansion, or commercial). - Separate from the residential sales grade, which is about the sellable system. - """ - if max_kw >= 30: return "P1" # very high — large/commercial-scale roof - elif max_kw >= 20: return "P2" - elif max_kw >= 12: return "P3" - elif max_kw >= 7: return "P4" - else: return "P5" # limited — barely fits a full residential set - - -# ── Step 8: Database ────────────────────────────────────────────────────────── - -def setup_db(path): - con = duckdb.connect(path) - # Fresh table each run so test results reflect only the current tile/bbox - # (the real multi-tile pipeline would accumulate instead of dropping). - con.execute("DROP TABLE IF EXISTS homes") - con.execute(""" - CREATE TABLE IF NOT EXISTS homes ( - osm_id BIGINT PRIMARY KEY, - lat DOUBLE, - lon DOUBLE, - footprint_m2 DOUBLE, - usable_roof_area_m2 DOUBLE, - primary_tilt_deg DOUBLE, - primary_azimuth_deg DOUBLE, - roof_plane_count INTEGER, - res_annual_kwh DOUBLE, - res_system_kw DOUBLE, - res_panel_count INTEGER, - solar_score INTEGER, - solar_grade VARCHAR(2), - max_annual_kwh DOUBLE, - max_system_kw DOUBLE, - max_panel_count INTEGER, - potential_grade VARCHAR(2), - shade_loss_pct DOUBLE, - processed_at TIMESTAMP - ) - """) - return con - - -def save_result(con, building, result): - con.execute( - """ - INSERT OR REPLACE INTO homes VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?) - """, - [ - building["osm_id"], - building["lat"], - building["lon"], - building["footprint_m2"], - result["usable_area_m2"], - result["primary_tilt"], - result["primary_azimuth"], - result["plane_count"], - result["res_kwh"], - result["res_kw"], - result["res_panels"], - result["score"], - result["grade"], - result["max_kwh"], - result["max_kw"], - result["max_panels"], - result["potential_grade"], - result["shade_loss_pct"], - pd.Timestamp.now(), - ], - ) - - -# ── Main ────────────────────────────────────────────────────────────────────── - -def main(): - print("=" * 60) - print("SOLAR GRADER — Pipeline Test") - print(f"Test area bbox: {TEST_BBOX}") - print("=" * 60) - - # 1. Find tile - print("\n[1/5] Finding LiDAR tile...") - tile_info = find_lidar_tile(TEST_BBOX) - - # 2. Download tile - print("\n[2/5] Downloading tile...") - tile_path = download_tile(tile_info["downloadURL"], TILE_CACHE_DIR) - - # Build a transformer from lon/lat (OSM) into the tile's projected CRS (meters) - tile_srs = get_tile_srs(tile_path) - transformer = Transformer.from_crs("EPSG:4326", tile_srs, always_xy=True) - print(f" Tile CRS: {tile_srs.split(',')[0].split('[')[-1].strip(chr(34))}") - - # Load the whole tile into memory once (filtered per-building in RAM below) - print(" Loading tile points into memory...") - tile_points = load_tile_points(tile_path) - print(f" Loaded {len(tile_points[0]):,} points.") - - # 3. Get buildings - print("\n[3/5] Fetching building footprints from OpenStreetMap...") - buildings = get_buildings_from_osm(TEST_BBOX) - if not buildings: - print("ERROR: No buildings found. Try adjusting TEST_BBOX.") - return - - # 4. Solar irradiance data (one fetch covers all buildings in the area) - print("\n[4/5] Fetching solar irradiance data...") - tmy = get_tmy_data(TEST_LAT, TEST_LON) - - # 5. Process each building - print(f"\n[5/5] Processing {len(buildings)} buildings...") - con = setup_db(DB_PATH) - - counts = {"scored": 0, "skipped_no_points": 0, "skipped_no_planes": 0} - grade_counts = {} - timings = [] - density_samples = [] # (n_points, footprint_m2) for buildings that had LiDAR - t_total = time.time() - - for i, building in enumerate(buildings): - t0 = time.time() - print( - f"\r [{i+1}/{len(buildings)}] scored={counts['scored']} " - f"skipped={counts['skipped_no_points']+counts['skipped_no_planes']} ", - end="", - ) - - # Clip LiDAR - points = clip_lidar_to_building(tile_points, building["geometry"], transformer) - if points is not None: - density_samples.append((len(points), building["footprint_m2"])) - if points is None: - counts["skipped_no_points"] += 1 - continue - - # Extract roof planes - planes = extract_roof_planes(points) - if not planes: - counts["skipped_no_planes"] += 1 - continue - - # Skyline shading from surrounding trees/buildings (observer = roof centroid), - # excluding the building's own footprint so it can't shade itself - observer = points.mean(axis=0) - own_footprint = shapely_transform(transformer.transform, building["geometry"]).buffer(1.0) - horizon = compute_horizon(tile_points, observer, own_footprint) - - # Solar calc — returns both a residential-capped and a max-roof system - sysres = calculate_annual_kwh( - building["lat"], building["lon"], tmy, planes, horizon - ) - if sysres is None or sysres["res_panels"] == 0: - counts["skipped_no_planes"] += 1 - continue - - # Grade: residential = sellable system (lead quality); potential = max roof - best_plane = sysres["best_plane"] - total_area = sum(p["area_m2"] for p in planes) - grade, score = grade_home( - sysres["res_kwh"], sysres["res_kw"], best_plane["azimuth_deg"] - ) - potential_grade = grade_potential(sysres["max_kw"]) - - save_result( - con, - building, - { - "usable_area_m2": total_area, - "primary_tilt": best_plane["tilt_deg"], - "primary_azimuth": best_plane["azimuth_deg"], - "plane_count": len(planes), - "res_kwh": sysres["res_kwh"], - "res_kw": sysres["res_kw"], - "res_panels": sysres["res_panels"], - "max_kwh": sysres["max_kwh"], - "max_kw": sysres["max_kw"], - "max_panels": sysres["max_panels"], - "potential_grade": potential_grade, - "shade_loss_pct": sysres["shade_loss_pct"], - "grade": grade, - "score": score, - }, - ) - - counts["scored"] += 1 - grade_counts[grade] = grade_counts.get(grade, 0) + 1 - timings.append(time.time() - t0) - - elapsed = time.time() - t_total - - # ── Results ─────────────────────────────────────────────────────────────── - print(f"\n\n{'=' * 60}") - print("RESULTS") - print(f"{'=' * 60}") - print(f" Buildings found: {len(buildings)}") - print(f" Successfully scored: {counts['scored']}") - print(f" Skipped (no LiDAR): {counts['skipped_no_points']}") - print(f" Skipped (no roof found): {counts['skipped_no_planes']}") - if timings: - print(f" Avg time per building: {np.mean(timings):.2f}s") - print(f" Total time: {elapsed:.1f}s") - - print(f"\n Grade breakdown:") - for g in ["A+", "A", "B+", "B", "C", "D"]: - n = grade_counts.get(g, 0) - if n: - print(f" {g:2s} {'█' * n} ({n})") - - print(f"\n Top 10 A/A+ leads (res = sellable ~{MAX_SYSTEM_KW:.0f}kW system, max = full-roof upside):") - df = con.execute(""" - SELECT lat, lon, - solar_grade AS grade, - ROUND(res_system_kw, 1) AS res_kw, - ROUND(res_annual_kwh) AS res_kwh, - ROUND(shade_loss_pct) AS shade_pct, - potential_grade AS pot, - ROUND(max_system_kw, 1) AS max_kw, - ROUND(max_annual_kwh) AS max_kwh, - ROUND(primary_tilt_deg) AS tilt, - ROUND(primary_azimuth_deg) AS azimuth - FROM homes - WHERE solar_grade IN ('A+', 'A') - ORDER BY res_annual_kwh DESC, max_annual_kwh DESC - LIMIT 10 - """).fetchdf() - - if df.empty: - print(" No A/A+ homes found in this tile.") - print(" Try a denser suburban area, or check LiDAR point density.") - else: - print(df.to_string(index=False)) - - # Realism check: validate the PHYSICS on an UNSHADED basis (~1,100–1,300 in PA), - # and report real-world shading separately so it doesn't look like a model error. - print(f"\n Realism check (UNSHADED kWh per kW — PA expected ~1,100–1,300):") - yields = con.execute(""" - SELECT (res_annual_kwh / NULLIF(1 - shade_loss_pct / 100.0, 0)) / res_system_kw - AS unshaded_yield, - shade_loss_pct - FROM homes WHERE res_system_kw > 0 - """).fetchdf() - if not yields.empty: - med_yield = float(yields["unshaded_yield"].median()) - med_shade = float(yields["shade_loss_pct"].median()) - print(f" Median unshaded yield: ~{med_yield:,.0f} kWh/kW") - if 950 <= med_yield <= 1450: - print(" OK — the panel/energy physics is in the realistic range.") - else: - print(" WARNING — outside the expected band; the panel/energy model needs review.") - print(f" Median shading loss: {med_shade:.0f}% (real reduction from trees/neighbours)") - - # Diagnostic: LiDAR point density check (averaged over buildings that had LiDAR) - print(f"\n LiDAR density check:") - if density_samples: - densities = [n / max(fp, 1) for n, fp in density_samples] - density = float(np.median(densities)) - print(f" Buildings with LiDAR: {len(density_samples)} of {len(buildings)}") - print(f" Median point density: ~{density:.2f} pts/m²") - if density < 0.5: - print(" WARNING: Low density (<0.5 pts/m²) — RANSAC results may be unreliable.") - print(" Consider finding a higher-density tile for this area.") - elif density >= 4: - print(" Density is good — RANSAC results should be reliable.") - else: - print(" Density is acceptable.") - else: - print(" No buildings overlapped the tile — density unknown.") - print(" The test bbox is likely larger than this single tile's footprint.") - - print(f"\n Results saved to: {DB_PATH}") - print(f" Open with: python -c \"import duckdb; print(duckdb.connect('{DB_PATH}').execute('SELECT * FROM homes LIMIT 5').fetchdf())\"") - con.close() - - -if __name__ == "__main__": - main() diff --git a/tests/test_config.py b/tests/test_config.py new file mode 100644 index 0000000..d697f49 --- /dev/null +++ b/tests/test_config.py @@ -0,0 +1,19 @@ +"""Config derived quantities and overrides.""" + +import pytest + +from solargrader.config import Config + + +def test_derived_quantities(): + c = Config(panel_watts=400, panel_width_m=1.05, panel_height_m=1.74) + assert c.panel_area_m2 == pytest.approx(1.827, abs=1e-3) + assert c.module_efficiency == pytest.approx(0.2189, abs=1e-3) + assert c.residential_panel_cap == 25 + + +def test_overrides_only_touch_named_fields(): + c = Config(building_source="ms", n_workers=8) + assert c.building_source == "ms" + assert c.n_workers == 8 + assert c.db_path == "solar_grader.duckdb" # untouched default diff --git a/tests/test_geometry_and_sources.py b/tests/test_geometry_and_sources.py new file mode 100644 index 0000000..857be1f --- /dev/null +++ b/tests/test_geometry_and_sources.py @@ -0,0 +1,51 @@ +"""Pure geometry / source helpers that need no network.""" + +import pytest +from shapely.geometry import box + +from solargrader.config import Config +from solargrader.geometry import polygon_area_m2 +from solargrader.solar import panels_on_plane +from solargrader.sources.buildings import _quadkey +from solargrader.sources.regions import county_bbox, resolve_targets +from solargrader.sources.tiles import _covering_cells + + +def test_polygon_area_positive_and_monotonic(): + small = polygon_area_m2(box(-76.0, 40.0, -75.999, 40.001)) + big = polygon_area_m2(box(-76.0, 40.0, -75.998, 40.002)) + assert small > 0 + assert big > small + + +def test_panels_on_plane_scales_with_area(): + cfg = Config() + assert panels_on_plane(1.0, cfg) == 0 # tiny plane erodes to nothing + forty = panels_on_plane(100.0, cfg) + assert forty > 0 + assert panels_on_plane(200.0, cfg) > forty + + +def test_quadkey_known_value(): + # Harrisburg point at zoom 9 (verified against the MS partition scheme) + assert _quadkey(-76.88, 40.27, 9) == "032010012" + + +def test_covering_cells_derives_utm_zone_from_bbox(): + cells, zone, epsg = _covering_cells((-76.93, 40.24, -76.83, 40.31)) + assert zone == 18 + assert epsg == 26918 + assert len(cells) > 0 + # Columbus, OH is UTM zone 17 — proves it isn't hardcoded. + _, oh_zone, oh_epsg = _covering_cells((-83.02, 39.95, -82.98, 39.99)) + assert oh_zone == 17 and oh_epsg == 26917 + + +def test_county_name_injection_is_rejected(): + with pytest.raises(ValueError): + county_bbox("Dauphin' OR '1'='1", Config()) + + +def test_resolve_targets_bbox_needs_no_network(): + targets = resolve_targets(Config(), bbox=(-77.0, 40.0, -76.0, 41.0)) + assert targets == [("bbox", (-77.0, 40.0, -76.0, 41.0), False)] diff --git a/tests/test_grading.py b/tests/test_grading.py new file mode 100644 index 0000000..bbae68c --- /dev/null +++ b/tests/test_grading.py @@ -0,0 +1,34 @@ +"""Grading is a pure policy function — the cheapest, most important thing to pin.""" + +from solargrader.grading import grade_home, grade_potential + + +def test_prime_south_roof_is_top_grade(): + grade, score = grade_home(res_annual_kwh=14000, res_system_kw=10, primary_azimuth_deg=180) + assert grade == "A+" + assert score >= 85 + + +def test_poorly_oriented_low_yield_roof_grades_low(): + grade, _ = grade_home(res_annual_kwh=9000, res_system_kw=10, primary_azimuth_deg=0) + assert grade in ("C", "D") + + +def test_yield_dominates_when_orientation_and_size_fixed(): + _, low = grade_home(10000, 10, 180) + _, high = grade_home(14000, 10, 180) + assert high > low + + +def test_sub_residential_system_is_penalized(): + _, big = grade_home(14000, 10, 180) # 10 kW, +10 size + _, tiny = grade_home(1400, 1, 180) # 1 kW, same yield, -10 size + assert tiny < big + + +def test_grade_potential_tiers(): + assert grade_potential(35) == "P1" + assert grade_potential(25) == "P2" + assert grade_potential(15) == "P3" + assert grade_potential(8) == "P4" + assert grade_potential(3) == "P5" diff --git a/tests/test_storage.py b/tests/test_storage.py new file mode 100644 index 0000000..3fb02ae --- /dev/null +++ b/tests/test_storage.py @@ -0,0 +1,39 @@ +"""ResultStore round-trip + resume ledgers, against a temporary DuckDB file.""" + +from shapely.geometry import box + +from solargrader.models import Building, ScoredHome +from solargrader.storage import ResultStore + + +def _home(home_id: int, grade: str = "A+") -> ScoredHome: + b = Building(id=home_id, geometry=box(0, 0, 1, 1), lat=40.0, lon=-76.0, footprint_m2=120.0) + return ScoredHome(building=b, usable_area_m2=80, primary_tilt_deg=25, + primary_azimuth_deg=180, plane_count=2, res_kwh=13000, res_kw=10, + res_panels=25, max_kwh=15000, max_kw=12, max_panels=30, + potential_grade="P3", shade_loss_pct=5.0, grade=grade, score=90) + + +def test_store_roundtrip_and_ledgers(tmp_path): + db = str(tmp_path / "t.duckdb") + with ResultStore(db) as s: + s.save(_home(1)) + s.save(_home(2, grade="B")) + s.mark_tile_done("tileA", 2) + s.mark_region_done("regionA") + assert s.count_homes() == 2 + assert s.done_tiles() == {"tileA"} + assert s.done_regions() == {"regionA"} + assert s.grade_breakdown() == {"A+": 1, "B": 1} + # State persists across reopen — the basis of resumability. + with ResultStore(db) as s: + assert s.count_homes() == 2 + assert s.done_tiles() == {"tileA"} + + +def test_save_is_idempotent(tmp_path): + db = str(tmp_path / "t.duckdb") + with ResultStore(db) as s: + s.save(_home(1)) + s.save(_home(1)) # same id → INSERT OR REPLACE, not a duplicate + assert s.count_homes() == 1 From 5bff88d1ea82dcc4642d3068b183001841cfcd83 Mon Sep 17 00:00:00 2001 From: WagnerJ-Dev Date: Sat, 4 Jul 2026 19:24:57 -0400 Subject: [PATCH 12/13] cleanup --- map.html | 48 ------------------------------------------------ 1 file changed, 48 deletions(-) delete mode 100644 map.html diff --git a/map.html b/map.html deleted file mode 100644 index 907a137..0000000 --- a/map.html +++ /dev/null @@ -1,48 +0,0 @@ - - -Solar Grader — Leads - - - - -
- \ No newline at end of file From dbcdfabc8f96d5277c2fab07e7e3af332cd4e8a5 Mon Sep 17 00:00:00 2001 From: WagnerJ-Dev Date: Sat, 4 Jul 2026 19:29:52 -0400 Subject: [PATCH 13/13] CI/CD fix --- .github/workflows/ci.yml | 3 +++ pyproject.toml | 2 ++ 2 files changed, 5 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 414b508..c3dd48f 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -25,6 +25,9 @@ jobs: run: | python -m pip install --upgrade pip pip install ruff pytest numpy scipy pandas requests shapely pyproj duckdb + # Install the package itself (no deps — the pure-function tests don't need + # the heavy optional deps, which are imported lazily) so `solargrader` imports. + pip install -e . --no-deps - name: Lint run: ruff check solargrader tests diff --git a/pyproject.toml b/pyproject.toml index 10a4d5a..52d1cf6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -49,6 +49,8 @@ include = ["solargrader*"] [tool.pytest.ini_options] testpaths = ["tests"] +# Import the package from the repo root, so `pytest` works even without an install. +pythonpath = ["."] [tool.ruff] line-length = 100