"""Diagnose why a specific parcel was/wasn't included in the buildable list. Usage: scripts/diagnose_parcel.py [GISLINK] Default: 090063 00203 (BrightRidge / Innovation Dr 300) """ from __future__ import annotations import sys from pathlib import Path import geopandas as gpd import pandas as pd from shapely import make_valid, unary_union from shapely.ops import polylabel from shapely.geometry import MultiPolygon, Polygon ROOT = Path(__file__).resolve().parent.parent CACHE = ROOT / "cache" RESIDENTIAL_DISTRICT_ZONECLASSES = { "R-1", "R-2", "R-2A", "R-2B", "R-2C", "R-3", "R-4", "R-5", "R-6", "RP-2", "RP-3", "RP-4", "RP-5", "RP-6", "RO-1", "RO-2", "A-1", "RM-3", "RM-4", "RM-5", } RESIDENTIAL_LANDUSE_CODES = {"11", "12", "13", "14"} RESIDENTIAL_PLACETYPES = { "SF RESIDENCE", "APARTMENT", "APARTMENTS", "MOBILE HOME", "CONDOMINIUM", "TOWNHOUSE", "DUPLEX", "PUBLIC HOUSING", "DORMITORY", } VACANT_PLACETYPES = {"VACANT LOT", "VACANT", "VACANT UNIT"} DC_INTERNAL_SETBACK_FT = 50 def fix(g): if g is None or g.is_empty: return g return g if g.is_valid else make_valid(g) def landuse_is_residential(x): if x is None or (isinstance(x, float) and pd.isna(x)): return False prefix = str(x).split(" - ", 1)[0] return bool({c.strip() for c in prefix.split(":") if c.strip()} & RESIDENTIAL_LANDUSE_CODES) def inscribed_radius(geom, tol=2.0): if geom is None or geom.is_empty: return 0.0 polys = [] if isinstance(geom, Polygon): polys = [geom] elif isinstance(geom, MultiPolygon): polys = list(geom.geoms) best = 0.0 for p in polys: if p.is_empty or p.area < 1: continue try: pt = polylabel(p, tolerance=tol) r = p.exterior.distance(pt) if r > best: best = r except Exception: continue return best def main(): target_gid = sys.argv[1] if len(sys.argv) > 1 else "090063 00203" print(f"Diagnosing parcel: '{target_gid}'\n") zoning = gpd.read_parquet(CACHE / "zoning.parquet") parcels = gpd.read_parquet(CACHE / "parcels.parquet") addresses = gpd.read_parquet(CACHE / "addresses.parquet") zoning["geometry"] = zoning.geometry.apply(fix) parcels["geometry"] = parcels.geometry.apply(fix) parcels = parcels[parcels.geometry.notna() & ~parcels.geometry.is_empty].copy() target = parcels[parcels["GISLINK"].astype(str) == target_gid] if target.empty: print(f"No parcel found with GISLINK '{target_gid}'") return tgt = target.iloc[0] tg = tgt.geometry print(f"Owner: {tgt.OWNER} Address: {tgt.ADDRESS}") print(f"Zone: {tgt.ZONECLASS} LANDUSE: {tgt.LANDUSE}") print(f"Parcel acres: {tg.area/43560:.2f}") # Build residential exclusion (vacant-override rule) addr_pt = addresses["PLACETYPE"].astype(str).str.upper() res_addrs = addresses[addr_pt.isin(RESIDENTIAL_PLACETYPES)] vac_addrs = addresses[addr_pt.isin(VACANT_PLACETYPES)] sj_res = gpd.sjoin(res_addrs[["geometry"]], parcels[["GISLINK","geometry"]], how="inner", predicate="within") res_gids = set(sj_res["GISLINK"].dropna().unique()) sj_vac = gpd.sjoin(vac_addrs[["geometry"]], parcels[["GISLINK","geometry"]], how="inner", predicate="within") vac_gids = set(sj_vac["GISLINK"].dropna().unique()) override = vac_gids - res_gids res_lu_mask = parcels["LANDUSE"].apply(landuse_is_residential) & ~parcels["GISLINK"].isin(override) res_lu_parcels = parcels[res_lu_mask] res_addr_parcels = parcels[parcels["GISLINK"].isin(res_gids)] zc = zoning["ZONECLASS"].fillna("").astype(str) res_zones = zoning[zc.isin(RESIDENTIAL_DISTRICT_ZONECLASSES)] i2_union = unary_union(zoning[zoning["ZONECLASS"]=="I-2"].geometry.values) # I-2 portion of this parcel i2_part = fix(tg.intersection(i2_union)) print(f"\nI-2 portion of parcel: {i2_part.area/43560:.2f} acres ({i2_part.area:.0f} sqft)") # Residential exclusion sources WITHIN 1500ft of the parcel buf_test = tg.buffer(1500) nearby_res_zones = res_zones[res_zones.intersects(buf_test)] nearby_res_lu = res_lu_parcels[res_lu_parcels.intersects(buf_test)] nearby_res_addr = res_addr_parcels[res_addr_parcels.intersects(buf_test)] print(f"\nResidential sources within 1500ft of parcel:") print(f" residential zone polygons: {len(nearby_res_zones)}") print(f" parcels with residential LU: {len(nearby_res_lu)}") print(f" parcels with dwelling addr: {len(nearby_res_addr)}") # Compute the residential exclusion union limited to nearby (faster) nearby_geoms = list(nearby_res_zones.geometry.values) + \ list(nearby_res_lu.geometry.values) + \ list(nearby_res_addr.geometry.values) if nearby_geoms: nearby_union = unary_union(nearby_geoms) else: nearby_union = None for setback in (200, 500, 1000): if nearby_union is None: buf = None remnant = i2_part else: buf = nearby_union.buffer(setback) remnant = fix(i2_part.difference(buf)) eroded = fix(remnant.buffer(-DC_INTERNAL_SETBACK_FT)) if not remnant.is_empty else remnant ir_remnant = inscribed_radius(remnant) ir_eroded = inscribed_radius(eroded) # Compute frontage ps_idx = parcels.sindex nearby_ids = list(ps_idx.query(tg.buffer(10))) nearby_parcels = [parcels.iloc[i].geometry for i in nearby_ids if parcels.iloc[i].geometry is not None and not parcels.iloc[i].geometry.equals(tg)] if nearby_parcels: nbr_union = unary_union([n.buffer(5) for n in nearby_parcels]) parent_road = tg.boundary.difference(nbr_union) else: parent_road = tg.boundary patch_frontage = 0.0 if not remnant.is_empty and not parent_road.is_empty: front = remnant.boundary.intersection(parent_road.buffer(5)) patch_frontage = front.length if not front.is_empty else 0.0 print(f"\n@ {setback} ft setback:") print(f" compliant remnant: {remnant.area/43560:.3f} ac ({remnant.area:.0f} sqft)") print(f" inscribed diam: {2*ir_remnant:.1f} ft") print(f" eroded (siteable): {eroded.area/43560:.3f} ac ({eroded.area:.0f} sqft)") print(f" inscribed diam: {2*ir_eroded:.1f} ft") print(f" patch road frontage: {patch_frontage:.1f} ft (parent has {parent_road.length:.1f} ft total)") passes_filters = (remnant.area >= 5000 and 2*ir_eroded >= 50 and patch_frontage >= 50) print(f" PASSES filters? {passes_filters}") if __name__ == "__main__": main()