"""Boston BPDA OBJ collections with explicit archive, status and coordinate identity.
The official source is tiled and version-pinned. A local, hash-pinned catalog
can select downloaded/historical archives without changing the Scene API.
It is never an implicit fallback from a failed current-source request.
"""
from __future__ import annotations
import asyncio
import csv
import hashlib
import io
import json
import math
import re
import zipfile
from pathlib import Path
import httpx
import numpy as np
from pyproj import Transformer
from shapely.geometry import Polygon, box, shape
from shapely.ops import transform
from openworld_radio_twin.geodesy import LocalFrame
from openworld_radio_twin.providers.base import BuildingProviderError, BuildingQueryResult
from openworld_radio_twin.providers.cityjson import triangulate_surface
from openworld_radio_twin.providers.native_mesh import mesh_feature
ATTRIBUTION = "Boston Planning and Development Agency / Boston Planning Department"
LICENSE_URL = "https://www.bostonplans.org/3d-data-maps/3d-smart-model/3d-data-download"
RELEASE = "20250128"
TILE_QUERY = (
"https://services.arcgis.com/sFnw0xNflSi8J0uh/ArcGIS/rest/services/"
"Bos3D_TileGrid_PopUps/FeatureServer/70/query"
)
DOWNLOAD_ROOT = f"https://maps.bostonplans.org/3d/Bos3d_BldgModels_{RELEASE}_OBJ"
FOOT_M = 1200 / 3937 # EPSG:9003 US survey foot, including NAVD88 height.
MAX_ARCHIVE_BYTES = 256 * 1024 * 1024
MAX_EXPANDED_BYTES = 1024 * 1024 * 1024
EXISTING_STATUSES = {"Current", "Construction Complete", "Approved Demo", "Permitted Demo"}
TO_STATE_PLANE = Transformer.from_crs(4326, 2249, always_xy=True)
[docs]
def configuration_identity(settings):
path = settings.boston_source_catalog
if path is None:
return {"release": RELEASE}
return {"catalog_sha256": hashlib.sha256(path.read_bytes()).hexdigest()}
[docs]
def json_document(content: bytes):
"""Read JSON or a single JS variable assignment, never execute JavaScript."""
text = content.decode("utf-8-sig").strip()
text = re.sub(r"^(?:var|let|const)\s+\w+\s*=\s*", "", text)
return json.loads(text.rstrip(";\r\n "))
def _checked_archive(content):
archive = zipfile.ZipFile(io.BytesIO(content))
if sum(member.file_size for member in archive.infolist()) > MAX_EXPANDED_BYTES:
archive.close()
raise ValueError("Boston archive exceeds the expanded size limit")
return archive
def _model_obj(archive, members, identifier):
package_name = members.get(f"{identifier}_OBJ.zip")
if package_name:
with _checked_archive(archive.read(package_name)) as package:
objs = [n for n in package.namelist() if n.lower().endswith(".obj")]
if len(objs) != 1:
raise ValueError("Building package must contain exactly one OBJ")
return package.read(objs[0])
return archive.read(members[f"{identifier}.obj"])
[docs]
def catalog_offset(archive, catalog, members):
"""Resolve the common translation from paired source footprint/OBJ bounds.
Historical Boston tile-frame metadata can contradict the model coordinates.
Require agreement of three independent complete model extents within 0.25 ft;
do not move buildings independently or infer scale/rotation from a picture.
"""
candidates = []
identifiers = []
for feature in catalog["features"]:
attributes = feature["properties"]
if attributes.get("Model_LOD") not in {1, "1"}:
continue
identifier = str(attributes["Model_ID"])
try:
obj = _model_obj(archive, members, identifier)
points = np.asarray(
[
list(map(float, line.split()[1:4]))
for line in obj.decode("utf-8-sig").splitlines()
if line.startswith("v ")
]
)
footprint = transform(TO_STATE_PLANE.transform, shape(feature["geometry"]))
target = np.asarray(footprint.bounds)
original = np.r_[points[:, :2].min(axis=0), points[:, :2].max(axis=0)]
delta = target - original
candidate = np.round((delta[:2] + delta[2:]) / 2)
if np.max(np.abs(delta - np.tile(candidate, 2))) > 0.25:
continue
candidates.append(candidate)
identifiers.append(identifier)
except (KeyError, ValueError, zipfile.BadZipFile):
continue
if len(candidates) == 3:
break
if len(candidates) < 3 or not np.allclose(candidates, candidates[0], atol=0.01, rtol=0):
raise ValueError("Boston catalog cannot establish a consistent model-coordinate origin")
return candidates[0].tolist(), identifiers
[docs]
def decode_obj(content: bytes, metadata: dict, offset_xy_ft, *, source_url: str):
"""Preserve OBJ faces; record normal-derived proxies, not source semantics."""
points, polygons = [], []
for line in content.decode("utf-8-sig").splitlines():
words = line.split("#", 1)[0].split()
if not words:
continue
if words[0] == "v":
if len(words) != 4:
raise ValueError("Boston OBJ requires three-coordinate vertices")
points.append(tuple(map(float, words[1:])))
elif words[0] == "f":
face = []
for word in words[1:]:
value = int(word.split("/", 1)[0])
if value == 0:
raise ValueError("OBJ vertex indices cannot be zero")
face.append(value - 1 if value > 0 else len(points) + value)
polygons.append(face)
vertices = np.asarray(points, dtype=float)
if vertices.ndim != 2 or vertices.shape[1] != 3 or not np.isfinite(vertices).all():
raise ValueError("Boston OBJ has invalid coordinates")
if len(offset_xy_ft) != 2 or not np.isfinite(offset_xy_ft).all():
raise ValueError("Boston coordinate offset must contain two finite values")
# XY stays in ftUS (EPSG:2249); BuildingGeometry Z is always metric.
vertices[:, :2] += offset_xy_ft
metric = vertices * FOOT_M
vertices[:, 2] *= FOOT_M
triangles, kinds, omitted = [], [], 0
for polygon in polygons:
faces = triangulate_surface(metric, [polygon])
omitted += not faces
for face in faces:
a, b, c = metric[list(face)]
normal = np.cross(b - a, c - a)
slope = normal[2] / np.linalg.norm(normal)
kinds.append(
"WallSurface"
if abs(slope) < 0.1
else "RoofSurface"
if slope > 0
else "GroundSurface"
)
triangles.append(face)
identifier = str(metadata["Model_ID"])
building = mesh_feature(
identifier=identifier,
source="boston",
parent_id=identifier,
vertices=vertices,
triangles=triangles,
surface_types=kinds,
crs="EPSG:2249",
vertical_datum="NAVD88",
lod=str(metadata.get("Model_LOD", "unspecified")),
attributes={
**metadata,
"attribution": ATTRIBUTION,
"license_url": LICENSE_URL,
"source_url": source_url,
"source_offset_xy_ft": list(offset_xy_ft),
"source_vertical_unit": "US survey foot",
"surface_semantics": "derived_from_source_face_normals_nz_threshold_0.1",
"omitted_zero_area_surface_count": omitted,
},
)
ground = metadata.get("Gnd_El_Ft")
if ground is not None and str(ground).strip():
ground = float(ground) * FOOT_M
if not math.isfinite(ground):
raise ValueError("Boston source ground elevation must be finite")
building.geometry.ground_elevation_m = ground
building.min_height_m = max(0.0, vertices[:, 2].min() - ground)
building.source_attributes["ground_reference"] = "source_Gnd_El_Ft"
else:
building.source_attributes["ground_reference"] = "lowest_downward_facing_source_surface"
building.name = metadata.get("Name")
return building
[docs]
def decode_collection(content: bytes, *, domain, offset_xy_ft, source_url):
"""Filter existing buildings, then decode their OBJ packages without extracting ZIPs."""
records, warnings = [], []
with _checked_archive(content) as archive:
catalogs = [
name
for name in archive.namelist()
if Path(name).name in {"ModelCatalog_geojson.js", "ModelCatalog.geojson"}
]
members = {Path(name).name: name for name in archive.namelist() if not name.endswith("/")}
if catalogs:
catalog = json_document(archive.read(catalogs[0]))
entries = [(f["properties"], shape(f["geometry"])) for f in catalog["features"]]
# Validate rather than silently trusting legacy frame offsets.
if len(entries) >= 3:
actual_offset, anchors = catalog_offset(archive, catalog, members)
if not np.allclose(actual_offset, offset_xy_ft, atol=0.01, rtol=0):
raise ValueError(
f"Boston supplied origin {offset_xy_ft} conflicts with model catalog "
f"origin {actual_offset}; regenerate the explicit source catalog"
)
else:
anchors = []
else:
csv_names = [name for name in archive.namelist() if name.lower().endswith(".csv")]
if len(csv_names) != 1:
raise ValueError("Boston archive requires a model GeoJSON or CSV catalog")
rows = csv.DictReader(io.StringIO(archive.read(csv_names[0]).decode("utf-8-sig")))
entries = [(dict(row), None) for row in rows]
for attributes, footprint in entries:
if attributes.get("Status") not in EXISTING_STATUSES:
continue
if str(attributes.get("StructType", "Building")).lower() not in {"building", "bldg"}:
continue
if footprint is not None and not footprint.intersects(domain):
continue
identifier = str(attributes.get("Model_ID", ""))
try:
obj = _model_obj(archive, members, identifier)
building = decode_obj(obj, attributes, offset_xy_ft, source_url=source_url)
if catalogs:
building.source_attributes["origin_validation_models"] = anchors
if Polygon(building.coordinates, building.holes).intersects(domain):
records.append(building)
except (ValueError, KeyError, TypeError, zipfile.BadZipFile) as exc:
warnings.append(f"Skipped Boston model {identifier}: {exc}")
return records, warnings
[docs]
class BostonProvider:
def __init__(self, settings, client: httpx.AsyncClient):
self.client = client
self.catalog_path = settings.boston_source_catalog
self._lock = asyncio.Lock()
async def _download(self, url, maximum=MAX_ARCHIVE_BYTES):
chunks, size = [], 0
async with self.client.stream("GET", url, timeout=120) as response:
response.raise_for_status()
async for chunk in response.aiter_bytes():
size += len(chunk)
if size > maximum:
raise ValueError("Boston download exceeds the size limit")
chunks.append(chunk)
return b"".join(chunks)
async def _tiles(self, domain):
if self.catalog_path is not None:
document = json.loads(self.catalog_path.read_text())
if document.get("format") != "owrt-boston-archives-v1":
raise ValueError("Unsupported Boston source catalog format")
tiles = [tile for tile in document["tiles"] if box(*tile["bbox"]).intersects(domain)]
return tiles, str(document["release"])
response = await self.client.get(
TILE_QUERY,
params={
"f": "geojson",
"where": "1=1",
"outFields": "*",
"outSR": "4326",
"geometry": ",".join(map(str, domain.bounds)),
"geometryType": "esriGeometryEnvelope",
"inSR": "4326",
"spatialRel": "esriSpatialRelIntersects",
"resultRecordCount": 1000,
},
timeout=60,
)
response.raise_for_status()
data = response.json()
if "error" in data or data.get("exceededTransferLimit"):
raise ValueError("Boston tile index query failed or was truncated")
tiles = []
for feature in data["features"]:
attributes = feature["properties"]
tile_id = attributes["tile_id"]
if not re.fullmatch(r"BOS_[A-Z]+_\d+", tile_id):
raise ValueError("Invalid Boston tile identifier")
tiles.append(
{
"id": tile_id,
"url": f"{DOWNLOAD_ROOT}/{tile_id}_BldgModels_OBJ.zip",
"offset_xy_ft": [
attributes["MASP_X"] - attributes["BosShift_X"],
attributes["MASP_Y"] - attributes["BosShift_Y"],
],
}
)
return tiles, RELEASE
[docs]
async def buildings(self, latitude, longitude, radius_m, limit):
if limit < 1 or not math.isfinite(radius_m) or not 0 < radius_m <= 5000:
raise ValueError("Boston queries require a positive limit and radius at most 5000 m")
if not (-90 <= latitude <= 90 and -180 <= longitude <= 180):
raise ValueError("Invalid latitude/longitude")
domain = Polygon(
LocalFrame.at(longitude, latitude).corners(-radius_m, -radius_m, radius_m, radius_m)
)
buildings, warnings, seen = [], [], set()
async with self._lock:
try:
tiles, release = await self._tiles(domain)
if not tiles:
raise ValueError("Query does not intersect the selected Boston source coverage")
if len(tiles) > 64:
raise ValueError("Boston query exceeds 64 tiles")
for tile in tiles:
if "path" in tile:
path = self.catalog_path.parent / tile["path"]
if path.stat().st_size > MAX_ARCHIVE_BYTES:
raise ValueError("Boston local archive exceeds size limit")
content = await asyncio.to_thread(path.read_bytes)
else:
content = await self._download(tile["url"])
digest = hashlib.sha256(content).hexdigest()
if self.catalog_path is not None and digest != tile.get("sha256"):
raise ValueError("Boston archive SHA-256 differs from its source catalog")
records, rejected = await asyncio.to_thread(
decode_collection,
content,
domain=domain,
offset_xy_ft=tile["offset_xy_ft"],
source_url=tile.get("url", tile.get("path")),
)
warnings.extend(rejected)
for record in records:
if record.feature_id in seen:
continue
seen.add(record.feature_id)
record.source_release = release
record.source_attributes.update(source_sha256=digest, tile_id=tile["id"])
buildings.append(record)
if len(buildings) >= limit:
warnings.append(f"Boston building limit ({limit}) reached")
break
if len(buildings) >= limit:
break
except (
httpx.HTTPError,
ValueError,
KeyError,
TypeError,
OSError,
zipfile.BadZipFile,
) as exc:
raise BuildingProviderError(f"Boston building acquisition failed: {exc}") from exc
return BuildingQueryResult(
buildings=buildings, provider="Boston BPDA OBJ", release=release, warnings=warnings
)
