153 lines
6.6 KiB
Python
153 lines
6.6 KiB
Python
"""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()
|