.zip` with the
+ GeoPackage inside. The largest `.gpkg` member wins if there are several;
+ anything else in the archive (metadata PDFs, licence texts) is ignored.
+ """
+ import io
+ import zipfile
+
+ try:
+ with zipfile.ZipFile(io.BytesIO(payload)) as zf:
+ members = [i for i in zf.infolist() if i.filename.lower().endswith(".gpkg")]
+ if not members:
+ raise BuildingLookupError(
+ "the building tile ZIP holds no GeoPackage"
+ )
+ member = max(members, key=lambda i: i.file_size)
+ return zf.read(member)
+ except zipfile.BadZipFile as exc:
+ raise BuildingLookupError(
+ f"the building tile was announced as ZIP but did not open: {exc}"
+ ) from exc
+
+
+def _features_from_item(
+ item: StacItem,
+ client: GeotorgetClient | None = None,
+ window_boxes: list[tuple[float, float, float, float]] | None = None,
+) -> list[dict[str, Any]]:
+ """Every building-like feature an item carries.
+
+ A STAC item may *be* the building -- geometry inline, no download -- or it
+ may be a tile whose asset holds thousands of them. Lantmaeteriet's live
+ catalogue publishes `byggnader` as one zipped **GeoPackage** per
+ municipality, so the asset path is the normal one and the inline path the
+ exception — and `window_boxes` (the search window, in each frame the file
+ could be in) keeps a whole city from being decoded for a 150 m search.
+ """
+ # A usable data asset wins over the item's own geometry: a tile item's
+ # geometry is the TILE's outline — for Lantmaeteriet, the whole
+ # municipality — and reading it as a building both invents a footprint
+ # nobody has and skips the real ones in the asset. Inline geometry is the
+ # fallback for catalogues whose items *are* the buildings.
+ asset = item.pick(MEDIA_GEOPACKAGE, MEDIA_ZIP, MEDIA_GEOJSON)
+ if asset is None or not asset.href or client is None:
+ geom = (item.raw or {}).get("geometry")
+ if geom:
+ return [{"geometry": geom, "properties": (item.raw or {}).get("properties") or {},
+ "id": item.item_id}]
+ return []
+ media = asset.effective_media_type
+ payload = client.download(asset.href)
+ if media == MEDIA_GEOJSON:
+ return _features_from_geojson(payload)
+ if payload[:4] == b"PK\x03\x04":
+ payload = _gpkg_from_zip(payload)
+ try:
+ return read_features(payload, bboxes=window_boxes)
+ except GeoPackageError as exc:
+ raise BuildingLookupError(
+ f"the building tile for this site could not be read: {exc}"
+ ) from exc
+
+
+def _features_from_geojson(payload: bytes) -> list[dict[str, Any]]:
+ """A GeoJSON FeatureCollection asset, for catalogues that publish one."""
+ import json
+
+ try:
+ doc = json.loads(payload.decode("utf-8"))
+ except (UnicodeDecodeError, ValueError) as exc:
+ raise BuildingLookupError(
+ f"the building tile was announced as GeoJSON but did not parse: {exc}"
+ ) from exc
+ if isinstance(doc, dict) and doc.get("type") == "FeatureCollection":
+ return list(doc.get("features") or [])
+ if isinstance(doc, dict) and doc.get("type") == "Feature":
+ return [doc]
+ return []
+
+
+def buildings_from_features(
+ features: Iterable[dict[str, Any]],
+ *,
+ latitude: float,
+ longitude: float,
+ fallback_id: str = "building",
+) -> list[Building]:
+ """Turn GeoJSON-ish features into ranked Building candidates."""
+ site_n, site_e = sweref.wgs84_to_sweref99tm(latitude, longitude)
+ out: list[Building] = []
+ for i, feat in enumerate(features):
+ for j, ring in enumerate(_rings_from_geometry(feat.get("geometry") or {})):
+ ring_sweref = _to_sweref(ring)
+ area = _shoelace_area(ring_sweref)
+ if area < MIN_FOOTPRINT_AREA_M2 or area > MAX_FOOTPRINT_AREA_M2:
+ continue
+ cx, cy = _centroid(ring_sweref)
+ props = dict(feat.get("properties") or {})
+ bid = str(feat.get("id") or props.get("objektidentitet") or f"{fallback_id}-{i}")
+ if j:
+ bid = f"{bid}-{j}"
+ out.append(Building(
+ building_id=bid,
+ ring_sweref=ring_sweref,
+ area_m2=area,
+ distance_m=math.hypot(cx - site_e, cy - site_n),
+ properties={k: v for k, v in props.items() if isinstance(v, (str, int, float))},
+ ))
+ out.sort(key=lambda b: b.distance_m)
+ return out
+
+
+def search_buildings(
+ client: GeotorgetClient,
+ *,
+ latitude: float,
+ longitude: float,
+ radius_m: float = DEFAULT_SEARCH_RADIUS_M,
+ limit: int = 50,
+ collection: str = COLLECTION_BUILDINGS,
+ bbox_epsg: int = 4326,
+) -> list[Building]:
+ """Building footprints near a site, nearest first."""
+ bbox = sweref.stac_search_bbox(latitude, longitude, radius_m, bbox_epsg)
+ # The same window in every frame and axis order a tile could be stored
+ # in. Coordinate magnitude keeps the frames apart (degrees never look
+ # like metres), and carrying both axis orders costs nothing but a few
+ # extra decoded rows — a box only ever widens what is kept.
+ def swapped(b):
+ return (b[1], b[0], b[3], b[2])
+
+ wgs = sweref.stac_search_bbox(latitude, longitude, radius_m, 4326)
+ swe = sweref.stac_search_bbox(latitude, longitude, radius_m, 3006)
+ window_boxes = [wgs, swapped(wgs), swe, swapped(swe)]
+ items = client.search(collection, bbox, limit=limit)
+ features: list[dict[str, Any]] = []
+ for item in items:
+ features.extend(_features_from_item(item, client, window_boxes))
+ if not features:
+ raise BuildingLookupError(
+ "no building footprints were returned for this site. The Geotorget "
+ "account needs access to 'Byggnad Nedladdning, vektor', and the data "
+ "covers Sweden only."
+ )
+ return buildings_from_features(
+ features, latitude=latitude, longitude=longitude
+ )
+
+
+def point_in_ring(x: float, y: float, ring: list[tuple[float, float]]) -> bool:
+ """Ray-casting point-in-polygon, in the ring's own projected frame."""
+ inside = False
+ n = len(ring)
+ for i in range(n):
+ x1, y1 = ring[i]
+ x2, y2 = ring[(i + 1) % n]
+ if (y1 > y) != (y2 > y):
+ if y2 != y1 and x < x1 + (y - y1) * (x2 - x1) / (y2 - y1):
+ inside = not inside
+ return inside
+
+
+def inflate_ring(ring: list[tuple[float, float]], metres: float) -> list[tuple[float, float]]:
+ """Push a ring outward from its centroid by roughly `metres`.
+
+ Approximate on purpose: a true polygon offset needs a geometry library, and
+ this only has to catch the eaves. For the compact, roughly convex outlines
+ houses actually have, scaling about the centroid is within a few centimetres
+ of a proper offset; for a long thin wing it over-buffers the short sides,
+ which costs a little neighbouring ground rather than losing roof.
+ """
+ if metres <= 0 or len(ring) < 3:
+ return ring
+ cx, cy = _centroid(ring)
+ out = []
+ for x, y in ring:
+ dx, dy = x - cx, y - cy
+ d = math.hypot(dx, dy)
+ if d < 1e-9:
+ out.append((x, y))
+ continue
+ scale = (d + metres) / d
+ out.append((cx + dx * scale, cy + dy * scale))
+ return out
+
+
+def clip_to_footprint(points, ring: list[tuple[float, float]],
+ *, buffer_m: float = DEFAULT_EAVES_BUFFER_M):
+ """Keep only the returns standing on one building.
+
+ `points` is (N, 3) in SWEREF 99 TM metres, the frame the LiDAR arrives in.
+ """
+ import numpy as np
+
+ pts = np.asarray(points, dtype=float)
+ if len(pts) == 0 or len(ring) < 3:
+ return pts
+ outline = inflate_ring(ring, buffer_m)
+
+ # Cheap rejection first: most of a tile is nowhere near the building.
+ xs = [p[0] for p in outline]
+ ys = [p[1] for p in outline]
+ box = (
+ (pts[:, 0] >= min(xs)) & (pts[:, 0] <= max(xs))
+ & (pts[:, 1] >= min(ys)) & (pts[:, 1] <= max(ys))
+ )
+ candidates = np.nonzero(box)[0]
+ keep = [i for i in candidates if point_in_ring(pts[i, 0], pts[i, 1], outline)]
+ return pts[keep]
diff --git a/roofmodel/ftw_roofmodel/geopackage.py b/roofmodel/ftw_roofmodel/geopackage.py
new file mode 100644
index 000000000..76e3c757c
--- /dev/null
+++ b/roofmodel/ftw_roofmodel/geopackage.py
@@ -0,0 +1,303 @@
+"""Read polygon features out of a GeoPackage, using only the standard library.
+
+Lantmaeteriet publishes *Byggnad Nedladdning, vektor* as GeoPackage, so a STAC
+item for a building tile carries a `.gpkg` asset rather than inline GeoJSON.
+
+A GeoPackage is a SQLite database with an agreed set of metadata tables, and
+geometry stored as a small binary header followed by standard WKB. Both are
+published specifications with fixed layouts, so decoding them here costs about
+a hundred lines of `sqlite3` and `struct` -- against a GDAL/fiona/geopandas
+dependency that would not install on a Pi without a compiler and pulls in a
+second projection stack we already decided not to carry (see sweref.py, which
+implements SWEREF 99 TM directly for the same reason).
+
+Only what a roof model needs is read: polygon and multipolygon rings, in the
+file's own coordinate reference system. Curves, triangulated surfaces and the
+extended (`GPB` "extended geometry") binary types are rejected explicitly rather
+than mis-parsed, because a wrong ring silently clips the wrong LiDAR.
+
+References
+----------
+GeoPackage Encoding Standard (OGC 12-128r19), clause 2.1.3 "BLOB Format".
+OpenGIS Simple Features (OGC 06-103r4), clause 8.2 "Well-known Binary".
+"""
+
+from __future__ import annotations
+
+import os
+import sqlite3
+import struct
+import tempfile
+from typing import Any, Iterator
+
+# "GP" -- the two magic bytes every GeoPackage geometry blob starts with.
+GPKG_MAGIC = b"GP"
+
+# Envelope sizes in doubles, indexed by the header's envelope indicator.
+# 0 = absent, 1 = xy, 2 = xyz, 3 = xym, 4 = xyzm. 5-7 are reserved.
+_ENVELOPE_DOUBLES = {0: 0, 1: 4, 2: 6, 3: 6, 4: 8}
+
+# WKB geometry type codes we can use. The ISO variants add 1000 for Z, 2000 for
+# M and 3000 for ZM, so the base code is recovered with % 1000.
+_WKB_POLYGON = 3
+_WKB_MULTIPOLYGON = 6
+
+# The EWKB flag bits PostGIS adds to the type word. GeoPackage forbids them, but
+# files written by other tools do turn up, and reading a flagged type as a
+# geometry code would silently produce nonsense.
+_EWKB_Z = 0x80000000
+_EWKB_M = 0x40000000
+_EWKB_SRID = 0x20000000
+
+
+class GeoPackageError(ValueError):
+ """The file is not a GeoPackage, or holds geometry we will not guess at."""
+
+
+def _unpack(fmt: str, data: bytes, offset: int) -> tuple[Any, ...]:
+ size = struct.calcsize(fmt)
+ if offset + size > len(data):
+ raise GeoPackageError("geometry blob ended mid-value")
+ return struct.unpack_from(fmt, data, offset)
+
+
+def _dimensions(type_word: int) -> tuple[int, int]:
+ """(base geometry code, coordinates per point) for a WKB type word."""
+ has_z = bool(type_word & _EWKB_Z)
+ has_m = bool(type_word & _EWKB_M)
+ code = type_word & ~(_EWKB_Z | _EWKB_M | _EWKB_SRID)
+ # ISO style: 1000/2000/3000 offsets carry the same information.
+ if code >= 3000:
+ code, has_z, has_m = code - 3000, True, True
+ elif code >= 2000:
+ code, has_m = code - 2000, True
+ elif code >= 1000:
+ code, has_z = code - 1000, True
+ return code, 2 + int(has_z) + int(has_m)
+
+
+def _read_ring(data: bytes, offset: int, endian: str, coords: int) -> tuple[list[tuple[float, float]], int]:
+ (count,) = _unpack(endian + "I", data, offset)
+ offset += 4
+ stride = 8 * coords
+ ring: list[tuple[float, float]] = []
+ for _ in range(count):
+ x, y = _unpack(endian + "dd", data, offset)
+ ring.append((x, y))
+ offset += stride
+ return ring, offset
+
+
+def _read_polygon(data: bytes, offset: int, endian: str, coords: int) -> tuple[list[list[tuple[float, float]]], int]:
+ (n_rings,) = _unpack(endian + "I", data, offset)
+ offset += 4
+ rings = []
+ for _ in range(n_rings):
+ ring, offset = _read_ring(data, offset, endian, coords)
+ rings.append(ring)
+ return rings, offset
+
+
+def _read_geometry(data: bytes, offset: int) -> tuple[list[list[list[tuple[float, float]]]], int]:
+ """Read one WKB geometry, returning it as a list of polygons."""
+ (byte_order,) = _unpack("B", data, offset)
+ endian = "<" if byte_order == 1 else ">"
+ (type_word,) = _unpack(endian + "I", data, offset + 1)
+ offset += 5
+ if type_word & _EWKB_SRID:
+ offset += 4 # embedded SRID, which we take from the header instead
+ code, coords = _dimensions(type_word)
+ if code == _WKB_POLYGON:
+ rings, offset = _read_polygon(data, offset, endian, coords)
+ return [rings], offset
+ if code == _WKB_MULTIPOLYGON:
+ (n,) = _unpack(endian + "I", data, offset)
+ offset += 4
+ polys = []
+ for _ in range(n):
+ # Each part carries its own byte order and type word.
+ part, offset = _read_geometry(data, offset)
+ polys.extend(part)
+ return polys, offset
+ raise GeoPackageError(
+ f"WKB geometry type {code} is not a polygon; a building footprint must be "
+ "a Polygon or MultiPolygon"
+ )
+
+
+def parse_geometry_blob(blob: bytes) -> dict[str, Any] | None:
+ """Decode a GeoPackage geometry BLOB into a GeoJSON-shaped geometry.
+
+ Returns None for the empty geometry, which GeoPackage represents with a flag
+ rather than an absent row.
+ """
+ if len(blob) < 8 or blob[:2] != GPKG_MAGIC:
+ raise GeoPackageError("not a GeoPackage geometry blob (bad magic)")
+ flags = blob[3]
+ if flags & 0x20:
+ raise GeoPackageError(
+ "extended (ExtendedGeoPackageBinary) geometry is not supported"
+ )
+ envelope_indicator = (flags >> 1) & 0x07
+ if envelope_indicator not in _ENVELOPE_DOUBLES:
+ raise GeoPackageError(f"reserved envelope indicator {envelope_indicator}")
+ if flags & 0x10: # empty geometry
+ return None
+ offset = 8 + 8 * _ENVELOPE_DOUBLES[envelope_indicator]
+ polygons, _ = _read_geometry(blob, offset)
+ if not polygons:
+ return None
+ if len(polygons) == 1:
+ return {"type": "Polygon", "coordinates": [[list(p) for p in r] for r in polygons[0]]}
+ return {
+ "type": "MultiPolygon",
+ "coordinates": [[[list(p) for p in r] for r in poly] for poly in polygons],
+ }
+
+
+def _feature_tables(conn: sqlite3.Connection) -> list[tuple[str, str]]:
+ """(table, geometry column) for every feature table in the file."""
+ try:
+ rows = conn.execute(
+ "SELECT c.table_name, g.column_name FROM gpkg_contents c "
+ "JOIN gpkg_geometry_columns g ON g.table_name = c.table_name "
+ "WHERE c.data_type = 'features'"
+ ).fetchall()
+ except sqlite3.DatabaseError as exc:
+ raise GeoPackageError(f"not a readable GeoPackage: {exc}") from exc
+ return [(str(t), str(c)) for t, c in rows]
+
+
+def _blob_envelope(blob: bytes) -> tuple[float, float, float, float] | None:
+ """The header envelope of a geometry blob as (minx, miny, maxx, maxy).
+
+ GeoPackage stores it as [minx, maxx, miny, maxy] doubles right after the
+ 8-byte header (OGC 12-128r19, clause 2.1.3), in the header's own byte
+ order. None when the writer chose not to include one.
+ """
+ if len(blob) < 8 or blob[:2] != GPKG_MAGIC:
+ return None
+ flags = blob[3]
+ if ((flags >> 1) & 0x07) == 0 or flags & 0x10:
+ return None
+ endian = "<" if flags & 0x01 else ">"
+ try:
+ minx, maxx, miny, maxy = _unpack(endian + "dddd", blob, 8)
+ except GeoPackageError:
+ return None
+ return (minx, miny, maxx, maxy)
+
+
+def _intersects(a: tuple[float, float, float, float], b: tuple[float, float, float, float]) -> bool:
+ return a[0] <= b[2] and a[2] >= b[0] and a[1] <= b[3] and a[3] >= b[1]
+
+
+def _geometry_bounds(geometry: dict[str, Any]) -> tuple[float, float, float, float] | None:
+ """(minx, miny, maxx, maxy) over every ring of a parsed geometry.
+
+ The fallback when a writer omitted the header envelope — Lantmaeteriet's
+ municipality files do (flags 0x00 on every blob), so without this the
+ bbox filter would keep all 90k+ rows and the row limit would truncate the
+ table before it ever reached the site.
+ """
+ polys = geometry.get("coordinates") or []
+ if geometry.get("type") == "Polygon":
+ polys = [polys]
+ xs: list[float] = []
+ ys: list[float] = []
+ for poly in polys:
+ for ring in poly:
+ for point in ring:
+ xs.append(point[0])
+ ys.append(point[1])
+ if not xs:
+ return None
+ return (min(xs), min(ys), max(xs), max(ys))
+
+
+def read_features(
+ data: bytes,
+ *,
+ limit: int = 5000,
+ bboxes: list[tuple[float, float, float, float]] | None = None,
+) -> list[dict[str, Any]]:
+ """Every polygon feature in a GeoPackage, as GeoJSON-shaped dicts.
+
+ Attributes travel alongside the geometry so the picker can label a building
+ with whatever the source calls it.
+
+ `bboxes` filters by the geometry blobs' header envelopes: a row is kept
+ when its envelope intersects ANY of the boxes (each (minx, miny, maxx,
+ maxy)). More than one box exists because the file's CRS isn't declared to
+ this reader — the caller passes the same window in every frame the file
+ could be in, and the frames' coordinate magnitudes are so far apart
+ (degrees vs. six-figure metres) that only the matching one can intersect.
+ Lantmaeteriet ships one GeoPackage per *municipality*, so without a filter
+ a 150 m search would decode a whole city.
+ """
+ return list(iter_features(data, limit=limit, bboxes=bboxes))
+
+
+def iter_features(
+ data: bytes,
+ *,
+ limit: int = 5000,
+ bboxes: list[tuple[float, float, float, float]] | None = None,
+) -> Iterator[dict[str, Any]]:
+ if not data.startswith(b"SQLite format 3\x00"):
+ raise GeoPackageError(
+ "asset is not a GeoPackage (missing the SQLite file header)"
+ )
+ # sqlite3 opens paths, not buffers, and a GeoPackage is random-access by
+ # design, so the bytes land in a temp file for the life of the read.
+ fd, path = tempfile.mkstemp(suffix=".gpkg")
+ try:
+ with os.fdopen(fd, "wb") as fh:
+ fh.write(data)
+ conn = sqlite3.connect(path)
+ try:
+ conn.row_factory = sqlite3.Row
+ yielded = 0
+ for table, geom_col in _feature_tables(conn):
+ cursor = conn.execute(f'SELECT * FROM "{table}"')
+ for row in cursor:
+ if yielded >= limit:
+ return
+ blob = row[geom_col]
+ if not isinstance(blob, (bytes, bytearray)):
+ continue
+ env = None
+ if bboxes:
+ env = _blob_envelope(bytes(blob))
+ if env is not None and not any(_intersects(env, b) for b in bboxes):
+ continue
+ try:
+ geometry = parse_geometry_blob(bytes(blob))
+ except GeoPackageError:
+ # One unreadable row must not lose the other buildings
+ # in the tile.
+ continue
+ if geometry is None:
+ continue
+ if bboxes and env is None:
+ bounds = _geometry_bounds(geometry)
+ if bounds is not None and not any(_intersects(bounds, b) for b in bboxes):
+ continue
+ props = {
+ k: row[k]
+ for k in row.keys()
+ if k != geom_col and isinstance(row[k], (str, int, float))
+ }
+ yield {
+ "id": str(props.get("objektidentitet") or f"{table}-{yielded}"),
+ "geometry": geometry,
+ "properties": props,
+ }
+ yielded += 1
+ finally:
+ conn.close()
+ finally:
+ try:
+ os.unlink(path)
+ except OSError: # pragma: no cover - Windows may hold the handle briefly
+ pass
diff --git a/roofmodel/ftw_roofmodel/geotorget.py b/roofmodel/ftw_roofmodel/geotorget.py
new file mode 100644
index 000000000..bbe4a13e6
--- /dev/null
+++ b/roofmodel/ftw_roofmodel/geotorget.py
@@ -0,0 +1,326 @@
+"""STAC catalog access: authentication, search, and asset download.
+
+The default catalog is Lantmaeteriet's Geotorget, where two products are used,
+both free open data (CC BY 4.0) but both gated behind a Geotorget account the
+operator orders themselves:
+
+ * *Byggnad Nedladdning, vektor* -- building footprint polygons.
+ * *Laserdata Nedladdning, Skog* -- airborne LiDAR, 1-2 points/m2, from 2018.
+
+Authentication is HTTP Basic with the operator's own Geotorget account
+username and password: Lantmaeteriet provides no OAuth for its STAC download
+APIs, so the account credential is the only door. Credentials are the
+operator's own and are never shipped, logged or echoed back through the API.
+FTW stores them the same way it stores `weather.api_key`, and redacts them in
+config responses.
+
+Nothing below is Lantmaeteriet-specific beyond the defaults: search is the
+standard `POST {base}/search` of the STAC API spec and downloads follow asset
+hrefs, so any STAC-conformant catalog behind Basic auth (or none) works by
+pointing `base_url` and the collection ids elsewhere. The one non-standard
+wrinkle is the search bbox CRS -- the spec mandates WGS84, Lantmaeteriet
+expects SWEREF 99 TM -- which is why callers choose the bbox they send (see
+`--bbox-epsg` in __main__).
+
+Only `requests` is used. The STAC API is plain JSON over HTTP, so pulling in
+pystac-client would add a dependency for a search body we can write in six
+lines -- and a thinner surface is easier to keep working when Lantmaeteriet
+moves an endpoint.
+"""
+
+from __future__ import annotations
+
+import dataclasses
+import datetime as dt
+from typing import Any, Iterable
+
+# Roots of the STAC APIs, as verified against the live service (2026-09-02):
+# Lantmaeteriet does not serve one catalogue — vector products and elevation
+# products have separate STAC roots. The standard search endpoint is
+# POST {base}/search on each, it takes a WGS84 (lon/lat) bbox per the STAC
+# spec, and the catalogue metadata is anonymously readable; the Geotorget
+# credentials are enforced on the asset downloads (dl1.lantmateriet.se
+# answers 401 without them).
+DEFAULT_BASE_URL = "https://api.lantmateriet.se/stac-vektor/v1"
+DEFAULT_LIDAR_BASE_URL = "https://api.lantmateriet.se/stac-hojd/v1"
+
+# Collection ids as published in the live catalogues. Buildings are one item
+# per municipality whose asset is a ZIP holding a GeoPackage; "Laserdata Skog"
+# is published as `dsm-skoglig-copc` — a surface-model point cloud in
+# LAZ/COPC — on the elevation root.
+COLLECTION_BUILDINGS = "byggnader"
+COLLECTION_LIDAR = "dsm-skoglig-copc"
+
+# Media types, so an asset is chosen by *what it is* rather than by hoping the
+# publisher named the key "data". Byggnader delivers a zipped GeoPackage;
+# Laserdata Skog delivers LAZ organised as COPC (Cloud Optimized Point Cloud).
+MEDIA_GEOPACKAGE = "application/geopackage+sqlite3"
+MEDIA_COPC = "application/vnd.laszip+copc"
+MEDIA_LAZ = "application/vnd.laszip"
+MEDIA_LAS = "application/vnd.las"
+MEDIA_GEOJSON = "application/geo+json"
+MEDIA_ZIP = "application/zip"
+
+# Longest suffix first: a COPC file is also a .laz, and reading it as a plain
+# one would download the whole tile instead of the part we asked for.
+_EXTENSION_MEDIA: tuple[tuple[str, str], ...] = (
+ (".copc.laz", MEDIA_COPC),
+ (".gpkg", MEDIA_GEOPACKAGE),
+ (".geojson", MEDIA_GEOJSON),
+ (".laz", MEDIA_LAZ),
+ (".las", MEDIA_LAS),
+ (".zip", MEDIA_ZIP),
+)
+
+
+class GeotorgetError(RuntimeError):
+ """Any failure talking to Geotorget."""
+
+
+class MissingCredentials(GeotorgetError):
+ """No usable credentials were supplied."""
+
+
+@dataclasses.dataclass(frozen=True)
+class Credentials:
+ username: str
+ password: str
+
+ def validate(self) -> None:
+ if not self.username or not self.password:
+ raise MissingCredentials(
+ "a STAC username and password are both required; for Lantmäteriet, "
+ "order access at https://geotorget.lantmateriet.se and set "
+ "roofmodel.stac_username and roofmodel.stac_password to your "
+ "Geotorget account credentials"
+ )
+
+
+def media_type_for(href: str) -> str | None:
+ """Media type implied by a URL's extension, or None if it says nothing."""
+ path = href.split("?", 1)[0].split("#", 1)[0].lower()
+ for suffix, media in _EXTENSION_MEDIA:
+ if path.endswith(suffix):
+ return media
+ return None
+
+
+@dataclasses.dataclass(frozen=True)
+class Asset:
+ """One STAC asset: where it is, and what it is."""
+
+ href: str
+ media_type: str | None = None
+ roles: tuple[str, ...] = ()
+ title: str = ""
+
+ @property
+ def effective_media_type(self) -> str | None:
+ """Declared media type, or the one the extension implies.
+
+ Catalogues are inconsistent about `type`, and an asset with no declared
+ type is common enough that refusing to guess would mean refusing most
+ real items. The extension is only consulted when nothing was declared.
+ """
+ return self.media_type or media_type_for(self.href)
+
+
+@dataclasses.dataclass
+class StacItem:
+ """One STAC item, reduced to what the pipeline needs."""
+
+ item_id: str
+ collection: str
+ assets: dict[str, Asset]
+ captured_at: dt.datetime | None
+ raw: dict[str, Any] = dataclasses.field(default_factory=dict, repr=False)
+
+ def __post_init__(self) -> None:
+ # A bare href is accepted wherever an Asset is, so callers and tests can
+ # write {"data": "http://.../tile.copc.laz"} without losing the typing
+ # that selection depends on -- the extension supplies it.
+ self.assets = {
+ name: value if isinstance(value, Asset) else Asset(href=str(value))
+ for name, value in (self.assets or {}).items()
+ }
+
+ def pick(self, *media_types: str) -> Asset | None:
+ """Best asset for a wanted media type, most preferred type first.
+
+ Where nothing declares a usable type the search widens: an asset with
+ the `data` role, then a lone asset, since an item carrying exactly one
+ asset is unambiguous however it is labelled.
+
+ Both fallbacks consider only assets of *unknown* type. Guessing in the
+ absence of information is reasonable; guessing against it is not, and an
+ item whose single asset is a thumbnail must not be handed back as a
+ point cloud.
+ """
+ for wanted in media_types:
+ for asset in self.assets.values():
+ if asset.effective_media_type == wanted:
+ return asset
+ untyped = [a for a in self.assets.values() if a.effective_media_type is None]
+ for asset in untyped:
+ if "data" in asset.roles:
+ return asset
+ if len(untyped) == 1:
+ return untyped[0]
+ return None
+
+ def asset_url(self, *preferred: str) -> str | None:
+ """First matching asset href, trying each preferred key in order."""
+ for key in preferred:
+ if key in self.assets:
+ return self.assets[key].href
+ first = next(iter(self.assets.values()), None)
+ return first.href if first else None
+
+
+def _parse_datetime(value: str | None) -> dt.datetime | None:
+ """Parse a STAC RFC 3339 timestamp.
+
+ Lantmaeteriet is backfilling `properties.datetime` across 2026, so it is
+ routinely absent. That is a missing provenance date, not an error -- the UI
+ degrades to "capture date unknown" rather than refusing the model.
+ """
+ if not value:
+ return None
+ try:
+ return dt.datetime.fromisoformat(value.replace("Z", "+00:00"))
+ except ValueError:
+ return None
+
+
+def _item_from_feature(feature: dict[str, Any]) -> StacItem:
+ props = feature.get("properties") or {}
+ assets = {
+ name: Asset(
+ href=asset.get("href", ""),
+ media_type=asset.get("type") or None,
+ roles=tuple(asset.get("roles") or ()),
+ title=str(asset.get("title") or ""),
+ )
+ for name, asset in (feature.get("assets") or {}).items()
+ if asset.get("href")
+ }
+ captured = _parse_datetime(props.get("datetime")) or _parse_datetime(
+ props.get("start_datetime")
+ )
+ if captured is None:
+ # Laser strips carry an acquisition date as `datum` (e.g. "20180301")
+ # even where the STAC datetime has not been backfilled yet.
+ datum = props.get("datum")
+ if datum:
+ try:
+ captured = dt.datetime.strptime(str(datum), "%Y%m%d").replace(
+ tzinfo=dt.timezone.utc
+ )
+ except ValueError:
+ captured = None
+ return StacItem(
+ item_id=feature.get("id", ""),
+ collection=feature.get("collection", ""),
+ assets=assets,
+ captured_at=captured,
+ raw=feature,
+ )
+
+
+class StacClient:
+ """Thin client for a STAC search-and-download API over HTTP Basic auth.
+
+ Defaults target Lantmaeteriet's Geotorget catalog, but nothing here
+ depends on it: base_url and the collection ids passed to `search` are the
+ whole coupling.
+ """
+
+ def __init__(
+ self,
+ credentials: Credentials,
+ session: Any = None,
+ base_url: str = DEFAULT_BASE_URL,
+ timeout: float = 60.0,
+ ) -> None:
+ self._base_url = base_url.rstrip("/")
+ # Lantmäteriet always demands the operator's Geotorget account, so
+ # missing credentials there deserve the ordering instructions early.
+ # A custom catalog may be open: no credentials means anonymous access,
+ # while half a credential is still an error either way.
+ self._anonymous = not (credentials.username or credentials.password)
+ if not self._anonymous or self._base_url in (DEFAULT_BASE_URL, DEFAULT_LIDAR_BASE_URL):
+ credentials.validate()
+ self._credentials = credentials
+ self._timeout = timeout
+ if session is None:
+ import requests # imported lazily so tests can inject a fake session
+
+ session = requests.Session()
+ if not self._anonymous:
+ session.auth = (credentials.username, credentials.password)
+ self._session = session
+
+ @property
+ def session(self) -> Any:
+ """The authenticated session, for readers that stream their own ranges."""
+ return self._session
+
+ def search(
+ self,
+ collection: str,
+ bbox: tuple[float, float, float, float],
+ limit: int = 20,
+ ) -> list[StacItem]:
+ """POST {base}/search for one collection over a bbox.
+
+ The bbox is (min_x, min_y, max_x, max_y) in whatever CRS the catalog
+ expects -- the STAC spec says WGS84 lon/lat, and the live Lantmaeteriet
+ service follows it -- so the caller chooses what to send and no
+ reprojection happens here.
+ """
+ body = {
+ "collections": [collection],
+ "bbox": list(bbox),
+ "limit": limit,
+ }
+ url = f"{self._base_url}/search"
+ try:
+ resp = self._session.post(url, json=body, timeout=self._timeout)
+ except Exception as exc: # network, DNS, TLS
+ raise GeotorgetError(f"STAC search failed: {exc}") from exc
+ if resp.status_code in (401, 403):
+ if self._anonymous:
+ raise MissingCredentials(
+ f"the STAC catalog requires credentials for {collection} "
+ f"(HTTP {resp.status_code}). Set roofmodel.stac_username "
+ "and roofmodel.stac_password for this catalog."
+ )
+ raise MissingCredentials(
+ f"the STAC catalog rejected the credentials for {collection} "
+ f"(HTTP {resp.status_code}). Check the username and password, and "
+ "that the account has ordered access to this product."
+ )
+ if resp.status_code != 200:
+ raise GeotorgetError(f"STAC search returned HTTP {resp.status_code}")
+ payload = resp.json()
+ return [_item_from_feature(f) for f in payload.get("features", [])]
+
+ def download(self, url: str) -> bytes:
+ """Fetch one asset."""
+ try:
+ resp = self._session.get(url, timeout=self._timeout)
+ except Exception as exc:
+ raise GeotorgetError(f"asset download failed: {exc}") from exc
+ if resp.status_code != 200:
+ raise GeotorgetError(f"asset download returned HTTP {resp.status_code}")
+ return resp.content
+
+
+def newest_capture(items: Iterable[StacItem]) -> dt.datetime | None:
+ """Most recent known capture date across items, or None if none carry one."""
+ dates = [i.captured_at for i in items if i.captured_at is not None]
+ return max(dates) if dates else None
+
+
+# The client predates its generalization; the old name stays importable.
+GeotorgetClient = StacClient
diff --git a/roofmodel/ftw_roofmodel/pipeline.py b/roofmodel/ftw_roofmodel/pipeline.py
new file mode 100644
index 000000000..f8f764e3c
--- /dev/null
+++ b/roofmodel/ftw_roofmodel/pipeline.py
@@ -0,0 +1,310 @@
+"""Derive roof geometry for a site: STAC search -> LiDAR -> planes -> arrays.
+
+The output contract is a versioned `roof_model.json`, which is the whole reason
+this lives in a separate module rather than inside core: core only ever reads
+that document, so the heavy geospatial dependencies, their failure modes and
+their update cadence stay on this side of the boundary.
+
+Nothing here is authoritative. The derived arrays *pre-fill* the operator's
+editable `weather.pv_arrays`; they are hints, and the numeric editor stays the
+final word. A failure produces a clean error and leaves the existing config
+untouched.
+"""
+
+from __future__ import annotations
+
+import dataclasses
+import datetime as dt
+from typing import Any
+
+from . import geotorget, pointcloud, sweref
+from .buildings import (
+ DEFAULT_EAVES_BUFFER_M,
+ DEFAULT_SEARCH_RADIUS_M,
+ Building,
+ building_from_drawn_footprint,
+ clip_to_footprint,
+ search_buildings,
+)
+from .geotorget import (
+ COLLECTION_BUILDINGS,
+ COLLECTION_LIDAR,
+ Credentials,
+ GeotorgetClient,
+ GeotorgetError,
+ StacItem,
+ newest_capture,
+)
+from .segment import (
+ DEFAULT_MODULE_W_PER_M2,
+ DEFAULT_PACKING_FACTOR,
+ RoofPlane,
+ segment_roof,
+)
+
+SCHEMA_VERSION = 1
+
+# How far around the site to pull LiDAR. 40 m comfortably contains a detached
+# house and its outbuildings without dragging in the neighbours' roofs.
+DEFAULT_RADIUS_M = 40.0
+
+# Roof faces smaller than this are dormers, porches and sheds: real surfaces,
+# but not worth proposing as a PV array.
+MIN_ARRAY_AREA_M2 = 8.0
+
+# Lantmaeteriet's Laserdata Skog is specified at 1-2 points/m2.
+NOMINAL_POINT_DENSITY = 1.5
+
+# Below this, a clipped footprint cannot support a plane fit -- segment_roof
+# needs 40 points for a single face, and a roof has at least two.
+MIN_POINTS_AFTER_CLIP = 80
+
+
+class RoofModelError(RuntimeError):
+ """Derivation failed."""
+
+
+@dataclasses.dataclass
+class DerivedArray:
+ name: str
+ rated_w: float
+ tilt_deg: float
+ azimuth_deg: float
+ area_m2: float
+ segment_id: str
+
+ def to_json(self) -> dict[str, Any]:
+ return {
+ "name": self.name,
+ "rated_w": round(self.rated_w),
+ "tilt_deg": round(self.tilt_deg, 1),
+ "azimuth_deg": round(self.azimuth_deg, 1),
+ "area_m2": round(self.area_m2, 1),
+ "segment_id": self.segment_id,
+ }
+
+
+def _compass_name(azimuth_deg: float, tilt_deg: float) -> str:
+ """Human-readable face name, e.g. "Roof south"."""
+ if tilt_deg < 5.0:
+ return "Roof flat"
+ points = [
+ (0, "north"), (45, "north-east"), (90, "east"), (135, "south-east"),
+ (180, "south"), (225, "south-west"), (270, "west"), (315, "north-west"),
+ (360, "north"),
+ ]
+ best = min(points, key=lambda p: abs(p[0] - azimuth_deg))
+ return f"Roof {best[1]}"
+
+
+def planes_to_arrays(
+ planes: list[RoofPlane],
+ *,
+ packing_factor: float = DEFAULT_PACKING_FACTOR,
+ module_w_per_m2: float = DEFAULT_MODULE_W_PER_M2,
+ min_area_m2: float = MIN_ARRAY_AREA_M2,
+) -> list[DerivedArray]:
+ """Convert roof planes into candidate PV arrays.
+
+ North-facing pitched roofs are dropped: at Swedish latitudes a north face
+ at any real pitch yields so little that proposing it as an array would be
+ noise in the operator's config. Flat roofs are kept -- they are mounted to
+ face south regardless of which way the building points.
+ """
+ arrays: list[DerivedArray] = []
+ used: dict[str, int] = {}
+ for idx, plane in enumerate(planes):
+ if plane.area_m2 < min_area_m2:
+ continue
+ if plane.tilt_deg >= 5.0 and (plane.azimuth_deg <= 45.0 or plane.azimuth_deg >= 315.0):
+ continue
+ name = _compass_name(plane.azimuth_deg, plane.tilt_deg)
+ used[name] = used.get(name, 0) + 1
+ if used[name] > 1:
+ name = f"{name} {used[name]}"
+ arrays.append(
+ DerivedArray(
+ name=name,
+ rated_w=plane.rated_w(packing_factor, module_w_per_m2),
+ tilt_deg=plane.tilt_deg,
+ azimuth_deg=plane.azimuth_deg,
+ area_m2=plane.area_m2,
+ segment_id=f"seg-{idx}",
+ )
+ )
+ return arrays
+
+
+def load_points(data: bytes) -> Any:
+ """Decode a whole LAZ/LAS payload into (N, 3) SWEREF 99 TM metres.
+
+ Re-exported from pointcloud so that laspy stays lazily imported: everything
+ above -- projection, segmentation, array derivation -- is importable and
+ testable without the geospatial stack installed.
+ """
+ return pointcloud.load_points(data)
+
+
+def _read_lidar(
+ client: GeotorgetClient,
+ items: list[StacItem],
+ chosen: Building | None,
+) -> tuple[Any, str]:
+ """Points for the first readable LiDAR asset, and how they were fetched.
+
+ Laserdata Skog is LAZ organised as COPC, so when the operator has already
+ picked a building there is no reason to move the rest of a 2.5 km tile
+ across the network: the footprint's bounding box is exactly the query COPC
+ is built to answer. Everything about that is best-effort -- a plain LAZ
+ asset, a host that ignores `Range`, or a laspy without COPC support all
+ fall back to reading the tile whole.
+ """
+ for item in items:
+ asset = item.pick(
+ geotorget.MEDIA_COPC, geotorget.MEDIA_LAZ, geotorget.MEDIA_LAS
+ )
+ if asset is None or not asset.href:
+ continue
+ if chosen is not None and asset.effective_media_type == geotorget.MEDIA_COPC:
+ session = getattr(client, "session", None)
+ if session is not None:
+ bounds = pointcloud.bounds_of(chosen.ring_sweref, DEFAULT_EAVES_BUFFER_M)
+ try:
+ return pointcloud.read_copc_window(session, asset.href, bounds), "copc-window"
+ except pointcloud.PointCloudError:
+ # A slow success beats a failure; the operator gets their
+ # roof either way, and `fetch` records which path ran.
+ pass
+ return load_points(client.download(asset.href)), "whole-tile"
+ raise RoofModelError("LiDAR tiles carried no readable point data")
+
+
+def derive(
+ *,
+ latitude: float,
+ longitude: float,
+ credentials: Credentials,
+ client: GeotorgetClient | None = None,
+ radius_m: float = DEFAULT_RADIUS_M,
+ packing_factor: float = DEFAULT_PACKING_FACTOR,
+ module_w_per_m2: float = DEFAULT_MODULE_W_PER_M2,
+ building_id: str | None = None,
+ footprint: list[Any] | None = None,
+ now: dt.datetime | None = None,
+ base_url: str = geotorget.DEFAULT_BASE_URL,
+ buildings_collection: str = COLLECTION_BUILDINGS,
+ lidar_collection: str = COLLECTION_LIDAR,
+ bbox_epsg: int = 4326,
+) -> dict[str, Any]:
+ """Derive a roof model for one site and return it as a JSON-ready dict.
+
+ Pass `building_id` -- one the operator picked from `search_buildings` -- to
+ clip the LiDAR to that footprint before segmenting, or `footprint` -- a
+ hand-drawn [lon, lat] ring -- where the catalog publishes no building
+ dataset to pick from. A drawn footprint wins over a building id. Without
+ either the whole radius is segmented, which will happily return the
+ neighbour's roof and lets coplanar buildings steal each other's points;
+ see buildings.py.
+
+ The catalog parameters default to Lantmaeteriet; any STAC-conformant
+ catalog can stand in (see geotorget.py), as long as its data arrives in
+ SWEREF 99 TM -- the segmentation works in that frame.
+ """
+ # Lantmaeteriet splits its STAC service per product family: buildings live
+ # on the vector root, the point clouds on the elevation root. An injected
+ # client (tests, the demo) serves both; so does a custom single-root
+ # catalog. Only the Lantmaeteriet default needs the second client.
+ lidar_client = client
+ if client is None:
+ client = GeotorgetClient(credentials, base_url=base_url)
+ lidar_client = client
+ if base_url == geotorget.DEFAULT_BASE_URL:
+ lidar_client = GeotorgetClient(
+ credentials, base_url=geotorget.DEFAULT_LIDAR_BASE_URL
+ )
+
+ chosen: Building | None = None
+ if footprint:
+ chosen = building_from_drawn_footprint(footprint)
+ elif building_id:
+ # Re-find the picked footprint with the same reach the picker's own
+ # search had. `radius_m` is the LiDAR/segmentation radius (tens of
+ # metres); a building the operator picked from the 150 m search would
+ # vanish here if it alone bounded the re-search.
+ candidates = search_buildings(
+ client, latitude=latitude, longitude=longitude,
+ radius_m=max(radius_m, DEFAULT_SEARCH_RADIUS_M),
+ collection=buildings_collection, bbox_epsg=bbox_epsg,
+ )
+ chosen = next((b for b in candidates if b.building_id == building_id), None)
+ if chosen is None:
+ raise RoofModelError(
+ f"building {building_id!r} was not found near this site; it may "
+ "have been picked against a different coordinate"
+ )
+
+ # The LiDAR lookup centres on the roof being derived: for a picked
+ # building that is its centroid, not the site pin — a barn at the edge of
+ # the search radius must find *its* tile, not the pin's.
+ lidar_lat, lidar_lon = latitude, longitude
+ if chosen is not None:
+ lidar_lat, lidar_lon = chosen.centroid_wgs84()
+ bbox = sweref.stac_search_bbox(lidar_lat, lidar_lon, radius_m, bbox_epsg)
+
+ try:
+ lidar_items: list[StacItem] = lidar_client.search(lidar_collection, bbox)
+ except GeotorgetError:
+ raise
+ if not lidar_items:
+ hint = "; Lantmaeteriet data is Sweden only" if lidar_collection == COLLECTION_LIDAR else ""
+ raise RoofModelError(
+ f"no LiDAR tiles cover ({lidar_lat:.5f}, {lidar_lon:.5f}) in "
+ f"collection {lidar_collection!r}{hint}"
+ )
+
+ points, fetch = _read_lidar(lidar_client, lidar_items, chosen)
+ if points is None or len(points) == 0:
+ raise RoofModelError("LiDAR tiles carried no readable point data")
+
+ total_returns = len(points)
+ if chosen is not None:
+ points = clip_to_footprint(points, chosen.ring_sweref)
+ if len(points) < MIN_POINTS_AFTER_CLIP:
+ raise RoofModelError(
+ f"only {len(points)} LiDAR returns fall on building "
+ f"{chosen.building_id!r}. The footprint and the point cloud may "
+ "be from different years, or the building is newer than the scan."
+ )
+
+ planes = segment_roof(points, point_density=NOMINAL_POINT_DENSITY)
+ arrays = planes_to_arrays(
+ planes, packing_factor=packing_factor, module_w_per_m2=module_w_per_m2
+ )
+
+ captured = newest_capture(lidar_items)
+ stamp = now or dt.datetime.now(dt.timezone.utc)
+ return {
+ "schema_version": SCHEMA_VERSION,
+ "site": {"latitude": latitude, "longitude": longitude, "radius_m": radius_m},
+ "source": {
+ "provider": "lantmateriet" if base_url == geotorget.DEFAULT_BASE_URL else base_url,
+ "collection": lidar_collection,
+ "item_count": len(lidar_items),
+ "dataset_datetime": captured.isoformat() if captured else None,
+ # "copc-window" means only the footprint's neighbourhood was moved
+ # across the network, which also makes returns_in_radius below a
+ # count over that window rather than over the whole search radius.
+ "fetch": fetch,
+ },
+ "building": {
+ "building_id": chosen.building_id,
+ "area_m2": round(chosen.area_m2, 1),
+ "footprint": chosen.to_geojson()["geometry"],
+ "returns_used": len(points),
+ "returns_in_radius": total_returns,
+ } if chosen is not None else None,
+ "arrays": [a.to_json() for a in arrays],
+ "planes_found": len(planes),
+ "captured_at_ms": int(captured.timestamp() * 1000) if captured else None,
+ "derived_at_ms": int(stamp.timestamp() * 1000),
+ }
diff --git a/roofmodel/ftw_roofmodel/pointcloud.py b/roofmodel/ftw_roofmodel/pointcloud.py
new file mode 100644
index 000000000..c49927dea
--- /dev/null
+++ b/roofmodel/ftw_roofmodel/pointcloud.py
@@ -0,0 +1,245 @@
+"""Read Lantmaeteriet LiDAR, fetching only the part of the tile we need.
+
+*Laserdata Nedladdning, Skog* is delivered as LAZ organised as **COPC** (Cloud
+Optimized Point Cloud): the points are ordered into an octree and the node index
+lives in a VLR at a known offset, so a reader that can issue HTTP range requests
+can pull the handful of octree nodes covering one building instead of the whole
+2.5 km tile.
+
+That is worth real money on a Pi. A Laserdata Skog tile is hundreds of megabytes;
+a detached house is a few tens of metres across. Since the operator has already
+told us *which building* they mean, the bounding box of that footprint is exactly
+the query COPC exists to answer.
+
+The fallbacks are deliberate and ordered, because none of the preconditions are
+guaranteed:
+
+ 1. COPC asset + a bounding box + a server that honours `Range` -> spatial query.
+ 2. Anything else -> download the asset and read it whole.
+
+A server that ignores `Range` returns 200 with the entire body, which would
+otherwise be mistaken for a successful partial read, so that case is detected on
+the status code rather than assumed away.
+"""
+
+from __future__ import annotations
+
+import io
+from typing import Any
+
+__all__ = [
+ "PointCloudError",
+ "HttpRangeFile",
+ "bounds_of",
+ "copc_query_z_range",
+ "load_points",
+ "read_copc_window",
+]
+
+# Read this much per range request. COPC chunks are small, and a request per
+# chunk would spend more time in round trips than in transfer.
+DEFAULT_CHUNK_BYTES = 1 << 20
+
+
+class PointCloudError(RuntimeError):
+ """The LiDAR asset could not be read."""
+
+
+class HttpRangeFile(io.RawIOBase):
+ """A seekable read-only file over HTTP `Range` requests.
+
+ laspy's COPC reader needs `seek`/`read` and nothing else, so this is the
+ whole adapter: it turns an HTTP URL into something that behaves like an open
+ file without ever holding the tile in memory.
+ """
+
+ def __init__(self, session: Any, url: str, *, timeout: float = 60.0,
+ chunk_bytes: int = DEFAULT_CHUNK_BYTES) -> None:
+ self._session = session
+ self._url = url
+ self._timeout = timeout
+ self._chunk = max(1, chunk_bytes)
+ self._pos = 0
+ self._size: int | None = None
+ # One cached chunk. COPC reads are clustered -- header, then index, then
+ # the nodes -- so a single block absorbs most of the repeat traffic.
+ self._cache: tuple[int, bytes] | None = None
+ self.requests = 0
+ self.bytes_fetched = 0
+
+ # -- io.RawIOBase ----------------------------------------------------
+ def readable(self) -> bool:
+ return True
+
+ def seekable(self) -> bool:
+ return True
+
+ def tell(self) -> int:
+ return self._pos
+
+ def seek(self, offset: int, whence: int = io.SEEK_SET) -> int:
+ if whence == io.SEEK_SET:
+ self._pos = offset
+ elif whence == io.SEEK_CUR:
+ self._pos += offset
+ elif whence == io.SEEK_END:
+ self._pos = self.size + offset
+ else: # pragma: no cover - io module only defines the three
+ raise ValueError(f"invalid whence {whence}")
+ self._pos = max(0, self._pos)
+ return self._pos
+
+ def read(self, size: int = -1) -> bytes:
+ if size is None or size < 0:
+ size = max(0, self.size - self._pos)
+ if size == 0:
+ return b""
+ data = self._read_at(self._pos, size)
+ self._pos += len(data)
+ return data
+
+ def readall(self) -> bytes:
+ return self.read(-1)
+
+ def readinto(self, buffer) -> int: # type: ignore[override]
+ data = self.read(len(buffer))
+ buffer[: len(data)] = data
+ return len(data)
+
+ # -- range plumbing --------------------------------------------------
+ @property
+ def size(self) -> int:
+ if self._size is None:
+ self._size = self._head_size()
+ return self._size
+
+ def _head_size(self) -> int:
+ try:
+ resp = self._session.head(self._url, timeout=self._timeout)
+ except Exception as exc:
+ raise PointCloudError(f"could not stat {self._url}: {exc}") from exc
+ length = (getattr(resp, "headers", None) or {}).get("Content-Length")
+ if getattr(resp, "status_code", 0) != 200 or not length:
+ raise PointCloudError(
+ "the LiDAR host did not report a size, so it cannot be read in ranges"
+ )
+ return int(length)
+
+ def _read_at(self, offset: int, size: int) -> bytes:
+ cached = self._from_cache(offset, size)
+ if cached is not None:
+ return cached
+ want = max(size, self._chunk)
+ end = offset + want - 1
+ try:
+ resp = self._session.get(
+ self._url,
+ headers={"Range": f"bytes={offset}-{end}"},
+ timeout=self._timeout,
+ )
+ except Exception as exc:
+ raise PointCloudError(f"range request failed: {exc}") from exc
+ status = getattr(resp, "status_code", 0)
+ if status == 200:
+ # The server ignored Range and sent everything. Honest failure: the
+ # caller falls back to a whole-tile read rather than silently
+ # paying for the full download on every seek.
+ raise PointCloudError("the LiDAR host does not support range requests")
+ if status != 206:
+ raise PointCloudError(f"range request returned HTTP {status}")
+ body = resp.content
+ self.requests += 1
+ self.bytes_fetched += len(body)
+ self._cache = (offset, body)
+ return body[:size]
+
+ def _from_cache(self, offset: int, size: int) -> bytes | None:
+ if self._cache is None:
+ return None
+ start, body = self._cache
+ if offset < start or offset + size > start + len(body):
+ return None
+ rel = offset - start
+ return body[rel : rel + size]
+
+
+def bounds_of(ring: list[tuple[float, float]], pad_m: float = 2.0) -> tuple[float, float, float, float]:
+ xs = [p[0] for p in ring]
+ ys = [p[1] for p in ring]
+ return (min(xs) - pad_m, min(ys) - pad_m, max(xs) + pad_m, max(ys) + pad_m)
+
+
+def _points_from_las(las: Any) -> Any:
+ import numpy as np
+
+ return np.column_stack([np.asarray(las.x), np.asarray(las.y), np.asarray(las.z)])
+
+
+def load_points(data: bytes) -> Any:
+ """Decode a whole LAZ/LAS payload into (N, 3) SWEREF 99 TM metres."""
+ laspy = _import_laspy()
+ with laspy.open(io.BytesIO(data)) as reader:
+ las = reader.read()
+ return _points_from_las(las)
+
+
+def _import_laspy():
+ try:
+ import laspy
+ except ImportError as exc: # pragma: no cover - depends on the install
+ raise PointCloudError(
+ "laspy is required to read Lantmaeteriet LiDAR. Install the module's "
+ "extras: pip install -e roofmodel[geo]"
+ ) from exc
+ return laspy
+
+
+def copc_query_z_range(center_z: float, halfsize: float) -> tuple[float, float]:
+ """The vertical range a windowed COPC query must span: the whole cube.
+
+ Lantmäteriet's COPC writer emits octree z-keys measured from some other
+ origin than the file's own cube (observed live on Laserdata Skog: a level-6
+ node keyed z=19, implying a -1698..-1542 m slab, holds points at +18..+42 m
+ while its x/y keys are exact). A 2D query lets laspy fill z from the header
+ and prune nodes by those broken slabs, which silently discards every dense
+ deep level and leaves only the sparse preview points. Spanning the full
+ cube keeps z from ever pruning; x/y pruning and the exact post-filter still
+ bound the read.
+ """
+ pad = abs(halfsize) + 1.0
+ return (center_z - 2.0 * pad, center_z + 2.0 * pad)
+
+
+def read_copc_window(
+ session: Any,
+ url: str,
+ bounds: tuple[float, float, float, float],
+ *,
+ timeout: float = 60.0,
+) -> Any:
+ """Points inside `bounds` from a COPC file, over HTTP range requests.
+
+ Raises PointCloudError if the file or the host cannot support it, so the
+ caller can fall back to a whole-tile read.
+ """
+ laspy = _import_laspy()
+ try:
+ from laspy.copc import Bounds, CopcReader
+ except ImportError as exc:
+ raise PointCloudError(
+ "this laspy build has no COPC support; install laspy[lazrs] 2.5 or newer"
+ ) from exc
+
+ min_x, min_y, max_x, max_y = bounds
+ handle = HttpRangeFile(session, url, timeout=timeout)
+ try:
+ with CopcReader.open(handle) as reader:
+ info = reader.copc_info
+ z_lo, z_hi = copc_query_z_range(float(info.center[2]), float(info.halfsize))
+ query = Bounds(mins=[min_x, min_y, z_lo], maxs=[max_x, max_y, z_hi])
+ points = reader.query(query)
+ except PointCloudError:
+ raise
+ except Exception as exc:
+ raise PointCloudError(f"COPC read failed: {exc}") from exc
+ return _points_from_las(points)
diff --git a/roofmodel/ftw_roofmodel/segment.py b/roofmodel/ftw_roofmodel/segment.py
new file mode 100644
index 000000000..5a2bc0a11
--- /dev/null
+++ b/roofmodel/ftw_roofmodel/segment.py
@@ -0,0 +1,251 @@
+"""Roof-plane segmentation from a LiDAR point cloud.
+
+Implements the method described in the SPAN paper (Yavuzdogan, *Renewable
+Energy* 2023, doi:10.1016/j.renene.2023.119022): iterative RANSAC plane fitting
+to pull one roof surface at a time out of the cloud, then DBSCAN over each
+plane's inliers to split faces that share a plane but not a location -- the two
+halves of a gable on opposite wings of a building fit the same equation and are
+not the same roof face.
+
+The method is reimplemented from its description. No code is taken from SPAN's
+QGIS plugin, which is GPL; everything here uses numpy and scikit-learn, both
+BSD, so this module carries no copyleft obligation.
+
+Coordinates are SWEREF 99 TM metres as (easting, northing, height). Azimuth
+follows FTW's convention: 0 = north, 90 = east, 180 = south, 270 = west.
+"""
+
+from __future__ import annotations
+
+import dataclasses
+import math
+
+import numpy as np
+
+# A roof plane must hold at least this many returns to be believed. Below this
+# a "plane" is usually a chimney, an aerial, or three points of noise that
+# happen to be collinear.
+MIN_PLANE_POINTS = 40
+
+# Surfaces flatter than this are treated as flat roofs: their azimuth is
+# meaningless (the normal is essentially vertical, so its horizontal component
+# is numerical noise that can point anywhere).
+FLAT_TILT_DEG = 5.0
+
+# Roofs steeper than this are walls, dormer cheeks or mis-fits.
+MAX_TILT_DEG = 80.0
+
+# Fraction of a roof plane's area that can carry modules once you subtract
+# ridges, eaves, chimneys, vents and walkways.
+DEFAULT_PACKING_FACTOR = 0.70
+
+# Module DC rating per square metre of module. ~20% efficiency at 1000 W/m2.
+DEFAULT_MODULE_W_PER_M2 = 200.0
+
+
+@dataclasses.dataclass
+class RoofPlane:
+ """One contiguous roof surface."""
+
+ tilt_deg: float
+ azimuth_deg: float
+ area_m2: float
+ point_count: int
+ mean_height_m: float
+
+ def rated_w(
+ self,
+ packing_factor: float = DEFAULT_PACKING_FACTOR,
+ module_w_per_m2: float = DEFAULT_MODULE_W_PER_M2,
+ ) -> float:
+ """Installable DC capacity for this surface, in watts."""
+ return self.area_m2 * packing_factor * module_w_per_m2
+
+ def kwp(
+ self,
+ packing_factor: float = DEFAULT_PACKING_FACTOR,
+ module_w_per_m2: float = DEFAULT_MODULE_W_PER_M2,
+ ) -> float:
+ """Installable DC capacity for this surface, in kWp (test helper)."""
+ return self.rated_w(packing_factor, module_w_per_m2) / 1000.0
+
+
+def _fit_plane(points: np.ndarray) -> np.ndarray:
+ """Least-squares plane through points; returns a unit normal pointing up.
+
+ Uses the smallest singular vector of the mean-centred points, which is the
+ total-least-squares fit -- it minimises perpendicular distance rather than
+ vertical distance, so a steep roof is not biased the way an ordinary
+ z = ax + by + c regression would bias it.
+ """
+ centred = points - points.mean(axis=0)
+ _, _, vh = np.linalg.svd(centred, full_matrices=False)
+ normal = vh[-1]
+ if normal[2] < 0:
+ normal = -normal
+ return normal / np.linalg.norm(normal)
+
+
+def _tilt_azimuth(normal: np.ndarray) -> tuple[float, float]:
+ """Convert an upward unit normal to (tilt_deg, azimuth_deg)."""
+ tilt = math.degrees(math.acos(max(-1.0, min(1.0, float(normal[2])))))
+ if tilt < FLAT_TILT_DEG:
+ # Horizontal component is noise at this point; report due south, which
+ # is what a flat-roof array is normally mounted to face anyway.
+ return tilt, 180.0
+ east, north = float(normal[0]), float(normal[1])
+ azimuth = math.degrees(math.atan2(east, north)) % 360.0
+ return tilt, azimuth
+
+
+def _convex_hull_area(xy: np.ndarray) -> float:
+ """Area of the convex hull of 2D points (monotone chain, shoelace).
+
+ Hand-rolled to avoid a scipy dependency for one small routine. The hull
+ overestimates a concave roof outline, so it is only used as a fallback
+ when the point count is too low for the density estimate to be stable.
+ """
+ pts = np.unique(xy, axis=0)
+ if len(pts) < 3:
+ return 0.0
+ order = np.lexsort((pts[:, 1], pts[:, 0]))
+ pts = pts[order]
+
+ def cross(o, a, b):
+ return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])
+
+ lower: list = []
+ for p in pts:
+ while len(lower) >= 2 and cross(lower[-2], lower[-1], p) <= 0:
+ lower.pop()
+ lower.append(p)
+ upper: list = []
+ for p in pts[::-1]:
+ while len(upper) >= 2 and cross(upper[-2], upper[-1], p) <= 0:
+ upper.pop()
+ upper.append(p)
+ hull = np.array(lower[:-1] + upper[:-1])
+ if len(hull) < 3:
+ return 0.0
+ x, y = hull[:, 0], hull[:, 1]
+ return 0.5 * abs(float(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1))))
+
+
+def _surface_area(points: np.ndarray, tilt_deg: float, point_density: float | None) -> float:
+ """Sloped surface area of a roof face, in m2.
+
+ LiDAR density is quoted per square metre of *ground*, so a known density
+ gives the horizontal footprint directly from the point count; dividing by
+ cos(tilt) lifts that onto the slope. Without a density we fall back to the
+ convex hull of the horizontal projection, which is looser -- it fills in
+ L-shapes and courtyards.
+ """
+ if point_density and point_density > 0:
+ horizontal = len(points) / point_density
+ else:
+ horizontal = _convex_hull_area(points[:, :2])
+ cos_t = math.cos(math.radians(min(tilt_deg, MAX_TILT_DEG)))
+ if cos_t <= 1e-6:
+ return horizontal
+ return horizontal / cos_t
+
+
+def _ransac_plane(
+ points: np.ndarray,
+ threshold_m: float,
+ iterations: int,
+ rng: np.random.Generator,
+) -> np.ndarray | None:
+ """Return a boolean inlier mask for the best plane found, or None.
+
+ Plain RANSAC over point triples. scikit-learn's RANSACRegressor is not used
+ because it regresses z on (x, y) and so cannot represent a vertical or
+ near-vertical surface, and weights errors vertically rather than
+ perpendicular to the plane.
+ """
+ n = len(points)
+ if n < 3:
+ return None
+ best_mask = None
+ best_count = 0
+ for _ in range(iterations):
+ idx = rng.choice(n, size=3, replace=False)
+ a, b, c = points[idx]
+ normal = np.cross(b - a, c - a)
+ norm = np.linalg.norm(normal)
+ if norm < 1e-9:
+ continue # degenerate (collinear) sample
+ normal = normal / norm
+ distances = np.abs((points - a) @ normal)
+ mask = distances < threshold_m
+ count = int(mask.sum())
+ if count > best_count:
+ best_count, best_mask = count, mask
+ if best_mask is None or best_count < 3:
+ return None
+ return best_mask
+
+
+def segment_roof(
+ points: np.ndarray,
+ *,
+ threshold_m: float = 0.25,
+ max_planes: int = 8,
+ min_plane_points: int = MIN_PLANE_POINTS,
+ cluster_eps_m: float = 1.5,
+ point_density: float | None = None,
+ ransac_iterations: int = 200,
+ seed: int = 0,
+) -> list[RoofPlane]:
+ """Segment a roof point cloud into planes.
+
+ `points` is an (N, 3) array of SWEREF 99 TM (easting, northing, height).
+ Returns planes ordered by descending area. Determinism is deliberate: the
+ same cloud must always yield the same arrays, or an operator re-running a
+ derive would see the geometry shuffle for no reason.
+ """
+ from sklearn.cluster import DBSCAN # imported here to keep import cost off the CLI path
+
+ pts = np.asarray(points, dtype=float)
+ if pts.ndim != 2 or pts.shape[1] != 3:
+ raise ValueError(f"points must be (N, 3), got {pts.shape}")
+
+ rng = np.random.default_rng(seed)
+ remaining = pts
+ planes: list[RoofPlane] = []
+
+ for _ in range(max_planes):
+ if len(remaining) < min_plane_points:
+ break
+ mask = _ransac_plane(remaining, threshold_m, ransac_iterations, rng)
+ if mask is None or int(mask.sum()) < min_plane_points:
+ break
+ inliers = remaining[mask]
+ remaining = remaining[~mask]
+
+ # One plane equation can describe several disjoint faces. Split them.
+ labels = DBSCAN(eps=cluster_eps_m, min_samples=10).fit(inliers[:, :2]).labels_
+ for label in sorted(set(labels)):
+ if label == -1:
+ continue # DBSCAN noise
+ cluster = inliers[labels == label]
+ if len(cluster) < min_plane_points:
+ continue
+ normal = _fit_plane(cluster)
+ tilt, azimuth = _tilt_azimuth(normal)
+ if tilt > MAX_TILT_DEG:
+ continue # a wall, not a roof
+ planes.append(
+ RoofPlane(
+ tilt_deg=round(tilt, 1),
+ # Round before normalising: 359.97 rounds to 360.0, which is
+ # the same direction as 0 but reads as an out-of-range value.
+ azimuth_deg=round(azimuth, 1) % 360.0,
+ area_m2=round(_surface_area(cluster, tilt, point_density), 1),
+ point_count=len(cluster),
+ mean_height_m=round(float(cluster[:, 2].mean()), 2),
+ )
+ )
+
+ planes.sort(key=lambda p: p.area_m2, reverse=True)
+ return planes
diff --git a/roofmodel/ftw_roofmodel/sweref.py b/roofmodel/ftw_roofmodel/sweref.py
new file mode 100644
index 000000000..2d5e8f296
--- /dev/null
+++ b/roofmodel/ftw_roofmodel/sweref.py
@@ -0,0 +1,191 @@
+"""SWEREF 99 TM <-> WGS84 conversion.
+
+Lantmaeteriet publishes everything in SWEREF 99 TM (EPSG:3006) while FTW stores
+site location as WGS84 latitude/longitude, so every bounding box we send and
+every point cloud we read has to cross this boundary.
+
+This is Lantmaeteriet's own published Gauss conformal projection algorithm
+(Krueger series), implemented directly rather than pulled in via pyproj. The
+reason is proportion: pyproj ships a full PROJ build for what is, here, exactly
+one projection with fixed parameters. The series below is accurate to well under
+a millimetre across Sweden, which is several orders of magnitude finer than the
+1-2 points/m2 LiDAR it is used to place.
+
+SWEREF 99 TM is a transverse Mercator on GRS 80 with central meridian 15 deg E,
+scale factor 0.9996, false easting 500 000 m and false northing 0.
+"""
+
+from __future__ import annotations
+
+import math
+
+# GRS 80 ellipsoid.
+_A = 6378137.0
+_F = 1.0 / 298.257222101
+
+# SWEREF 99 TM projection parameters.
+_CENTRAL_MERIDIAN = 15.0
+_SCALE = 0.9996
+_FALSE_EASTING = 500000.0
+_FALSE_NORTHING = 0.0
+
+# Derived constants, computed once.
+_E2 = _F * (2.0 - _F)
+_N = _F / (2.0 - _F)
+_A_HAT = _A / (1.0 + _N) * (1.0 + _N**2 / 4.0 + _N**4 / 64.0)
+
+
+def _forward_coefficients() -> tuple[float, float, float, float]:
+ n = _N
+ return (
+ n / 2.0 - 2.0 * n**2 / 3.0 + 5.0 * n**3 / 16.0 + 41.0 * n**4 / 180.0,
+ 13.0 * n**2 / 48.0 - 3.0 * n**3 / 5.0 + 557.0 * n**4 / 1440.0,
+ 61.0 * n**3 / 240.0 - 103.0 * n**4 / 140.0,
+ 49561.0 * n**4 / 161280.0,
+ )
+
+
+def _inverse_coefficients() -> tuple[float, float, float, float]:
+ n = _N
+ return (
+ n / 2.0 - 2.0 * n**2 / 3.0 + 37.0 * n**3 / 96.0 - n**4 / 360.0,
+ n**2 / 48.0 + n**3 / 15.0 - 437.0 * n**4 / 1440.0,
+ 17.0 * n**3 / 480.0 - 37.0 * n**4 / 840.0,
+ 4397.0 * n**4 / 161280.0,
+ )
+
+
+def wgs84_to_sweref99tm(lat: float, lon: float) -> tuple[float, float]:
+ """Convert WGS84 degrees to SWEREF 99 TM (northing, easting) in metres."""
+ phi = math.radians(lat)
+ lam = math.radians(lon)
+ lam0 = math.radians(_CENTRAL_MERIDIAN)
+
+ e2 = _E2
+ a_coef = e2
+ b_coef = (5.0 * e2**2 - e2**3) / 6.0
+ c_coef = (104.0 * e2**3 - 45.0 * e2**4) / 120.0
+ d_coef = 1237.0 * e2**4 / 1260.0
+
+ sin_phi = math.sin(phi)
+ phi_star = phi - sin_phi * math.cos(phi) * (
+ a_coef
+ + b_coef * sin_phi**2
+ + c_coef * sin_phi**4
+ + d_coef * sin_phi**6
+ )
+
+ dlam = lam - lam0
+ xi_p = math.atan2(math.tan(phi_star), math.cos(dlam))
+ eta_p = math.atanh(math.cos(phi_star) * math.sin(dlam))
+
+ b1, b2, b3, b4 = _forward_coefficients()
+ scaled = _SCALE * _A_HAT
+ northing = _FALSE_NORTHING + scaled * (
+ xi_p
+ + b1 * math.sin(2 * xi_p) * math.cosh(2 * eta_p)
+ + b2 * math.sin(4 * xi_p) * math.cosh(4 * eta_p)
+ + b3 * math.sin(6 * xi_p) * math.cosh(6 * eta_p)
+ + b4 * math.sin(8 * xi_p) * math.cosh(8 * eta_p)
+ )
+ easting = _FALSE_EASTING + scaled * (
+ eta_p
+ + b1 * math.cos(2 * xi_p) * math.sinh(2 * eta_p)
+ + b2 * math.cos(4 * xi_p) * math.sinh(4 * eta_p)
+ + b3 * math.cos(6 * xi_p) * math.sinh(6 * eta_p)
+ + b4 * math.cos(8 * xi_p) * math.sinh(8 * eta_p)
+ )
+ return northing, easting
+
+
+def sweref99tm_to_wgs84(northing: float, easting: float) -> tuple[float, float]:
+ """Convert SWEREF 99 TM (northing, easting) in metres to WGS84 degrees."""
+ scaled = _SCALE * _A_HAT
+ xi = (northing - _FALSE_NORTHING) / scaled
+ eta = (easting - _FALSE_EASTING) / scaled
+
+ d1, d2, d3, d4 = _inverse_coefficients()
+ xi_p = (
+ xi
+ - d1 * math.sin(2 * xi) * math.cosh(2 * eta)
+ - d2 * math.sin(4 * xi) * math.cosh(4 * eta)
+ - d3 * math.sin(6 * xi) * math.cosh(6 * eta)
+ - d4 * math.sin(8 * xi) * math.cosh(8 * eta)
+ )
+ eta_p = (
+ eta
+ - d1 * math.cos(2 * xi) * math.sinh(2 * eta)
+ - d2 * math.cos(4 * xi) * math.sinh(4 * eta)
+ - d3 * math.cos(6 * xi) * math.sinh(6 * eta)
+ - d4 * math.cos(8 * xi) * math.sinh(8 * eta)
+ )
+
+ phi_star = math.asin(math.sin(xi_p) / math.cosh(eta_p))
+ dlam = math.atan2(math.sinh(eta_p), math.cos(xi_p))
+
+ e2 = _E2
+ a_star = e2 + e2**2 + e2**3 + e2**4
+ b_star = -(7.0 * e2**2 + 17.0 * e2**3 + 30.0 * e2**4) / 6.0
+ c_star = (224.0 * e2**3 + 889.0 * e2**4) / 120.0
+ d_star = -(4279.0 * e2**4) / 1260.0
+
+ sin_ps = math.sin(phi_star)
+ phi = phi_star + sin_ps * math.cos(phi_star) * (
+ a_star
+ + b_star * sin_ps**2
+ + c_star * sin_ps**4
+ + d_star * sin_ps**6
+ )
+ lam = math.radians(_CENTRAL_MERIDIAN) + dlam
+ return math.degrees(phi), math.degrees(lam)
+
+
+def bbox_wgs84_to_sweref99tm(
+ min_lat: float, min_lon: float, max_lat: float, max_lon: float
+) -> tuple[float, float, float, float]:
+ """Project a WGS84 bounding box to a SWEREF 99 TM (minE, minN, maxE, maxN).
+
+ All four corners are projected and the extremes taken, rather than just the
+ two diagonal corners: the projection is not axis-aligned, so a box's
+ projected edges bow outward and the diagonal-only result would clip.
+ """
+ corners = [
+ wgs84_to_sweref99tm(min_lat, min_lon),
+ wgs84_to_sweref99tm(min_lat, max_lon),
+ wgs84_to_sweref99tm(max_lat, min_lon),
+ wgs84_to_sweref99tm(max_lat, max_lon),
+ ]
+ northings = [c[0] for c in corners]
+ eastings = [c[1] for c in corners]
+ return min(eastings), min(northings), max(eastings), max(northings)
+
+
+def stac_search_bbox(
+ lat: float, lon: float, radius_m: float, bbox_epsg: int = 3006
+) -> tuple[float, float, float, float]:
+ """Bounding box around a site, in the CRS a STAC catalog expects.
+
+ EPSG:3006 (the default) is what Lantmaeteriet's catalog takes; EPSG:4326
+ in lon/lat order is what the STAC spec itself mandates, for catalogs that
+ follow it. Anything else would need a projection stack this module
+ deliberately does not carry.
+ """
+ south, west, north, east = metre_box_around(lat, lon, radius_m)
+ if bbox_epsg == 4326:
+ return (west, south, east, north)
+ if bbox_epsg == 3006:
+ return bbox_wgs84_to_sweref99tm(south, west, north, east)
+ raise ValueError(f"unsupported bbox EPSG {bbox_epsg}; use 3006 or 4326")
+
+
+def metre_box_around(lat: float, lon: float, radius_m: float) -> tuple[float, float, float, float]:
+ """Return a WGS84 (min_lat, min_lon, max_lat, max_lon) box of +/- radius_m.
+
+ Built by projecting the centre, stepping in metres in SWEREF 99 TM, and
+ unprojecting: doing it that way keeps the box square on the ground instead
+ of stretching with latitude the way a naive degree offset would.
+ """
+ n, e = wgs84_to_sweref99tm(lat, lon)
+ south, west = sweref99tm_to_wgs84(n - radius_m, e - radius_m)
+ north, east = sweref99tm_to_wgs84(n + radius_m, e + radius_m)
+ return south, west, north, east
diff --git a/roofmodel/pyproject.toml b/roofmodel/pyproject.toml
new file mode 100644
index 000000000..4e138a246
--- /dev/null
+++ b/roofmodel/pyproject.toml
@@ -0,0 +1,27 @@
+[build-system]
+requires = ["setuptools>=68"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "ftw-roofmodel"
+version = "0.1.0"
+description = "Derive PV array geometry from Lantmateriet building and LiDAR open data"
+requires-python = ">=3.10"
+# Core dependencies are permissive-licensed and light enough to install on a Pi.
+dependencies = [
+ "numpy>=1.24",
+ "scikit-learn>=1.3",
+ "requests>=2.31",
+]
+
+[project.optional-dependencies]
+# LAZ decoding pulls a compiled backend, so it is opt-in. Everything except the
+# point-cloud read works without it, which keeps the module testable in CI.
+geo = ["laspy[lazrs]>=2.5"]
+test = ["pytest==8.4.2"]
+
+[tool.setuptools.packages.find]
+include = ["ftw_roofmodel*"]
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]
diff --git a/roofmodel/tests/__init__.py b/roofmodel/tests/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/roofmodel/tests/test_assets.py b/roofmodel/tests/test_assets.py
new file mode 100644
index 000000000..107a5baef
--- /dev/null
+++ b/roofmodel/tests/test_assets.py
@@ -0,0 +1,212 @@
+"""Choosing STAC assets by what they are, and reading them in ranges.
+
+Both Lantmaeteriet products are STAC APIs; they differ in what their items point
+at. Byggnad-vektor delivers GeoPackage, Laserdata Skog delivers LAZ organised as
+COPC. Selecting on media type rather than on an asset key is what keeps that
+difference from becoming a pile of special cases.
+"""
+
+from __future__ import annotations
+
+import io
+
+import pytest
+
+from ftw_roofmodel.geotorget import (
+ MEDIA_COPC,
+ MEDIA_GEOJSON,
+ MEDIA_GEOPACKAGE,
+ MEDIA_LAZ,
+ Asset,
+ StacItem,
+ _item_from_feature,
+ media_type_for,
+)
+from ftw_roofmodel.pointcloud import HttpRangeFile, PointCloudError, bounds_of
+
+
+def item(assets):
+ return StacItem("i", "c", assets, None, raw={})
+
+
+def test_a_bare_href_still_gets_a_type_from_its_extension():
+ """Tests and simple catalogues pass strings; selection must still work."""
+ it = item({"data": "https://x/tile.copc.laz"})
+ assert isinstance(it.assets["data"], Asset)
+ assert it.pick(MEDIA_COPC).href == "https://x/tile.copc.laz"
+
+
+def test_copc_is_recognised_before_plain_laz():
+ """A COPC file is also a .laz; reading it as one costs the whole tile."""
+ assert media_type_for("https://x/y/tile.copc.laz") == MEDIA_COPC
+ assert media_type_for("https://x/y/tile.laz") == MEDIA_LAZ
+
+
+def test_query_strings_do_not_hide_the_extension():
+ """Signed download URLs carry tokens after a '?'."""
+ assert media_type_for("https://x/tile.gpkg?token=abc&x=1") == MEDIA_GEOPACKAGE
+
+
+def test_a_declared_type_beats_the_extension():
+ """The catalogue knows better than the filename."""
+ a = Asset(href="https://x/download", media_type=MEDIA_GEOPACKAGE)
+ assert a.effective_media_type == MEDIA_GEOPACKAGE
+
+
+def test_preference_order_is_honoured():
+ it = item({
+ "laz": Asset("https://x/t.laz", MEDIA_LAZ),
+ "copc": Asset("https://x/t.copc.laz", MEDIA_COPC),
+ })
+ assert it.pick(MEDIA_COPC, MEDIA_LAZ).effective_media_type == MEDIA_COPC
+ assert it.pick(MEDIA_LAZ, MEDIA_COPC).effective_media_type == MEDIA_LAZ
+
+
+def test_falls_back_to_the_data_role_when_the_type_is_unknown():
+ it = item({
+ "thumbnail": Asset("https://x/preview.png", "image/png", roles=("thumbnail",)),
+ "mystery": Asset("https://x/blob", None, roles=("data",)),
+ })
+ assert it.pick(MEDIA_COPC).href == "https://x/blob"
+
+
+def test_a_lone_asset_is_unambiguous_whatever_it_is_called():
+ assert item({"whatever": Asset("https://x/blob")}).pick(MEDIA_COPC).href == "https://x/blob"
+
+
+def test_several_unlabelled_assets_are_refused_rather_than_guessed():
+ it = item({"a": Asset("https://x/a"), "b": Asset("https://x/b")})
+ assert it.pick(MEDIA_COPC) is None
+
+
+def test_stac_assets_keep_their_type_and_roles():
+ feature = {
+ "id": "tile-1",
+ "collection": "laserdata-nedladdning-skog",
+ "assets": {
+ "data": {
+ "href": "https://x/t.copc.laz",
+ "type": MEDIA_COPC,
+ "roles": ["data"],
+ "title": "Punktmoln",
+ }
+ },
+ "properties": {},
+ }
+ parsed = _item_from_feature(feature)
+ asset = parsed.pick(MEDIA_COPC)
+ assert asset.media_type == MEDIA_COPC
+ assert asset.roles == ("data",)
+ assert asset.title == "Punktmoln"
+
+
+def test_geojson_and_geopackage_are_both_selectable():
+ it = item({"gj": Asset("https://x/b.geojson"), "gp": Asset("https://x/b.gpkg")})
+ assert it.pick(MEDIA_GEOPACKAGE, MEDIA_GEOJSON).href == "https://x/b.gpkg"
+ assert it.pick(MEDIA_GEOJSON, MEDIA_GEOPACKAGE).href == "https://x/b.geojson"
+
+
+# --- range reads ------------------------------------------------------------
+
+
+class FakeResponse:
+ def __init__(self, status, content=b"", headers=None):
+ self.status_code = status
+ self.content = content
+ self.headers = headers or {}
+
+
+class RangeServer:
+ """Serves a byte string over Range, and counts what was actually moved."""
+
+ def __init__(self, body: bytes, *, supports_range: bool = True):
+ self.body = body
+ self.supports_range = supports_range
+ self.requests: list[str] = []
+
+ def head(self, url, timeout=None):
+ return FakeResponse(200, headers={"Content-Length": str(len(self.body))})
+
+ def get(self, url, headers=None, timeout=None):
+ rng = (headers or {}).get("Range")
+ if not self.supports_range or not rng:
+ self.requests.append("full")
+ return FakeResponse(200, self.body)
+ self.requests.append(rng)
+ spec = rng.split("=", 1)[1]
+ start, end = spec.split("-")
+ lo = int(start)
+ hi = min(int(end), len(self.body) - 1)
+ return FakeResponse(206, self.body[lo : hi + 1])
+
+
+BODY = bytes(range(256)) * 40 # 10 240 bytes, every offset distinguishable
+
+
+def test_reads_a_window_without_moving_the_whole_file():
+ server = RangeServer(BODY)
+ fh = HttpRangeFile(server, "https://x/t.copc.laz", chunk_bytes=512)
+ fh.seek(1000)
+ assert fh.read(16) == BODY[1000:1016]
+ assert fh.bytes_fetched == 512, "one chunk, not the whole file"
+ assert len(server.requests) == 1
+
+
+def test_seek_and_tell_track_the_position():
+ fh = HttpRangeFile(RangeServer(BODY), "https://x/t", chunk_bytes=64)
+ assert fh.seek(100) == 100 and fh.tell() == 100
+ fh.read(10)
+ assert fh.tell() == 110
+ assert fh.seek(-10, io.SEEK_END) == len(BODY) - 10
+ assert fh.seek(5, io.SEEK_CUR) == len(BODY) - 5
+
+
+def test_a_second_read_inside_the_chunk_costs_no_request():
+ server = RangeServer(BODY)
+ fh = HttpRangeFile(server, "https://x/t", chunk_bytes=1024)
+ fh.seek(0)
+ fh.read(8)
+ before = len(server.requests)
+ fh.seek(64)
+ assert fh.read(8) == BODY[64:72]
+ assert len(server.requests) == before, "the chunk was already held"
+
+
+def test_reading_across_the_chunk_boundary_fetches_again():
+ server = RangeServer(BODY)
+ fh = HttpRangeFile(server, "https://x/t", chunk_bytes=128)
+ fh.seek(0)
+ assert fh.read(8) == BODY[0:8]
+ fh.seek(4096)
+ assert fh.read(8) == BODY[4096:4104]
+ assert len(server.requests) == 2
+
+
+def test_a_host_that_ignores_range_is_detected_not_trusted():
+ """A 200 means the whole body arrived; treating it as partial corrupts."""
+ fh = HttpRangeFile(RangeServer(BODY, supports_range=False), "https://x/t")
+ fh.seek(10)
+ with pytest.raises(PointCloudError, match="range requests"):
+ fh.read(4)
+
+
+def test_a_host_that_will_not_report_a_size_is_refused():
+ class NoLength:
+ def head(self, url, timeout=None):
+ return FakeResponse(200, headers={})
+
+ with pytest.raises(PointCloudError, match="size"):
+ HttpRangeFile(NoLength(), "https://x/t").size
+
+
+def test_readinto_fills_the_buffer():
+ fh = HttpRangeFile(RangeServer(BODY), "https://x/t", chunk_bytes=256)
+ fh.seek(32)
+ buf = bytearray(16)
+ assert fh.readinto(buf) == 16
+ assert bytes(buf) == BODY[32:48]
+
+
+def test_bounds_pad_the_footprint_so_the_eaves_survive():
+ ring = [(100.0, 200.0), (110.0, 200.0), (110.0, 220.0), (100.0, 220.0)]
+ assert bounds_of(ring, 1.0) == (99.0, 199.0, 111.0, 221.0)
diff --git a/roofmodel/tests/test_buildings.py b/roofmodel/tests/test_buildings.py
new file mode 100644
index 000000000..90805abc7
--- /dev/null
+++ b/roofmodel/tests/test_buildings.py
@@ -0,0 +1,236 @@
+"""Building lookup, frame detection and footprint clipping."""
+
+from __future__ import annotations
+
+import math
+
+import numpy as np
+import pytest
+
+from ftw_roofmodel import sweref
+from ftw_roofmodel.buildings import (
+ Building,
+ BuildingLookupError,
+ buildings_from_features,
+ clip_to_footprint,
+ inflate_ring,
+ point_in_ring,
+ search_buildings,
+)
+from ftw_roofmodel.geotorget import COLLECTION_BUILDINGS, Credentials, GeotorgetClient
+
+STOCKHOLM = (59.33, 18.07)
+
+
+def square_ring(cx, cy, side):
+ h = side / 2.0
+ return [(cx - h, cy - h), (cx + h, cy - h), (cx + h, cy + h), (cx - h, cy + h)]
+
+
+class FakeResponse:
+ def __init__(self, payload, status_code=200):
+ self._payload = payload
+ self.status_code = status_code
+
+ def json(self):
+ return self._payload
+
+
+class FakeSession:
+ """Records what was asked for and replays a canned STAC response."""
+
+ def __init__(self, payload):
+ self.payload = payload
+ self.posts = []
+
+ def post(self, url, json=None, timeout=None):
+ self.posts.append((url, json))
+ return FakeResponse(self.payload)
+
+
+def stac_feature(ring, feature_id="bldg-1", **props):
+ return {
+ "id": feature_id,
+ "collection": COLLECTION_BUILDINGS,
+ "geometry": {"type": "Polygon", "coordinates": [[list(p) for p in ring] + [list(ring[0])]]},
+ "properties": props,
+ "assets": {},
+ }
+
+
+def test_area_and_centroid_of_a_known_square():
+ ring = square_ring(674000.0, 6580000.0, 10.0)
+ [b] = buildings_from_features(
+ [{"geometry": {"type": "Polygon", "coordinates": [ring]}, "id": "sq"}],
+ latitude=STOCKHOLM[0], longitude=STOCKHOLM[1],
+ )
+ assert b.area_m2 == pytest.approx(100.0)
+ cx, cy = b.centroid_sweref()
+ assert (cx, cy) == pytest.approx((674000.0, 6580000.0))
+
+
+def test_wgs84_rings_are_projected_before_measuring():
+ """A ring in degrees must be recognised and converted, not measured raw."""
+ lat, lon = STOCKHOLM
+ n, e = sweref.wgs84_to_sweref99tm(lat, lon)
+ ring_sweref = square_ring(e, n, 12.0)
+ ring_wgs84 = []
+ for x, y in ring_sweref:
+ blat, blon = sweref.sweref99tm_to_wgs84(y, x) # ring is (E, N)
+ ring_wgs84.append([blon, blat]) # GeoJSON is [lon, lat]
+
+ [b] = buildings_from_features(
+ [{"geometry": {"type": "Polygon", "coordinates": [ring_wgs84]}, "id": "deg"}],
+ latitude=lat, longitude=lon,
+ )
+ # 12 m square, recovered through a full round trip through degrees.
+ assert b.area_m2 == pytest.approx(144.0, abs=0.5)
+
+
+def test_tiles_and_slivers_are_not_offered_as_buildings():
+ lat, lon = STOCKHOLM
+ n, e = sweref.wgs84_to_sweref99tm(lat, lon)
+ feats = [
+ {"geometry": {"type": "Polygon", "coordinates": [square_ring(e, n, 2500.0)]}, "id": "tile"},
+ {"geometry": {"type": "Polygon", "coordinates": [square_ring(e, n, 1.0)]}, "id": "sliver"},
+ {"geometry": {"type": "Polygon", "coordinates": [square_ring(e, n, 11.0)]}, "id": "house"},
+ ]
+ got = buildings_from_features(feats, latitude=lat, longitude=lon)
+ assert [b.building_id for b in got] == ["house"]
+
+
+def test_candidates_come_back_nearest_first():
+ lat, lon = STOCKHOLM
+ n, e = sweref.wgs84_to_sweref99tm(lat, lon)
+ feats = [
+ {"geometry": {"type": "Polygon", "coordinates": [square_ring(e + 60, n, 10.0)]}, "id": "far"},
+ {"geometry": {"type": "Polygon", "coordinates": [square_ring(e + 5, n, 10.0)]}, "id": "near"},
+ {"geometry": {"type": "Polygon", "coordinates": [square_ring(e + 25, n, 10.0)]}, "id": "mid"},
+ ]
+ got = buildings_from_features(feats, latitude=lat, longitude=lon)
+ assert [b.building_id for b in got] == ["near", "mid", "far"]
+ assert got[0].distance_m < got[1].distance_m < got[2].distance_m
+
+
+def test_multipolygon_yields_one_candidate_per_part():
+ lat, lon = STOCKHOLM
+ n, e = sweref.wgs84_to_sweref99tm(lat, lon)
+ feat = {
+ "id": "pair",
+ "geometry": {
+ "type": "MultiPolygon",
+ "coordinates": [[square_ring(e, n, 10.0)], [square_ring(e + 30, n, 12.0)]],
+ },
+ }
+ got = buildings_from_features([feat], latitude=lat, longitude=lon)
+ assert len(got) == 2
+ assert len({b.building_id for b in got}) == 2, "parts must not share an id"
+
+
+def test_search_queries_the_building_collection_and_maps_results():
+ lat, lon = STOCKHOLM
+ n, e = sweref.wgs84_to_sweref99tm(lat, lon)
+ session = FakeSession({"features": [stac_feature(square_ring(e, n, 10.0), "b1")]})
+ client = GeotorgetClient(Credentials("u", "t"), session=session)
+
+ got = search_buildings(client, latitude=lat, longitude=lon)
+
+ assert [b.building_id for b in got] == ["b1"]
+ (_, body), = session.posts
+ assert body["collections"] == [COLLECTION_BUILDINGS]
+ # The bbox is WGS84 lon/lat per the STAC spec — verified against the live
+ # Lantmaeteriet service, which follows it.
+ assert 17.9 < body["bbox"][0] < 18.2, body["bbox"]
+ assert 59.2 < body["bbox"][1] < 59.5, body["bbox"]
+
+
+def test_search_says_what_to_do_when_nothing_comes_back():
+ client = GeotorgetClient(Credentials("u", "t"), session=FakeSession({"features": []}))
+ with pytest.raises(BuildingLookupError) as exc:
+ search_buildings(client, latitude=STOCKHOLM[0], longitude=STOCKHOLM[1])
+ assert "Byggnad" in str(exc.value)
+
+
+def test_point_in_ring_handles_edges_and_outside():
+ ring = square_ring(0.0, 0.0, 10.0)
+ assert point_in_ring(0.0, 0.0, ring)
+ assert point_in_ring(4.9, 4.9, ring)
+ assert not point_in_ring(5.1, 0.0, ring)
+ assert not point_in_ring(0.0, 99.0, ring)
+
+
+def test_inflate_ring_grows_the_outline():
+ ring = square_ring(0.0, 0.0, 10.0)
+ bigger = inflate_ring(ring, 1.0)
+ # Corners sit at radius 7.07; pushing 1 m out puts them at 8.07.
+ assert math.hypot(*bigger[0]) == pytest.approx(math.hypot(*ring[0]) + 1.0)
+ assert all(point_in_ring(x, y, bigger) for x, y in ring)
+
+
+def test_clip_keeps_the_building_and_drops_the_neighbours():
+ rng = np.random.default_rng(3)
+ mine = np.column_stack([
+ rng.uniform(-4, 4, 400), rng.uniform(-4, 4, 400), rng.uniform(0, 4, 400)])
+ theirs = np.column_stack([
+ rng.uniform(46, 54, 400), rng.uniform(-4, 4, 400), rng.uniform(0, 4, 400)])
+ cloud = np.vstack([mine, theirs])
+
+ kept = clip_to_footprint(cloud, square_ring(0.0, 0.0, 10.0), buffer_m=0.0)
+
+ assert len(kept) == len(mine)
+ assert kept[:, 0].max() < 10.0
+
+
+def test_clip_keeps_the_eaves():
+ """Roof returns overhang the wall line; clipping exactly would shave them."""
+ ring = square_ring(0.0, 0.0, 10.0)
+ eaves = np.array([[5.4, 0.0, 3.0], [-5.4, 0.0, 3.0], [0.0, 5.4, 3.0]])
+
+ assert len(clip_to_footprint(eaves, ring, buffer_m=0.0)) == 0
+ assert len(clip_to_footprint(eaves, ring, buffer_m=1.0)) == 3
+
+
+def test_clip_of_an_empty_cloud_is_empty_not_an_error():
+ assert len(clip_to_footprint(np.empty((0, 3)), square_ring(0, 0, 10))) == 0
+
+
+def test_geojson_feature_is_wgs84_and_closed():
+ lat, lon = STOCKHOLM
+ n, e = sweref.wgs84_to_sweref99tm(lat, lon)
+ b = Building("b1", square_ring(e, n, 10.0), 100.0, 0.0)
+ feat = b.to_geojson()
+
+ ring = feat["geometry"]["coordinates"][0]
+ assert ring[0] == ring[-1], "GeoJSON rings must close"
+ for x, y in ring:
+ assert -180 <= x <= 180 and -90 <= y <= 90
+ assert feat["properties"]["latitude"] == pytest.approx(lat, abs=1e-3)
+ assert feat["properties"]["longitude"] == pytest.approx(lon, abs=1e-3)
+
+
+def test_both_ring_orders_come_back_at_the_real_site():
+ """Axis order must not depend on where the ring came from.
+
+ A GeoPackage stores x=easting first; wgs84_to_sweref99tm returns northing
+ first; EPSG's registry declares EPSG:3006 north-first and some exports
+ follow it. The first mismatch shipped: GeoPackage-sourced buildings — the
+ normal Lantmäteriet case — reported their centroids in the Indian Ocean
+ (lat ≈ 4°) with 8 000 km distances, while every test asserted only areas
+ and SWEREF centroids, both of which are blind to a consistent swap.
+ """
+ lat, lon = STOCKHOLM
+ n, e = sweref.wgs84_to_sweref99tm(lat, lon)
+ ring_en = square_ring(e, n, 10.0) # as a GeoPackage stores it
+ ring_ne = [(y, x) for x, y in ring_en] # as the EPSG registry says
+
+ for ring in (ring_en, ring_ne):
+ [b] = buildings_from_features(
+ [{"geometry": {"type": "Polygon", "coordinates": [ring]}, "id": "b"}],
+ latitude=lat, longitude=lon,
+ )
+ feat = b.to_geojson()
+ assert feat["properties"]["latitude"] == pytest.approx(lat, abs=1e-3)
+ assert feat["properties"]["longitude"] == pytest.approx(lon, abs=1e-3)
+ assert b.distance_m < 50.0, (
+ f"a building drawn around the site is {b.distance_m:.0f} m away"
+ )
diff --git a/roofmodel/tests/test_derive_footprint.py b/roofmodel/tests/test_derive_footprint.py
new file mode 100644
index 000000000..466b8fae5
--- /dev/null
+++ b/roofmodel/tests/test_derive_footprint.py
@@ -0,0 +1,232 @@
+"""Deriving against a picked building footprint."""
+
+from __future__ import annotations
+
+import math
+
+import numpy as np
+import pytest
+
+from ftw_roofmodel import pipeline, sweref
+from ftw_roofmodel.buildings import BuildingLookupError, clip_to_footprint
+from ftw_roofmodel.geotorget import COLLECTION_BUILDINGS, Credentials, StacItem
+from ftw_roofmodel.pipeline import RoofModelError, derive
+from ftw_roofmodel.segment import segment_roof
+
+STOCKHOLM = (59.33, 18.07)
+
+
+def roof_face(tilt, azimuth, w, d, origin, density=8, noise=0.04, seed=1):
+ """Sample a tilted rectangle. Written as the inverse of what segment.py
+ computes, so recovering the tilt is a real test rather than a tautology."""
+ rng = np.random.default_rng(seed)
+ n = int(w * d * density)
+ x = rng.uniform(0, w, n)
+ y = rng.uniform(0, d, n)
+ az = math.radians(azimuth)
+ s = math.tan(math.radians(tilt))
+ z = -(x * math.sin(az) + y * math.cos(az)) * s + rng.normal(0, noise, n)
+ return np.column_stack([x + origin[0], y + origin[1], z + origin[2]])
+
+
+def ring(cx, cy, w, d):
+ return [(cx, cy), (cx + w, cy), (cx + w, cy + d), (cx, cy + d)]
+
+
+class FakeClient:
+ """Stands in for Geotorget: canned buildings, canned LiDAR."""
+
+ def __init__(self, buildings_payload, points):
+ self._buildings = buildings_payload
+ self._points = points
+ self.searched = []
+ self.searched_boxes = []
+
+ def search(self, collection, bbox, limit=20):
+ self.searched.append(collection)
+ self.searched_boxes.append((collection, bbox))
+ if collection == COLLECTION_BUILDINGS:
+ return [StacItem(f["id"], collection, {}, None, raw=f) for f in self._buildings]
+ return [StacItem("lidar-1", collection, {"data": "http://x/tile.laz"}, None, raw={})]
+
+ def download(self, url):
+ return b"laz-bytes"
+
+
+@pytest.fixture
+def scene(monkeypatch):
+ """A house and a neighbour that share a ridge orientation.
+
+ The neighbour is the point: an azimuth-180 plane is z = f(y) with no x term,
+ so it extends across the whole tile and the two buildings compete for each
+ other's returns unless the cloud is clipped first.
+ """
+ n, e = sweref.wgs84_to_sweref99tm(*STOCKHOLM)
+ mine = np.vstack([
+ roof_face(35, 180, 12, 6, (e, n, 0), seed=2),
+ roof_face(35, 0, 12, 6, (e, n + 6, 4.2), seed=3),
+ ])
+ neighbour = np.vstack([
+ roof_face(35, 180, 12, 6, (e + 40, n, 0), seed=4),
+ roof_face(35, 0, 12, 6, (e + 40, n + 6, 4.2), seed=5),
+ ])
+ cloud = np.vstack([mine, neighbour])
+
+ buildings = [
+ {"id": "mine", "geometry": {"type": "Polygon",
+ "coordinates": [[list(p) for p in ring(e, n, 12, 12)]]}, "properties": {}},
+ {"id": "neighbour", "geometry": {"type": "Polygon",
+ "coordinates": [[list(p) for p in ring(e + 40, n, 12, 12)]]}, "properties": {}},
+ ]
+ monkeypatch.setattr(pipeline, "load_points", lambda data: cloud)
+ return FakeClient(buildings, cloud), cloud, (e, n)
+
+
+def test_derive_clips_to_the_picked_building(scene):
+ client, cloud, _ = scene
+ model = derive(
+ latitude=STOCKHOLM[0], longitude=STOCKHOLM[1],
+ credentials=Credentials("u", "t"), client=client, building_id="mine",
+ )
+
+ b = model["building"]
+ assert b["building_id"] == "mine"
+ assert b["returns_in_radius"] == len(cloud)
+ # Roughly half the tile is the neighbour's, and it must be gone.
+ assert b["returns_used"] < b["returns_in_radius"] * 0.6
+ assert model["arrays"], "a clipped house still has a south roof"
+
+
+def test_clipping_recovers_the_true_area_the_neighbour_would_have_stolen(scene):
+ """The measurable payoff: coplanar buildings stop eating each other."""
+ _, cloud, (e, n) = scene
+ truth = 12 * 6 / math.cos(math.radians(35))
+
+ whole_tile = [p for p in segment_roof(cloud, point_density=8.0)
+ if abs(p.azimuth_deg - 180) < 5]
+ clipped = [p for p in segment_roof(clip_to_footprint(cloud, ring(e, n, 12, 12)),
+ point_density=8.0)
+ if abs(p.azimuth_deg - 180) < 5]
+
+ assert len(clipped) == 1, "one building has one south face"
+ assert clipped[0].area_m2 == pytest.approx(truth, rel=0.10)
+ # Unclipped, the two south faces are coplanar and merge into one oversized
+ # segment, so the site's own roof cannot be measured at all.
+ assert len(whole_tile) != 1 or whole_tile[0].area_m2 > truth * 1.5
+
+
+def test_derive_without_a_building_id_does_not_search_for_buildings(scene):
+ client, _, _ = scene
+ model = derive(
+ latitude=STOCKHOLM[0], longitude=STOCKHOLM[1],
+ credentials=Credentials("u", "t"), client=client,
+ )
+ assert COLLECTION_BUILDINGS not in client.searched
+ assert model["building"] is None
+
+
+def test_derive_rejects_a_building_id_that_is_not_there(scene):
+ client, _, _ = scene
+ with pytest.raises(RoofModelError) as exc:
+ derive(
+ latitude=STOCKHOLM[0], longitude=STOCKHOLM[1],
+ credentials=Credentials("u", "t"), client=client, building_id="not-a-building",
+ )
+ assert "not found" in str(exc.value)
+
+
+def test_derive_explains_a_footprint_with_no_returns_on_it(monkeypatch, scene):
+ """A building newer than the scan is a real case and needs a real message."""
+ client, _, (e, n) = scene
+ empty = np.empty((0, 3))
+ monkeypatch.setattr(pipeline, "load_points", lambda data: np.vstack([
+ roof_face(35, 180, 12, 6, (e + 400, n, 0), seed=9)]))
+
+ with pytest.raises(RoofModelError) as exc:
+ derive(
+ latitude=STOCKHOLM[0], longitude=STOCKHOLM[1],
+ credentials=Credentials("u", "t"), client=client, building_id="mine",
+ )
+ msg = str(exc.value)
+ assert "fall on building" in msg and "newer than the scan" in msg
+
+
+def test_derive_clips_to_a_drawn_footprint(scene):
+ """Where the catalog has no building dataset, the operator traces the
+ outline by hand — and that must clip exactly like a picked building,
+ without any building search happening at all."""
+ client, cloud, (e, n) = scene
+ corners = []
+ for ee, nn in ring(e, n, 12, 12):
+ lat, lon = sweref.sweref99tm_to_wgs84(nn, ee)
+ corners.append([lon, lat])
+ corners.append(list(corners[0])) # GeoJSON rings close themselves
+
+ model = derive(
+ latitude=STOCKHOLM[0], longitude=STOCKHOLM[1],
+ credentials=Credentials("u", "t"), client=client, footprint=corners,
+ )
+
+ assert COLLECTION_BUILDINGS not in client.searched
+ b = model["building"]
+ assert b["building_id"] == "drawn-footprint"
+ assert b["returns_used"] < b["returns_in_radius"] * 0.6
+ assert model["arrays"], "a traced house still has a south roof"
+
+
+def test_a_drawn_footprint_needs_three_corners(scene):
+ client, _, _ = scene
+ with pytest.raises(Exception) as exc:
+ derive(
+ latitude=STOCKHOLM[0], longitude=STOCKHOLM[1],
+ credentials=Credentials("u", "t"), client=client,
+ footprint=[[18.06, 59.33], [18.07, 59.33]],
+ )
+ assert "three corners" in str(exc.value)
+
+
+def test_the_re_search_reaches_as_far_as_the_picker_did(monkeypatch, scene):
+ """A barn 100 m out is inside the picker's 150 m search but outside the
+ 40 m LiDAR radius. The derive's re-find must use the picker's reach, or
+ every such pick dies with "not found near this site"."""
+ client, _, _ = scene
+ real_search = pipeline.search_buildings
+ radii = []
+
+ def windowed(client_, *, latitude, longitude, radius_m, **kw):
+ radii.append(radius_m)
+ # Model the real windowing: the barn only comes back when the search
+ # reaches at least as far as the picker's default.
+ if radius_m < 150.0:
+ raise BuildingLookupError("nothing inside this window")
+ return real_search(client_, latitude=latitude, longitude=longitude,
+ radius_m=radius_m, **kw)
+
+ monkeypatch.setattr(pipeline, "search_buildings", windowed)
+ model = derive(
+ latitude=STOCKHOLM[0], longitude=STOCKHOLM[1],
+ credentials=Credentials("u", "t"), client=client, building_id="mine",
+ )
+ assert model["building"]["building_id"] == "mine"
+ assert radii and min(radii) >= 150.0
+
+
+def test_the_lidar_lookup_centres_on_the_picked_building(scene):
+ """The tile search must cover the roof being derived, not the pin: a
+ building at the search edge can sit on a different tile."""
+ client, _, (e, n) = scene
+ derive(
+ latitude=STOCKHOLM[0], longitude=STOCKHOLM[1],
+ credentials=Credentials("u", "t"), client=client, building_id="neighbour",
+ )
+ lidar_bboxes = [b for c, b in getattr(client, "searched_boxes", [])
+ if c != COLLECTION_BUILDINGS]
+ assert lidar_bboxes, "the LiDAR collection was searched"
+ west, south, east, north = lidar_bboxes[-1]
+ centre_n, centre_e = sweref.wgs84_to_sweref99tm(
+ (south + north) / 2.0, (west + east) / 2.0
+ )
+ # The neighbour's centroid is 40 m east of the site; the bbox centre must
+ # follow it rather than stay on the pin.
+ assert centre_e == pytest.approx(e + 40 + 6, abs=15.0)
+ assert centre_n == pytest.approx(n + 6, abs=15.0)
diff --git a/roofmodel/tests/test_geopackage.py b/roofmodel/tests/test_geopackage.py
new file mode 100644
index 000000000..54c49e424
--- /dev/null
+++ b/roofmodel/tests/test_geopackage.py
@@ -0,0 +1,201 @@
+"""Decoding GeoPackage, the format Lantmaeteriet ships building vectors in.
+
+The fixtures are built byte by byte from the published layouts rather than by
+round-tripping the decoder, so a decoder that agrees with itself but not with
+the standard still fails here.
+"""
+
+from __future__ import annotations
+
+import os
+import sqlite3
+import struct
+import tempfile
+
+import pytest
+
+from ftw_roofmodel.geopackage import (
+ GeoPackageError,
+ parse_geometry_blob,
+ read_features,
+)
+
+SWEREF = 3006
+
+
+def wkb_polygon(rings, *, little=True, z=False):
+ """Standard WKB polygon, per OGC 06-103r4 clause 8.2."""
+ e = "<" if little else ">"
+ code = 1003 if z else 3
+ out = struct.pack("B", 1 if little else 0) + struct.pack(e + "I", code)
+ out += struct.pack(e + "I", len(rings))
+ for ring in rings:
+ out += struct.pack(e + "I", len(ring))
+ for point in ring:
+ out += struct.pack(e + "dd", point[0], point[1])
+ if z:
+ out += struct.pack(e + "d", point[2] if len(point) > 2 else 0.0)
+ return out
+
+
+def wkb_multipolygon(polygons, *, little=True):
+ e = "<" if little else ">"
+ out = struct.pack("B", 1 if little else 0) + struct.pack(e + "I", 6)
+ out += struct.pack(e + "I", len(polygons))
+ for rings in polygons:
+ out += wkb_polygon(rings, little=little)
+ return out
+
+
+def gpkg_blob(wkb, *, envelope=None, srs_id=SWEREF, little=True, empty=False):
+ """GeoPackage geometry BLOB, per OGC 12-128r19 clause 2.1.3.
+
+ Header is magic(2) + version(1) + flags(1) + srs_id(4), then the envelope,
+ then the WKB.
+ """
+ indicator = 0 if envelope is None else 1
+ flags = (1 if little else 0) | (indicator << 1) | (0x10 if empty else 0)
+ e = "<" if little else ">"
+ header = b"GP" + bytes([0, flags]) + struct.pack(e + "i", srs_id)
+ assert len(header) == 8, "the GeoPackage header is 8 bytes before the envelope"
+ body = b""
+ if envelope is not None:
+ body = struct.pack(e + "dddd", *envelope)
+ assert len(body) == 32, "an xy envelope is four doubles"
+ return header + body + wkb
+
+
+SQUARE = [[(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)]]
+
+
+def test_reads_a_plain_polygon():
+ geom = parse_geometry_blob(gpkg_blob(wkb_polygon(SQUARE)))
+ assert geom["type"] == "Polygon"
+ assert geom["coordinates"][0][0] == [0.0, 0.0]
+ assert geom["coordinates"][0][2] == [10.0, 10.0]
+ assert len(geom["coordinates"][0]) == 5
+
+
+def test_skips_the_envelope_when_one_is_present():
+ """The envelope sits between the header and the WKB and must be stepped over."""
+ with_env = parse_geometry_blob(
+ gpkg_blob(wkb_polygon(SQUARE), envelope=(0.0, 10.0, 0.0, 10.0))
+ )
+ without = parse_geometry_blob(gpkg_blob(wkb_polygon(SQUARE)))
+ assert with_env == without
+
+
+def test_reads_big_endian_geometry():
+ assert parse_geometry_blob(
+ gpkg_blob(wkb_polygon(SQUARE, little=False), little=False)
+ ) == parse_geometry_blob(gpkg_blob(wkb_polygon(SQUARE)))
+
+
+def test_reads_3d_polygons_by_stepping_the_z():
+ """Building footprints carry heights; the z must not shift the ring."""
+ ring = [[(0.0, 0.0, 12.5), (10.0, 0.0, 12.5), (10.0, 10.0, 12.5), (0.0, 0.0, 12.5)]]
+ geom = parse_geometry_blob(gpkg_blob(wkb_polygon(ring, z=True)))
+ assert geom["coordinates"][0] == [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 0.0]]
+
+
+def test_reads_a_multipolygon_as_several_rings():
+ other = [[(20.0, 20.0), (30.0, 20.0), (30.0, 30.0), (20.0, 20.0)]]
+ geom = parse_geometry_blob(gpkg_blob(wkb_multipolygon([SQUARE, other])))
+ assert geom["type"] == "MultiPolygon"
+ assert len(geom["coordinates"]) == 2
+
+
+def test_an_interior_ring_survives():
+ """A courtyard is a second ring, and dropping it would inflate the roof."""
+ hole = [(2.0, 2.0), (4.0, 2.0), (4.0, 4.0), (2.0, 2.0)]
+ geom = parse_geometry_blob(gpkg_blob(wkb_polygon(SQUARE + [hole])))
+ assert len(geom["coordinates"]) == 2
+
+
+def test_empty_geometry_is_none_not_an_error():
+ assert parse_geometry_blob(gpkg_blob(b"", empty=True)) is None
+
+
+def test_rejects_a_blob_that_is_not_a_geopackage_geometry():
+ with pytest.raises(GeoPackageError, match="magic"):
+ parse_geometry_blob(b"XX" + bytes(20))
+
+
+def test_refuses_extended_geometry_rather_than_guessing():
+ blob = bytearray(gpkg_blob(wkb_polygon(SQUARE)))
+ blob[3] |= 0x20
+ with pytest.raises(GeoPackageError, match="extended"):
+ parse_geometry_blob(bytes(blob))
+
+
+def test_refuses_a_non_polygon_rather_than_mis_clipping():
+ point = struct.pack("B", 1) + struct.pack(" arrays ------------------------------------------------------
+
+def test_north_facing_pitched_roofs_are_not_proposed():
+ """At Swedish latitudes a north pitch yields too little to be worth adding
+ to an operator's config."""
+ planes = [
+ RoofPlane(tilt_deg=35, azimuth_deg=0, area_m2=60, point_count=200, mean_height_m=6),
+ RoofPlane(tilt_deg=35, azimuth_deg=180, area_m2=60, point_count=200, mean_height_m=6),
+ ]
+ arrays = planes_to_arrays(planes)
+ assert len(arrays) == 1
+ assert arrays[0].azimuth_deg == 180
+
+
+def test_flat_roofs_are_kept_regardless_of_building_orientation():
+ planes = [RoofPlane(tilt_deg=1.0, azimuth_deg=180, area_m2=90, point_count=300, mean_height_m=9)]
+ assert len(planes_to_arrays(planes)) == 1
+
+
+def test_tiny_faces_are_dropped():
+ """Dormers and porches are real surfaces but not candidate arrays."""
+ planes = [RoofPlane(tilt_deg=35, azimuth_deg=180, area_m2=3.0, point_count=50, mean_height_m=5)]
+ assert planes_to_arrays(planes) == []
+
+
+def test_arrays_get_readable_and_unique_names():
+ planes = [
+ RoofPlane(tilt_deg=35, azimuth_deg=180, area_m2=60, point_count=200, mean_height_m=6),
+ RoofPlane(tilt_deg=35, azimuth_deg=182, area_m2=40, point_count=150, mean_height_m=6),
+ RoofPlane(tilt_deg=35, azimuth_deg=270, area_m2=30, point_count=120, mean_height_m=6),
+ ]
+ arrays = planes_to_arrays(planes)
+ names = [a.name for a in arrays]
+ assert len(set(names)) == len(names), names
+ assert names[0] == "Roof south"
+ assert "Roof west" in names
+
+
+def test_array_json_matches_the_config_field_names():
+ """The document pre-fills weather.pv_arrays, so the keys must line up."""
+ planes = [RoofPlane(tilt_deg=35, azimuth_deg=180, area_m2=60, point_count=200, mean_height_m=6)]
+ payload = planes_to_arrays(planes)[0].to_json()
+ assert set(payload) >= {"name", "rated_w", "tilt_deg", "azimuth_deg"}
+ assert payload["rated_w"] > 0
+
+
+# --- end to end ------------------------------------------------------------
+
+def _patched_points(monkeypatch, cloud):
+ monkeypatch.setattr(pipeline, "load_points", lambda data: cloud)
+
+
+def test_derive_produces_a_versioned_document(monkeypatch):
+ gable = np.vstack([
+ make_plane(tilt_deg=35, azimuth_deg=180, width=12, depth=6, seed=11),
+ make_plane(tilt_deg=35, azimuth_deg=0, width=12, depth=6, origin=(0, 6, -4.2), seed=12),
+ ])
+ _patched_points(monkeypatch, gable)
+ session = FakeSession(search={"features": [feature()]}, asset=b"laz-bytes")
+ client = GeotorgetClient(CREDS, session=session)
+
+ model = pipeline.derive(
+ latitude=59.33, longitude=18.07, credentials=CREDS, client=client,
+ now=dt.datetime(2026, 7, 31, tzinfo=dt.timezone.utc),
+ )
+
+ assert model["schema_version"] == pipeline.SCHEMA_VERSION
+ assert model["source"]["provider"] == "lantmateriet"
+ assert model["captured_at_ms"] > 0
+ assert model["derived_at_ms"] == 1785456000000
+ # The gable's north face is dropped, so exactly the south face survives.
+ assert len(model["arrays"]) == 1
+ south = model["arrays"][0]
+ assert south["azimuth_deg"] == pytest.approx(180.0, abs=3.0)
+ assert south["tilt_deg"] == pytest.approx(35.0, abs=3.0)
+ assert south["rated_w"] > 0
+ # Must be JSON-serialisable for the subprocess contract.
+ json.dumps(model)
+
+
+def test_derive_searches_the_spec_bbox(monkeypatch):
+ _patched_points(monkeypatch, make_plane(tilt_deg=35, azimuth_deg=180))
+ session = FakeSession(search={"features": [feature()]}, asset=b"x")
+ pipeline.derive(
+ latitude=59.33, longitude=18.07, credentials=CREDS,
+ client=GeotorgetClient(CREDS, session=session), radius_m=40.0,
+ )
+ (_, body), = session.posts
+ # WGS84 lon/lat per the STAC spec, the live service's frame — and still
+ # square on the ground: built by stepping metres in SWEREF and
+ # unprojecting, which is what stac_search_bbox does.
+ expected = sweref.stac_search_bbox(59.33, 18.07, 40.0, 4326)
+ assert list(body["bbox"]) == pytest.approx(list(expected))
+
+
+def test_derive_uses_both_lantmateriet_roots_by_default(monkeypatch):
+ """The live service splits its STAC per product family: buildings on the
+ vector root, point clouds on the elevation root. The default derive must
+ build a client for each."""
+ made = []
+
+ class RecordingClient:
+ def __init__(self, credentials, base_url=None, **kw):
+ made.append(base_url)
+
+ def search(self, collection, bbox, limit=20):
+ return []
+
+ monkeypatch.setattr(pipeline, "GeotorgetClient", RecordingClient)
+ with pytest.raises(RoofModelError):
+ pipeline.derive(latitude=59.33, longitude=18.07, credentials=CREDS)
+ assert made == [
+ pipeline.geotorget.DEFAULT_BASE_URL,
+ pipeline.geotorget.DEFAULT_LIDAR_BASE_URL,
+ ]
+
+
+def test_derive_uses_one_client_for_a_custom_catalog(monkeypatch):
+ """A custom single-root catalog serves both collections itself."""
+ made = []
+
+ class RecordingClient:
+ def __init__(self, credentials, base_url=None, **kw):
+ made.append(base_url)
+
+ def search(self, collection, bbox, limit=20):
+ return []
+
+ monkeypatch.setattr(pipeline, "GeotorgetClient", RecordingClient)
+ with pytest.raises(RoofModelError):
+ pipeline.derive(
+ latitude=59.33, longitude=18.07, credentials=CREDS,
+ base_url="https://stac.example.org",
+ )
+ assert made == ["https://stac.example.org"]
+
+
+def test_derive_outside_sweden_says_so(monkeypatch):
+ """No tiles come back for a site Lantmaeteriet does not cover."""
+ session = FakeSession(search={"features": []})
+ with pytest.raises(RoofModelError) as exc:
+ pipeline.derive(
+ latitude=-33.87, longitude=151.21, credentials=CREDS,
+ client=GeotorgetClient(CREDS, session=session),
+ )
+ assert "Sweden only" in str(exc.value)
+
+
+def test_derive_reports_unknown_capture_date_without_failing(monkeypatch):
+ _patched_points(monkeypatch, make_plane(tilt_deg=35, azimuth_deg=180))
+ session = FakeSession(search={"features": [feature(datetime_value=None)]}, asset=b"x")
+ model = pipeline.derive(
+ latitude=59.33, longitude=18.07, credentials=CREDS,
+ client=GeotorgetClient(CREDS, session=session),
+ )
+ assert model["captured_at_ms"] is None
+ assert model["source"]["dataset_datetime"] is None
+ assert model["arrays"], "a missing provenance date must not block the model"
diff --git a/roofmodel/tests/test_pointcloud.py b/roofmodel/tests/test_pointcloud.py
new file mode 100644
index 000000000..ea3f6a101
--- /dev/null
+++ b/roofmodel/tests/test_pointcloud.py
@@ -0,0 +1,29 @@
+"""The COPC window query must never prune on z.
+
+Lantmäteriet's COPC writer keys octree nodes with a z origin that does not
+match the file's own cube (seen live: a node keyed to a slab 1.7 km underground
+holding points at +18..+42 m). Pruning by those keys silently drops the dense
+deep levels, so the query's vertical range has to cover every slab the cube
+could describe.
+"""
+
+from ftw_roofmodel.pointcloud import copc_query_z_range
+
+
+def test_the_query_covers_the_whole_octree_cube():
+ for center, half in [(332.51, 5000.005), (0.0, 128.0), (-250.0, 4096.0)]:
+ z_lo, z_hi = copc_query_z_range(center, half)
+ assert z_lo < center - half
+ assert z_hi > center + half
+
+
+def test_the_query_stays_bounded_for_scaled_int_filters():
+ # laspy's post-filter casts scaled bounds to int32; the range must stay
+ # proportional to the cube, never an arbitrary huge sentinel.
+ z_lo, z_hi = copc_query_z_range(332.51, 5000.005)
+ assert z_hi - z_lo < 10 * 5000.005
+
+
+def test_a_degenerate_cube_still_yields_a_range():
+ z_lo, z_hi = copc_query_z_range(100.0, 0.0)
+ assert z_lo < 100.0 < z_hi
diff --git a/roofmodel/tests/test_segment.py b/roofmodel/tests/test_segment.py
new file mode 100644
index 000000000..e8ec35451
--- /dev/null
+++ b/roofmodel/tests/test_segment.py
@@ -0,0 +1,188 @@
+"""Roof segmentation tests.
+
+Every case is a synthetic roof whose tilt, azimuth and area are known exactly by
+construction, so the assertions check recovered geometry against ground truth
+rather than against a previous run's output.
+"""
+
+import math
+
+import numpy as np
+import pytest
+
+from ftw_roofmodel.segment import RoofPlane, segment_roof
+
+
+def make_plane(
+ *,
+ tilt_deg: float,
+ azimuth_deg: float,
+ width: float = 10.0,
+ depth: float = 8.0,
+ origin=(0.0, 0.0, 0.0),
+ density: float = 8.0,
+ noise_m: float = 0.0,
+ seed: int = 1,
+) -> np.ndarray:
+ """Sample a rectangular sloped surface with a known tilt and azimuth.
+
+ The surface is generated by taking a horizontal grid and raising it along
+ the downslope direction, which is the inverse of what segment_roof does, so
+ the test never shares an implementation with the code under test.
+ """
+ rng = np.random.default_rng(seed)
+ n = max(int(width * depth * density), 60)
+ x = rng.uniform(0, width, n)
+ y = rng.uniform(0, depth, n)
+
+ # Aspect points the way the surface faces; height falls off along it.
+ az = math.radians(azimuth_deg)
+ east_dir, north_dir = math.sin(az), math.cos(az)
+ slope = math.tan(math.radians(tilt_deg))
+ z = -(x * east_dir + y * north_dir) * slope
+
+ if noise_m:
+ z = z + rng.normal(0.0, noise_m, n)
+
+ return np.column_stack([x + origin[0], y + origin[1], z + origin[2]])
+
+
+@pytest.mark.parametrize(
+ "tilt,azimuth",
+ [
+ (35.0, 180.0), # classic south-facing pitched roof
+ (30.0, 90.0), # east
+ (30.0, 270.0), # west
+ (45.0, 0.0), # north
+ (20.0, 225.0), # south-west
+ (60.0, 135.0), # steep south-east
+ ],
+)
+def test_recovers_known_tilt_and_azimuth(tilt, azimuth):
+ cloud = make_plane(tilt_deg=tilt, azimuth_deg=azimuth)
+ planes = segment_roof(cloud, point_density=8.0)
+ assert planes, "expected at least one plane"
+ got = planes[0]
+ assert got.tilt_deg == pytest.approx(tilt, abs=1.5)
+ # Compare on the circle so 359.5 vs 0.5 is a 1-degree error, not 359.
+ delta = abs((got.azimuth_deg - azimuth + 180.0) % 360.0 - 180.0)
+ assert delta < 2.0, f"azimuth {got.azimuth_deg} vs {azimuth}"
+
+
+def test_flat_roof_reports_flat_and_does_not_invent_an_azimuth():
+ """A horizontal surface has no meaningful aspect: its normal is vertical,
+ so the horizontal component is pure noise and could point anywhere."""
+ cloud = make_plane(tilt_deg=0.0, azimuth_deg=180.0, noise_m=0.02)
+ planes = segment_roof(cloud, point_density=8.0)
+ assert planes
+ assert planes[0].tilt_deg < 5.0
+ assert planes[0].azimuth_deg == 180.0
+
+
+def test_gable_roof_splits_into_two_opposing_faces():
+ """The canonical case: one ridge, two faces 180 degrees apart."""
+ south = make_plane(tilt_deg=35, azimuth_deg=180, depth=6, origin=(0, 0, 0), seed=2)
+ north = make_plane(tilt_deg=35, azimuth_deg=0, depth=6, origin=(0, 6, -4.2), seed=3)
+ planes = segment_roof(np.vstack([south, north]), point_density=8.0)
+
+ assert len(planes) >= 2, f"expected two faces, got {len(planes)}"
+ azimuths = sorted(p.azimuth_deg for p in planes[:2])
+ opposed = abs((azimuths[1] - azimuths[0]) - 180.0)
+ assert opposed < 5.0, f"faces should oppose, got {azimuths}"
+ for p in planes[:2]:
+ assert p.tilt_deg == pytest.approx(35.0, abs=2.0)
+
+
+def test_two_separate_buildings_on_the_same_plane_are_split():
+ """Two roofs with identical pitch and aspect satisfy one plane equation.
+ RANSAC alone would merge them; the DBSCAN pass is what separates them."""
+ a = make_plane(tilt_deg=30, azimuth_deg=180, origin=(0, 0, 0), seed=4)
+ b = make_plane(tilt_deg=30, azimuth_deg=180, origin=(60, 0, 0), seed=5)
+ planes = segment_roof(np.vstack([a, b]), point_density=8.0, cluster_eps_m=1.5)
+
+ assert len(planes) >= 2, "spatially disjoint faces must not be merged"
+ for p in planes[:2]:
+ assert p.tilt_deg == pytest.approx(30.0, abs=2.0)
+
+
+def test_area_is_the_sloped_area_not_the_footprint():
+ """A 10x8 footprint at 60 degrees covers 80 / cos(60) = 160 m2 of roof."""
+ cloud = make_plane(tilt_deg=60.0, azimuth_deg=180.0, width=10, depth=8, density=10)
+ planes = segment_roof(cloud, point_density=10.0)
+ assert planes
+ assert planes[0].area_m2 == pytest.approx(160.0, rel=0.15)
+
+
+def test_area_of_a_flat_roof_matches_its_footprint():
+ cloud = make_plane(tilt_deg=0.0, azimuth_deg=180.0, width=12, depth=10, density=10)
+ planes = segment_roof(cloud, point_density=10.0)
+ assert planes
+ assert planes[0].area_m2 == pytest.approx(120.0, rel=0.15)
+
+
+def test_survives_realistic_measurement_noise():
+ """1-2 pts/m2 airborne LiDAR carries several centimetres of range noise."""
+ cloud = make_plane(tilt_deg=35.0, azimuth_deg=180.0, noise_m=0.05, density=6)
+ planes = segment_roof(cloud, point_density=6.0, threshold_m=0.3)
+ assert planes
+ assert planes[0].tilt_deg == pytest.approx(35.0, abs=3.0)
+
+
+def test_walls_are_rejected():
+ """A near-vertical surface is a wall or a dormer cheek, not a roof."""
+ cloud = make_plane(tilt_deg=88.0, azimuth_deg=180.0)
+ planes = segment_roof(cloud, point_density=8.0)
+ assert all(p.tilt_deg <= 80.0 for p in planes)
+
+
+def test_noise_alone_yields_nothing_believable():
+ """A handful of scattered returns must not become a roof."""
+ rng = np.random.default_rng(7)
+ cloud = rng.uniform(0, 10, size=(25, 3))
+ assert segment_roof(cloud, point_density=8.0) == []
+
+
+def test_is_deterministic():
+ """Re-running a derive must not shuffle an operator's arrays."""
+ cloud = make_plane(tilt_deg=35, azimuth_deg=180)
+ a = segment_roof(cloud, point_density=8.0)
+ b = segment_roof(cloud, point_density=8.0)
+ assert a == b
+
+
+def test_planes_are_ordered_by_area():
+ big = make_plane(tilt_deg=30, azimuth_deg=180, width=14, depth=12, origin=(0, 0, 0), seed=8)
+ small = make_plane(tilt_deg=30, azimuth_deg=180, width=5, depth=4, origin=(70, 0, 0), seed=9)
+ planes = segment_roof(np.vstack([big, small]), point_density=8.0)
+ assert len(planes) >= 2
+ assert planes[0].area_m2 >= planes[1].area_m2
+
+
+def test_rejects_malformed_input():
+ with pytest.raises(ValueError):
+ segment_roof(np.zeros((10, 2)))
+
+
+def test_kwp_derives_from_area_with_packing_losses():
+ """A 100 m2 face cannot carry 100 m2 of modules."""
+ plane = RoofPlane(tilt_deg=35, azimuth_deg=180, area_m2=100.0, point_count=800, mean_height_m=6.0)
+ kwp = plane.kwp()
+ assert kwp == pytest.approx(100 * 0.70 * 200 / 1000.0)
+ assert kwp < 100 * 200 / 1000.0, "packing factor must reduce the raw area"
+
+
+def test_kwp_scales_with_area():
+ a = RoofPlane(tilt_deg=35, azimuth_deg=180, area_m2=50.0, point_count=1, mean_height_m=1.0)
+ b = RoofPlane(tilt_deg=35, azimuth_deg=180, area_m2=100.0, point_count=1, mean_height_m=1.0)
+ assert b.kwp() == pytest.approx(2 * a.kwp())
+
+
+def test_azimuth_never_reports_360():
+ """A north face fits at ~359.97 deg, which rounds to 360.0 -- the same
+ direction as 0 but an out-of-range-looking value for anything consuming it.
+ Caught by the demo, so pinned here."""
+ cloud = make_plane(tilt_deg=35, azimuth_deg=0, seed=21)
+ planes = segment_roof(cloud, point_density=8.0)
+ assert planes
+ for p in planes:
+ assert 0.0 <= p.azimuth_deg < 360.0, p.azimuth_deg
diff --git a/roofmodel/tests/test_source_formats.py b/roofmodel/tests/test_source_formats.py
new file mode 100644
index 000000000..db41dff69
--- /dev/null
+++ b/roofmodel/tests/test_source_formats.py
@@ -0,0 +1,314 @@
+"""The two Lantmaeteriet products end to end, in the formats they ship in.
+
+Both are STAC APIs over one client and one credential. What differs is the
+payload: *Byggnad Nedladdning, vektor* is GeoPackage, *Laserdata Nedladdning,
+Skog* is LAZ organised as COPC.
+"""
+
+from __future__ import annotations
+
+import json
+import math
+
+import numpy as np
+import pytest
+
+from ftw_roofmodel import pipeline, sweref
+from ftw_roofmodel.buildings import BuildingLookupError, search_buildings
+from ftw_roofmodel.geotorget import (
+ COLLECTION_BUILDINGS,
+ MEDIA_COPC,
+ MEDIA_GEOJSON,
+ MEDIA_GEOPACKAGE,
+ MEDIA_LAZ,
+ Asset,
+ Credentials,
+ StacItem,
+)
+from ftw_roofmodel.pipeline import derive
+from ftw_roofmodel.pointcloud import PointCloudError
+
+from .test_geopackage import build_gpkg, gpkg_blob, wkb_polygon
+
+STOCKHOLM = (59.33, 18.07)
+N, E = sweref.wgs84_to_sweref99tm(*STOCKHOLM)
+
+
+def square(cx, cy, w, d):
+ return [[(cx, cy), (cx + w, cy), (cx + w, cy + d), (cx, cy + d), (cx, cy)]]
+
+
+def roof_face(tilt, azimuth, w, d, origin, density=8, noise=0.04, seed=1):
+ rng = np.random.default_rng(seed)
+ n = int(w * d * density)
+ x = rng.uniform(0, w, n)
+ y = rng.uniform(0, d, n)
+ az = math.radians(azimuth)
+ s = math.tan(math.radians(tilt))
+ z = -(x * math.sin(az) + y * math.cos(az)) * s + rng.normal(0, noise, n)
+ return np.column_stack([x + origin[0], y + origin[1], z + origin[2]])
+
+
+HOUSE = np.vstack([
+ roof_face(35, 180, 12, 6, (E, N, 0), seed=2),
+ roof_face(35, 0, 12, 6, (E, N + 6, 4.2), seed=3),
+])
+
+
+class FakeClient:
+ """A Geotorget stand-in that serves typed assets."""
+
+ def __init__(self, building_asset=None, lidar_asset=None, payloads=None,
+ session=None):
+ self._building_asset = building_asset
+ self._lidar_asset = lidar_asset
+ self._payloads = payloads or {}
+ self.session = session
+ self.downloaded: list[str] = []
+
+ def search(self, collection, bbox, limit=20):
+ if collection == COLLECTION_BUILDINGS:
+ if self._building_asset is None:
+ return []
+ return [StacItem("tile-b", collection, {"data": self._building_asset}, None,
+ raw={"id": "tile-b"})]
+ if self._lidar_asset is None:
+ return []
+ return [StacItem("tile-l", collection, {"data": self._lidar_asset}, None,
+ raw={"id": "tile-l"})]
+
+ def download(self, url):
+ self.downloaded.append(url)
+ return self._payloads[url]
+
+
+def a_geopackage_of_two_buildings():
+ return build_gpkg([
+ ("house-1", "Bostad", gpkg_blob(wkb_polygon(square(E - 6, N - 3, 12, 12)))),
+ ("shed-1", "Komplementbyggnad",
+ gpkg_blob(wkb_polygon(square(E + 40, N + 40, 8, 8)))),
+ ])
+
+
+def test_buildings_come_out_of_a_geopackage_asset():
+ """The normal Byggnad-vektor path: the item points at a .gpkg tile."""
+ url = "https://api.lantmateriet.se/x/byggnad.gpkg"
+ client = FakeClient(
+ building_asset=Asset(url, MEDIA_GEOPACKAGE),
+ payloads={url: a_geopackage_of_two_buildings()},
+ )
+ found = search_buildings(client, latitude=STOCKHOLM[0], longitude=STOCKHOLM[1])
+
+ assert [b.building_id for b in found] == ["house-1", "shed-1"], "nearest first"
+ assert found[0].area_m2 == pytest.approx(144.0, rel=0.01)
+ assert found[0].properties["andamal"] == "Bostad"
+ assert client.downloaded == [url]
+
+
+def test_buildings_come_out_of_a_zipped_geopackage_asset():
+ """The live Lantmaeteriet shape: one item per municipality whose only
+ asset is byggnad_kn.zip with the GeoPackage inside."""
+ import io
+ import zipfile
+
+ buf = io.BytesIO()
+ with zipfile.ZipFile(buf, "w") as zf:
+ zf.writestr("byggnad_kn0180.gpkg", a_geopackage_of_two_buildings())
+ zf.writestr("licens.txt", "CC BY 4.0")
+ url = "https://dl1.lantmateriet.se/byggnadsverk/byggnad_kn0180.zip"
+ client = FakeClient(
+ building_asset=Asset(url, "application/zip"),
+ payloads={url: buf.getvalue()},
+ )
+ found = search_buildings(client, latitude=STOCKHOLM[0], longitude=STOCKHOLM[1])
+
+ assert [b.building_id for b in found] == ["house-1", "shed-1"]
+ assert client.downloaded == [url]
+
+
+def test_a_tile_asset_wins_over_the_tile_outline():
+ """The live municipality items carry their own geometry — the TILE's
+ outline, not a building. With a data asset present, the asset must be
+ read and the outline ignored; treating the outline as a building both
+ invented a giant footprint and skipped the 90k real ones (found live:
+ 'found=0' with no error, because the outline failed the area filter)."""
+ url = "https://dl1.lantmateriet.se/byggnadsverk/byggnad_kn0180.zip"
+ import io
+ import zipfile
+
+ buf = io.BytesIO()
+ with zipfile.ZipFile(buf, "w") as zf:
+ zf.writestr("byggnad_kn0180.gpkg", a_geopackage_of_two_buildings())
+ kommun_outline = {
+ "type": "Polygon",
+ "coordinates": [[[E - 20000, N - 20000], [E + 20000, N - 20000],
+ [E + 20000, N + 20000], [E - 20000, N + 20000],
+ [E - 20000, N - 20000]]],
+ }
+
+ class TileClient(FakeClient):
+ def search(self, collection, bbox, limit=20):
+ return [StacItem("0180", collection, {"data": Asset(url, "application/zip")},
+ None, raw={"id": "0180", "geometry": kommun_outline})]
+
+ client = TileClient(payloads={url: buf.getvalue()})
+ found = search_buildings(client, latitude=STOCKHOLM[0], longitude=STOCKHOLM[1])
+
+ assert [b.building_id for b in found] == ["house-1", "shed-1"]
+ assert client.downloaded == [url]
+
+
+def test_a_municipality_tile_is_filtered_to_the_search_window():
+ """One GeoPackage covers a whole municipality — Stockholm's has 90k+
+ buildings — so rows outside the search window must be dropped, or the
+ reader's row limit truncates the table before it reaches the site.
+
+ Both drop paths matter: a stored header envelope skips the row before
+ its WKB is decoded, and the live Lantmaeteriet files write NO envelopes
+ (flags 0x00), where the parsed geometry's own bounds must do it."""
+ far_e, far_n = E + 5000, N + 5000
+ gpkg = build_gpkg([
+ # Near the site, no envelope: kept via its parsed bounds.
+ ("near-1", "Bostad", gpkg_blob(wkb_polygon(square(E - 6, N - 3, 12, 12)))),
+ # Far away with an envelope: dropped before the WKB is decoded.
+ ("far-1", "Bostad", gpkg_blob(wkb_polygon(square(far_e, far_n, 12, 12)),
+ envelope=(far_e, far_e + 12, far_n, far_n + 12))),
+ # Far away without an envelope: dropped via its parsed bounds.
+ ("far-2", "Bostad", gpkg_blob(wkb_polygon(square(far_e, far_n - 60, 12, 12)))),
+ ])
+ url = "https://api.lantmateriet.se/x/byggnad_kn0180.gpkg"
+ client = FakeClient(building_asset=Asset(url, MEDIA_GEOPACKAGE), payloads={url: gpkg})
+ found = search_buildings(client, latitude=STOCKHOLM[0], longitude=STOCKHOLM[1])
+
+ assert [b.building_id for b in found] == ["near-1"]
+
+
+def test_buildings_come_out_of_a_geojson_asset_too():
+ """Some catalogues publish GeoJSON; both are handled by media type."""
+ url = "https://api.lantmateriet.se/x/byggnad.geojson"
+ doc = {
+ "type": "FeatureCollection",
+ "features": [{
+ "type": "Feature",
+ "id": "gj-1",
+ "geometry": {"type": "Polygon",
+ "coordinates": [[list(p) for p in square(E - 6, N - 3, 12, 12)[0]]]},
+ "properties": {"andamal": "Bostad"},
+ }],
+ }
+ client = FakeClient(
+ building_asset=Asset(url, MEDIA_GEOJSON),
+ payloads={url: json.dumps(doc).encode()},
+ )
+ found = search_buildings(client, latitude=STOCKHOLM[0], longitude=STOCKHOLM[1])
+ assert [b.building_id for b in found] == ["gj-1"]
+
+
+def test_an_unreadable_building_tile_says_what_is_wrong():
+ url = "https://api.lantmateriet.se/x/byggnad.gpkg"
+ client = FakeClient(
+ building_asset=Asset(url, MEDIA_GEOPACKAGE),
+ payloads={url: b"this is not a database"},
+ )
+ with pytest.raises(BuildingLookupError, match="could not be read"):
+ search_buildings(client, latitude=STOCKHOLM[0], longitude=STOCKHOLM[1])
+
+
+def test_an_item_with_no_usable_asset_is_skipped_not_fatal():
+ client = FakeClient(building_asset=Asset("https://x/preview.png", "image/png"))
+ with pytest.raises(BuildingLookupError, match="no building footprints"):
+ search_buildings(client, latitude=STOCKHOLM[0], longitude=STOCKHOLM[1])
+
+
+# --- LiDAR: COPC windowing -------------------------------------------------
+
+
+@pytest.fixture
+def lidar_scene(monkeypatch):
+ """A client whose LiDAR asset is COPC, with both read paths instrumented."""
+ gpkg_url = "https://api.lantmateriet.se/x/byggnad.gpkg"
+ copc_url = "https://api.lantmateriet.se/x/tile.copc.laz"
+ client = FakeClient(
+ building_asset=Asset(gpkg_url, MEDIA_GEOPACKAGE),
+ lidar_asset=Asset(copc_url, MEDIA_COPC),
+ payloads={gpkg_url: a_geopackage_of_two_buildings(), copc_url: b"laz-bytes"},
+ session=object(),
+ )
+ calls: dict[str, object] = {}
+
+ def fake_window(session, url, bounds, timeout=60.0):
+ calls["window"] = {"url": url, "bounds": bounds}
+ return HOUSE
+
+ monkeypatch.setattr(pipeline.pointcloud, "read_copc_window", fake_window)
+ monkeypatch.setattr(pipeline, "load_points", lambda data: HOUSE)
+ return client, calls, copc_url
+
+
+def test_a_picked_building_is_read_as_a_copc_window(lidar_scene):
+ """The payoff: a footprint costs a window, not a 2.5 km tile."""
+ client, calls, copc_url = lidar_scene
+ model = derive(
+ latitude=STOCKHOLM[0], longitude=STOCKHOLM[1],
+ credentials=Credentials("u", "t"), client=client, building_id="house-1",
+ )
+
+ assert model["source"]["fetch"] == "copc-window"
+ assert calls["window"]["url"] == copc_url
+ min_x, min_y, max_x, max_y = calls["window"]["bounds"]
+ # The footprint is 12 m square around the site, plus the eaves buffer.
+ assert min_x == pytest.approx(E - 7, abs=0.5)
+ assert max_x == pytest.approx(E + 7, abs=0.5)
+ assert copc_url not in client.downloaded, "the tile was never downloaded whole"
+
+
+def test_without_a_picked_building_the_tile_is_read_whole(lidar_scene):
+ """No footprint means no window to ask for."""
+ client, calls, copc_url = lidar_scene
+ model = derive(
+ latitude=STOCKHOLM[0], longitude=STOCKHOLM[1],
+ credentials=Credentials("u", "t"), client=client,
+ )
+ assert model["source"]["fetch"] == "whole-tile"
+ assert "window" not in calls
+ assert copc_url in client.downloaded
+
+
+def test_a_host_without_range_support_falls_back_to_the_whole_tile(monkeypatch, lidar_scene):
+ """Best-effort: the operator still gets their roof, just slower."""
+ client, _, copc_url = lidar_scene
+
+ def refuses(session, url, bounds, timeout=60.0):
+ raise PointCloudError("the LiDAR host does not support range requests")
+
+ monkeypatch.setattr(pipeline.pointcloud, "read_copc_window", refuses)
+ model = derive(
+ latitude=STOCKHOLM[0], longitude=STOCKHOLM[1],
+ credentials=Credentials("u", "t"), client=client, building_id="house-1",
+ )
+ assert model["source"]["fetch"] == "whole-tile"
+ assert copc_url in client.downloaded
+ assert model["arrays"], "and the roof still comes out"
+
+
+def test_a_plain_laz_asset_is_read_whole_even_with_a_building(monkeypatch):
+ """Only COPC can be windowed; plain LAZ has no index to seek with."""
+ gpkg_url = "https://api.lantmateriet.se/x/byggnad.gpkg"
+ laz_url = "https://api.lantmateriet.se/x/tile.laz"
+ client = FakeClient(
+ building_asset=Asset(gpkg_url, MEDIA_GEOPACKAGE),
+ lidar_asset=Asset(laz_url, MEDIA_LAZ),
+ payloads={gpkg_url: a_geopackage_of_two_buildings(), laz_url: b"laz-bytes"},
+ session=object(),
+ )
+ monkeypatch.setattr(pipeline, "load_points", lambda data: HOUSE)
+ called = []
+ monkeypatch.setattr(pipeline.pointcloud, "read_copc_window",
+ lambda *a, **k: called.append(1))
+
+ model = derive(
+ latitude=STOCKHOLM[0], longitude=STOCKHOLM[1],
+ credentials=Credentials("u", "t"), client=client, building_id="house-1",
+ )
+ assert model["source"]["fetch"] == "whole-tile"
+ assert not called
diff --git a/roofmodel/tests/test_sweref.py b/roofmodel/tests/test_sweref.py
new file mode 100644
index 000000000..1cd6dfa35
--- /dev/null
+++ b/roofmodel/tests/test_sweref.py
@@ -0,0 +1,148 @@
+"""SWEREF 99 TM projection tests.
+
+The projection has exact analytic properties at the central meridian and the
+equator, which pin the parameters without needing a published coordinate table.
+Round-trip accuracy then pins the series expansion itself.
+"""
+
+import math
+
+import pytest
+
+from ftw_roofmodel.sweref import (
+ bbox_wgs84_to_sweref99tm,
+ metre_box_around,
+ sweref99tm_to_wgs84,
+ wgs84_to_sweref99tm,
+)
+
+
+@pytest.mark.parametrize("lat", [0.0, 55.0, 59.33, 63.0, 69.0])
+def test_central_meridian_maps_to_false_easting(lat):
+ """On the central meridian easting is exactly the false easting, by
+ definition. Any error in scale, ellipsoid or meridian shows up here."""
+ _, easting = wgs84_to_sweref99tm(lat, 15.0)
+ assert easting == pytest.approx(500000.0, abs=1e-6)
+
+
+def test_equator_on_central_meridian_is_the_projection_origin():
+ northing, easting = wgs84_to_sweref99tm(0.0, 15.0)
+ assert northing == pytest.approx(0.0, abs=1e-6)
+ assert easting == pytest.approx(500000.0, abs=1e-6)
+
+
+@pytest.mark.parametrize(
+ "lat,lon",
+ [
+ (55.34, 13.15), # Smygehuk, southernmost Sweden
+ (59.33, 18.07), # Stockholm
+ (57.71, 11.97), # Gothenburg
+ (63.83, 20.26), # Umea
+ (67.86, 20.23), # Kiruna
+ (69.06, 20.55), # Treriksroeset, northernmost
+ ],
+)
+def test_round_trip_is_sub_millimetre(lat, lon):
+ """Project and unproject across the full extent of Sweden."""
+ n, e = wgs84_to_sweref99tm(lat, lon)
+ back_lat, back_lon = sweref99tm_to_wgs84(n, e)
+ # 1e-8 degrees is about 1 mm of latitude.
+ assert back_lat == pytest.approx(lat, abs=1e-8)
+ assert back_lon == pytest.approx(lon, abs=1e-8)
+
+
+def test_coordinates_land_in_the_expected_range_for_sweden():
+ """Sanity-check magnitudes: Swedish SWEREF 99 TM eastings sit inside
+ 260-920 km and northings inside 6100-7700 km. A transposed or
+ wrongly-scaled result would fall far outside."""
+ n, e = wgs84_to_sweref99tm(59.33, 18.07)
+ assert 260_000 < e < 920_000, e
+ assert 6_100_000 < n < 7_700_000, n
+
+
+def test_east_of_the_meridian_increases_easting():
+ _, west = wgs84_to_sweref99tm(59.33, 14.0)
+ _, east = wgs84_to_sweref99tm(59.33, 16.0)
+ assert west < 500000.0 < east
+
+
+def test_north_increases_northing():
+ south, _ = wgs84_to_sweref99tm(55.0, 15.0)
+ north, _ = wgs84_to_sweref99tm(65.0, 15.0)
+ assert north > south
+
+
+def test_one_degree_of_latitude_is_about_111_km():
+ a, _ = wgs84_to_sweref99tm(59.0, 15.0)
+ b, _ = wgs84_to_sweref99tm(60.0, 15.0)
+ # Scaled by k0 = 0.9996 on the central meridian.
+ assert 110_000 < (b - a) < 112_000
+
+
+def test_bbox_uses_all_four_corners():
+ """The projection is not axis-aligned, so the projected box must be at
+ least as large as the one implied by the two diagonal corners."""
+ min_lat, min_lon, max_lat, max_lon = 59.0, 17.0, 60.0, 19.0
+ min_e, min_n, max_e, max_n = bbox_wgs84_to_sweref99tm(
+ min_lat, min_lon, max_lat, max_lon
+ )
+ sw_n, sw_e = wgs84_to_sweref99tm(min_lat, min_lon)
+ ne_n, ne_e = wgs84_to_sweref99tm(max_lat, max_lon)
+ assert min_e <= sw_e and min_n <= sw_n
+ assert max_e >= ne_e and max_n >= ne_n
+ assert min_e < max_e and min_n < max_n
+
+
+def test_metre_box_is_square_on_the_ground():
+ """A 100 m box must measure 200 m on both axes regardless of latitude --
+ the whole reason it is built in projected metres rather than in degrees."""
+ for lat in (55.5, 59.33, 68.0):
+ south, west, north, east = metre_box_around(lat, 15.0, 100.0)
+ sn, se = wgs84_to_sweref99tm(south, west)
+ nn, ne = wgs84_to_sweref99tm(north, east)
+ assert (ne - se) == pytest.approx(200.0, abs=0.5)
+ assert (nn - sn) == pytest.approx(200.0, abs=0.5)
+
+
+def test_metre_box_brackets_its_centre():
+ south, west, north, east = metre_box_around(59.33, 18.07, 50.0)
+ assert south < 59.33 < north
+ assert west < 18.07 < east
+
+
+def test_stac_search_bbox_sweref_matches_the_long_form():
+ from ftw_roofmodel.sweref import bbox_wgs84_to_sweref99tm, stac_search_bbox
+
+ south, west, north, east = metre_box_around(59.33, 18.07, 40.0)
+ assert stac_search_bbox(59.33, 18.07, 40.0) == bbox_wgs84_to_sweref99tm(
+ south, west, north, east
+ )
+
+
+def test_stac_search_bbox_wgs84_is_lon_lat_ordered():
+ """The STAC spec's bbox is [west, south, east, north] in degrees."""
+ from ftw_roofmodel.sweref import stac_search_bbox
+
+ west, south, east, north = stac_search_bbox(59.33, 18.07, 40.0, bbox_epsg=4326)
+ assert west < 18.07 < east
+ assert south < 59.33 < north
+
+
+def test_stac_search_bbox_refuses_a_crs_it_cannot_produce():
+ from ftw_roofmodel.sweref import stac_search_bbox
+
+ with pytest.raises(ValueError):
+ stac_search_bbox(59.33, 18.07, 40.0, bbox_epsg=3857)
+
+
+def test_degree_box_would_have_been_wrong_at_high_latitude():
+ """Guards the reason metre_box_around exists: a fixed degree offset gives
+ wildly different ground distances at Malmoe and Kiruna, so anyone tempted to
+ simplify it back to degrees has to defeat this test first."""
+ span = []
+ for lat in (55.5, 68.0):
+ # What a naive 0.001-degree longitude offset would span, in metres.
+ _, e0 = wgs84_to_sweref99tm(lat, 15.0)
+ _, e1 = wgs84_to_sweref99tm(lat, 15.001)
+ span.append(e1 - e0)
+ assert span[0] / span[1] > 1.5, span
diff --git a/web/components/ftw-bar-chart.js b/web/components/ftw-bar-chart.js
index d7465a6ba..6f0144d8e 100644
--- a/web/components/ftw-bar-chart.js
+++ b/web/components/ftw-bar-chart.js
@@ -188,6 +188,26 @@ class FtwBarChart extends FtwElement {
opacity: 0.85;
pointer-events: none;
}
+ /* Optional dashed overlay line (e.g. STRÅNG-expected PV vs the
+ produced bars). An SVG stretched to fill .bar-area, plotted in a
+ 0..100 viewBox so it needs no pixel math against the CSS grid;
+ non-scaling-stroke keeps the dash crisp despite the stretch. */
+ .overlay-line {
+ position: absolute;
+ inset: 0;
+ width: 100%;
+ height: 100%;
+ pointer-events: none;
+ overflow: visible;
+ }
+ .overlay-line polyline {
+ fill: none;
+ stroke: var(--ftw-overlay-color, var(--amber, #f59e0b));
+ stroke-width: 1.5;
+ stroke-dasharray: 4 3;
+ vector-effect: non-scaling-stroke;
+ opacity: 0.9;
+ }
`;
static get observedAttributes() {
@@ -197,6 +217,7 @@ class FtwBarChart extends FtwElement {
constructor() {
super();
this._data = [];
+ this._overlay = null;
}
attributeChangedCallback() { this.update(); }
@@ -209,6 +230,16 @@ class FtwBarChart extends FtwElement {
}
get data() { return this._data; }
+ // Optional dashed overlay line, index-aligned with .data. Shape:
+ // { values: (number|null)[], color?: string }
+ // values[i] is plotted above column i on the SAME axis as the bars
+ // (nulls / non-finite entries break the line). Set null to remove.
+ set overlay(o) {
+ this._overlay = o && Array.isArray(o.values) ? o : null;
+ this.update();
+ }
+ get overlay() { return this._overlay; }
+
render() {
const accent = this.getAttribute("accent");
const height = this.getAttribute("chart-height");
@@ -262,6 +293,17 @@ class FtwBarChart extends FtwElement {
}
const avg = count > 0 ? sum / count : 0;
+ // Fold overlay values into the axis max so the expected line and the
+ // produced bars share one scale (an expected level above the tallest
+ // bar must still fit in-frame).
+ const overlayVals = this._overlay ? this._overlay.values : null;
+ if (overlayVals) {
+ for (const ov of overlayVals) {
+ const n = Number(ov);
+ if (isFinite(n) && n > max) max = n;
+ }
+ }
+
const colsSvg = this._data.map((d) => {
const v = Number(d.value) || 0;
// 2% floor keeps tiny-but-nonzero values visible; gate on v>0 so
@@ -300,10 +342,34 @@ class FtwBarChart extends FtwElement {
`title="average ${display}">`;
}
+ // Dashed overlay line (expected series). Plotted in a 0..100 viewBox
+ // stretched to fill .bar-area: x centers each column, y is the value
+ // as a percentage of the shared max (inverted — SVG y grows downward).
+ let overlayLine = "";
+ if (overlayVals && max > 0) {
+ const n = this._data.length;
+ const pts = [];
+ for (let i = 0; i < n; i++) {
+ const val = Number(overlayVals[i]);
+ if (!isFinite(val)) continue;
+ const x = n > 1 ? (i + 0.5) / n * 100 : 50;
+ const y = Math.max(0, Math.min(100, 100 - (val / max) * 100));
+ pts.push(`${x.toFixed(2)},${y.toFixed(2)}`);
+ }
+ if (pts.length >= 2) {
+ const color = this._overlay.color;
+ if (color) this.style.setProperty("--ftw-overlay-color", color);
+ overlayLine =
+ ``;
+ }
+ }
+
return `
${lblsSvg}
diff --git a/web/components/ftw-history-card.js b/web/components/ftw-history-card.js
index ed659092e..1e1393c0c 100644
--- a/web/components/ftw-history-card.js
+++ b/web/components/ftw-history-card.js
@@ -28,7 +28,7 @@
import { FtwElement, ftwDebugDelay } from "./ftw-element.js";
import { apiFetch } from "./api-fetch.js";
-import "./ftw-bar-chart.js";
+import "./ftw-bar-chart.js?v=strang1";
const FIELD_BY_METRIC = {
import: "import_wh",
@@ -82,6 +82,36 @@ function fetchDailyEnergy(days) {
return promise;
}
+// STRÅNG-based PV performance (expected-vs-actual). Only the "Produced"
+// (metric="pv") tile fetches this, to overlay the weather-expected line.
+// Tolerant: a disabled/absent service or any error resolves to
+// { enabled:false } so the bars still render without an overlay.
+const pvPerfFetchCache = new Map(); // days -> { at, data?, promise? }
+
+function fetchPVPerformance(days) {
+ const now = Date.now();
+ const cached = pvPerfFetchCache.get(days);
+ if (cached && cached.data && now - cached.at < DAILY_CACHE_TTL_MS) {
+ return Promise.resolve(cached.data);
+ }
+ if (cached && cached.promise && now - cached.at < DAILY_CACHE_TTL_MS) {
+ return cached.promise;
+ }
+ const promise = apiFetch("/api/pv/performance?days=" + days)
+ .then((r) => (r.ok ? r.json() : { enabled: false }))
+ .then((resp) => {
+ const data = resp || { enabled: false };
+ pvPerfFetchCache.set(days, { at: Date.now(), data });
+ return data;
+ })
+ .catch(() => {
+ pvPerfFetchCache.delete(days);
+ return { enabled: false };
+ });
+ pvPerfFetchCache.set(days, { at: now, promise });
+ return promise;
+}
+
class FtwHistoryCard extends FtwElement {
static styles = `
:host { display: block; }
@@ -205,6 +235,27 @@ class FtwHistoryCard extends FtwElement {
margin-left: 6px;
letter-spacing: 0;
}
+ /* STRÅNG expected-vs-actual caption under the Produced chart. Hidden
+ until the pv tile has a performance overlay to describe. The dashed
+ swatch mirrors the overlay line so the legend reads at a glance. */
+ .strang-note {
+ margin-top: 6px;
+ font-size: 0.72rem;
+ color: var(--fg-muted);
+ font-family: var(--mono);
+ display: flex;
+ align-items: center;
+ gap: 6px;
+ }
+ .strang-note[hidden] { display: none; }
+ .strang-note .swatch {
+ display: inline-block;
+ width: 16px;
+ height: 0;
+ border-top: 2px dashed #fcd34d;
+ flex: 0 0 auto;
+ }
+ .strang-note .pr { color: var(--fg-label); font-weight: 600; }
@media (max-width: 900px) {
.card-inner { padding: var(--card-pad-tight, 12px 14px); }
}
@@ -302,6 +353,7 @@ class FtwHistoryCard extends FtwElement {
— kWh
— tiny 3D preview of a site's PV-array config.
//
-// Purpose: the settings Weather tab lets the operator list each PV
+// Purpose: the settings Control tab lets the operator list each PV
// plane (name + rated_w + tilt_deg + azimuth_deg). Those four numbers
// fully specify a panel's orientation but they are hard to cross-
// check by eye on a phone in the shed. This component turns the
@@ -12,7 +12,7 @@
//
// The component is self-contained: Three.js loads via the importmap
// declared in index.html ("three" + "three/addons/"); this file is
-// lazy-imported from settings.js the first time the Weather tab
+// lazy-imported from settings.js the first time the Control tab
// opens, so the dashboard's main thread never pays for three.js on
// pages that don't touch the settings modal.
//
diff --git a/web/components/index.js b/web/components/index.js
index 38c87bd81..b104323e4 100644
--- a/web/components/index.js
+++ b/web/components/index.js
@@ -23,8 +23,8 @@ import "./ftw-battery-control.js?v=apifetch1";
import "./ftw-pv-control.js?v=apifetch1";
import "./ftw-price-chart.js?v=zones1";
import "./ftw-energy-cake.js";
-import "./ftw-bar-chart.js";
-import "./ftw-history-card.js?v=apiread2";
+import "./ftw-bar-chart.js?v=strang1";
+import "./ftw-history-card.js?v=strang2";
import "./ftw-savings-card.js?v=zones1";
import "./ftw-update-check.js?v=apifetch1";
import "./ftw-notif-status.js?v=apifetch1";
diff --git a/web/components/pv-array-geometry.js b/web/components/pv-array-geometry.js
new file mode 100644
index 000000000..110f903e6
--- /dev/null
+++ b/web/components/pv-array-geometry.js
@@ -0,0 +1,208 @@
+// Turning a rectangle drawn on the map into a PV array.
+//
+// Drawing supplies two of the three numbers an array needs. Area comes from
+// the shape; azimuth from how the shape is turned. Tilt cannot be seen from
+// directly above at all, so it stays typed — and it is also what converts the
+// drawn outline into real panel area, because what you trace on a map is the
+// *horizontal projection* of a sloped rectangle, not the rectangle itself.
+//
+// Everything here is pure so it can be tested without a browser or a map.
+
+export const DEFAULT_PACKING_FACTOR = 0.7;
+export const DEFAULT_MODULE_W_PER_M2 = 200;
+export const DEFAULT_TILT_DEG = 35;
+
+// IUGG mean Earth radius. At the scale of one roof the radius matters far
+// less than the flat-Earth approximation below, which is exact enough for a
+// 20 m rectangle and meaningless across a county.
+const EARTH_RADIUS_M = 6371008.8;
+const DEG = Math.PI / 180;
+
+// Beyond this a roof is a wall: cos(tilt) approaches zero and the plan-area
+// division runs away. A wall has no horizontal projection to trace anyway, so
+// clamping here bounds the arithmetic instead of returning Infinity.
+const MAX_TILT_FOR_PROJECTION_DEG = 89;
+
+function stripClosingVertex(ring) {
+ if (ring.length > 1) {
+ const first = ring[0];
+ const last = ring[ring.length - 1];
+ if (first[0] === last[0] && first[1] === last[1]) return ring.slice(0, -1);
+ }
+ return ring.slice();
+}
+
+/**
+ * Project a WGS84 ring ([[lon, lat], …]) to metres about its own centroid.
+ *
+ * A local tangent plane, not a real projection: over one building the error
+ * is well under the precision anyone draws with, and it avoids carrying a
+ * projection library into the settings page.
+ */
+export function toLocalMetres(ring) {
+ const pts = stripClosingVertex(ring || []);
+ if (pts.length === 0) return [];
+ let lon0 = 0;
+ let lat0 = 0;
+ for (const [lon, lat] of pts) {
+ lon0 += lon;
+ lat0 += lat;
+ }
+ lon0 /= pts.length;
+ lat0 /= pts.length;
+ const mPerDegLat = EARTH_RADIUS_M * DEG;
+ const mPerDegLon = mPerDegLat * Math.cos(lat0 * DEG);
+ return pts.map(([lon, lat]) => [(lon - lon0) * mPerDegLon, (lat - lat0) * mPerDegLat]);
+}
+
+/** Area of the drawn outline in m², as seen from above. */
+export function planAreaM2(ring) {
+ const p = toLocalMetres(ring);
+ if (p.length < 3) return 0;
+ let twiceArea = 0;
+ for (let i = 0; i < p.length; i++) {
+ const [x1, y1] = p[i];
+ const [x2, y2] = p[(i + 1) % p.length];
+ twiceArea += x1 * y2 - x2 * y1;
+ }
+ return Math.abs(twiceArea) / 2;
+}
+
+/** Compass bearing of a local vector, 0 = north, 90 = east. */
+function bearingDeg(dx, dy) {
+ return (((Math.atan2(dx, dy) / DEG) % 360) + 360) % 360;
+}
+
+/**
+ * Direction of the ring's longest edge, as a line in [0, 180).
+ *
+ * For a panel rectangle that edge runs along the ridge, which is a line and
+ * not an arrow — calling it "north" rather than "south" would be a
+ * distinction the drawing does not contain.
+ */
+export function ridgeAzimuthDeg(ring) {
+ const p = toLocalMetres(ring);
+ if (p.length < 2) return null;
+ let longest = 0;
+ let bx = 0;
+ let by = 0;
+ for (let i = 0; i < p.length; i++) {
+ const [x1, y1] = p[i];
+ const [x2, y2] = p[(i + 1) % p.length];
+ const dx = x2 - x1;
+ const dy = y2 - y1;
+ const len = Math.hypot(dx, dy);
+ if (len > longest) {
+ longest = len;
+ bx = dx;
+ by = dy;
+ }
+ }
+ if (longest <= 0) return null;
+ return bearingDeg(bx, by) % 180;
+}
+
+/** Shortest angle between two compass bearings, in degrees. */
+export function angularDistanceDeg(a, b) {
+ const d = Math.abs((((a - b) % 360) + 360) % 360);
+ return d > 180 ? 360 - d : d;
+}
+
+/** Turn an azimuth to face the opposite way. */
+export function flipAzimuthDeg(azimuthDeg) {
+ return ((((azimuthDeg + 180) % 360) + 360) % 360);
+}
+
+/**
+ * The two directions the face could point: perpendicular to the ridge, either
+ * side of it. A flat outline genuinely does not say which.
+ */
+export function faceAzimuthCandidates(ring) {
+ const ridge = ridgeAzimuthDeg(ring);
+ if (ridge === null) return [];
+ return [(ridge + 90) % 360, (ridge + 270) % 360];
+}
+
+/**
+ * The candidate a panel is more likely to use: the equatorward one.
+ *
+ * This is a default, not a measurement. Both perpendiculars fit the drawing
+ * equally well, so the UI offers a flip rather than pretending to know.
+ */
+export function preferredAzimuthDeg(ring, latitudeDeg) {
+ const candidates = faceAzimuthCandidates(ring);
+ if (candidates.length === 0) return null;
+ const target = (latitudeDeg || 0) >= 0 ? 180 : 0;
+ const [a, b] = candidates;
+ return angularDistanceDeg(a, target) <= angularDistanceDeg(b, target) ? a : b;
+}
+
+/**
+ * Real panel area from the traced outline.
+ *
+ * A sloped rectangle of area A casts a shadow of A·cos(tilt) on the map, so
+ * recovering it divides that back out. A 35° roof carries about 22 % more
+ * panel than its outline suggests, which is the difference between a
+ * believable rating and a quietly low one.
+ */
+export function slopeAreaM2(planArea, tiltDeg) {
+ const tilt = Math.min(Math.max(tiltDeg || 0, 0), MAX_TILT_FOR_PROJECTION_DEG);
+ return planArea / Math.cos(tilt * DEG);
+}
+
+/** Installable DC capacity in watts for a roof area, matching the roof model's basis. */
+export function ratedWFromSlopeArea(areaM2, packingFactor, moduleWPerM2) {
+ const packing = packingFactor == null ? DEFAULT_PACKING_FACTOR : packingFactor;
+ const wPerM2 = moduleWPerM2 == null ? DEFAULT_MODULE_W_PER_M2 : moduleWPerM2;
+ return areaM2 * packing * wPerM2;
+}
+
+/** Human-readable face name, mirroring the roof model's naming. */
+export function compassName(azimuthDeg, tiltDeg) {
+ if (tiltDeg < 5) return "Roof flat";
+ const points = [
+ [0, "north"], [45, "north-east"], [90, "east"], [135, "south-east"],
+ [180, "south"], [225, "south-west"], [270, "west"], [315, "north-west"],
+ [360, "north"],
+ ];
+ let best = points[0];
+ for (const p of points) {
+ if (Math.abs(p[0] - azimuthDeg) < Math.abs(best[0] - azimuthDeg)) best = p;
+ }
+ return `Roof ${best[1]}`;
+}
+
+function round(value, places) {
+ const factor = 10 ** places;
+ return Math.round(value * factor) / factor;
+}
+
+/**
+ * Everything a drawn rectangle says about one array.
+ *
+ * Returns the config-shaped entry separately from the measurements, so only
+ * the four fields weather.pv_arrays actually defines are ever written back.
+ */
+export function arrayFromRing(ring, options) {
+ const opts = options || {};
+ const plan = planAreaM2(ring);
+ if (!(plan > 0)) return null;
+ const tiltDeg = opts.tiltDeg == null ? DEFAULT_TILT_DEG : opts.tiltDeg;
+ const candidates = faceAzimuthCandidates(ring);
+ const azimuth = opts.azimuthDeg == null
+ ? preferredAzimuthDeg(ring, opts.latitude)
+ : opts.azimuthDeg;
+ const azimuthDeg = azimuth == null ? 180 : Math.round(azimuth);
+ const slope = slopeAreaM2(plan, tiltDeg);
+ return {
+ array: {
+ name: opts.name || compassName(azimuthDeg, tiltDeg),
+ rated_w: Math.round(ratedWFromSlopeArea(slope, opts.packingFactor, opts.moduleWPerM2)),
+ tilt_deg: tiltDeg,
+ azimuth_deg: azimuthDeg,
+ },
+ planAreaM2: round(plan, 1),
+ slopeAreaM2: round(slope, 1),
+ azimuthCandidates: candidates.map((c) => Math.round(c)),
+ };
+}
diff --git a/web/components/pv-array-geometry.test.mjs b/web/components/pv-array-geometry.test.mjs
new file mode 100644
index 000000000..56b112cfb
--- /dev/null
+++ b/web/components/pv-array-geometry.test.mjs
@@ -0,0 +1,203 @@
+// node --test web/components/pv-array-geometry.test.mjs
+
+import assert from "node:assert/strict";
+import { describe, it } from "node:test";
+
+import {
+ DEFAULT_MODULE_W_PER_M2,
+ DEFAULT_PACKING_FACTOR,
+ angularDistanceDeg,
+ arrayFromRing,
+ compassName,
+ faceAzimuthCandidates,
+ flipAzimuthDeg,
+ ratedWFromSlopeArea,
+ planAreaM2,
+ preferredAzimuthDeg,
+ ridgeAzimuthDeg,
+ slopeAreaM2,
+} from "./pv-array-geometry.js";
+
+// Fixtures are built with the standard ellipsoidal metres-per-degree series,
+// deliberately *not* with the module's own spherical projection — a shape
+// round-tripped through the code under test would agree with itself and prove
+// nothing. The two disagree by roughly half a percent in area at this
+// latitude, which is the known bias of a spherical Earth against WGS84 and is
+// far below the precision anyone draws a roof with.
+const STOCKHOLM = { lat: 59.3293, lon: 18.0686 };
+
+function metresPerDegree(latDeg) {
+ const p = (latDeg * Math.PI) / 180;
+ return {
+ lat: 111132.92 - 559.82 * Math.cos(2 * p) + 1.175 * Math.cos(4 * p)
+ - 0.0023 * Math.cos(6 * p),
+ lon: 111412.84 * Math.cos(p) - 93.5 * Math.cos(3 * p) + 0.118 * Math.cos(5 * p),
+ };
+}
+
+/** Build a WGS84 ring from local east/north offsets in metres. */
+function ringFromMetres(offsets, origin) {
+ const m = metresPerDegree(origin.lat);
+ return offsets.map(([x, y]) => [origin.lon + x / m.lon, origin.lat + y / m.lat]);
+}
+
+/** Rotate local offsets counter-clockwise in the east/north plane. */
+function rotate(offsets, degrees) {
+ const r = (degrees * Math.PI) / 180;
+ const c = Math.cos(r);
+ const s = Math.sin(r);
+ return offsets.map(([x, y]) => [x * c - y * s, x * s + y * c]);
+}
+
+// 10 m along the ridge (east-west) by 6 m down the slope.
+const RECT_10x6 = [[-5, -3], [5, -3], [5, 3], [-5, 3]];
+
+describe("plan area", () => {
+ it("recovers the drawn size in square metres", () => {
+ const ring = ringFromMetres(RECT_10x6, STOCKHOLM);
+ const area = planAreaM2(ring);
+ assert.ok(Math.abs(area - 60) < 0.6, `area ${area} should be ~60 m²`);
+ });
+
+ it("does not care whether the ring repeats its first point", () => {
+ const ring = ringFromMetres(RECT_10x6, STOCKHOLM);
+ const closed = [...ring, ring[0]];
+ assert.ok(Math.abs(planAreaM2(ring) - planAreaM2(closed)) < 1e-9);
+ });
+
+ it("is unsigned, so winding order cannot produce a negative roof", () => {
+ const ring = ringFromMetres(RECT_10x6, STOCKHOLM);
+ assert.ok(Math.abs(planAreaM2(ring) - planAreaM2([...ring].reverse())) < 1e-9);
+ });
+
+ it("treats a shape with no area as no array", () => {
+ assert.equal(planAreaM2([]), 0);
+ assert.equal(planAreaM2([[18, 59], [18.001, 59]]), 0);
+ assert.equal(arrayFromRing([[18, 59], [18.001, 59]], {}), null);
+ });
+});
+
+describe("orientation", () => {
+ it("reads the ridge from the longest edge", () => {
+ const ring = ringFromMetres(RECT_10x6, STOCKHOLM);
+ // The 10 m edges run east-west: a ridge bearing of 90°.
+ assert.ok(Math.abs(ridgeAzimuthDeg(ring) - 90) < 0.5);
+ });
+
+ it("offers both faces the outline permits, and no others", () => {
+ const ring = ringFromMetres(RECT_10x6, STOCKHOLM);
+ const [a, b] = faceAzimuthCandidates(ring).map(Math.round);
+ assert.deepEqual([a, b].sort((x, y) => x - y), [0, 180]);
+ });
+
+ it("defaults to the equatorward face, per hemisphere", () => {
+ const north = ringFromMetres(RECT_10x6, STOCKHOLM);
+ assert.ok(Math.abs(preferredAzimuthDeg(north, STOCKHOLM.lat) - 180) < 0.5);
+
+ const south = ringFromMetres(RECT_10x6, { lat: -33.87, lon: 151.21 });
+ const picked = preferredAzimuthDeg(south, -33.87);
+ assert.ok(angularDistanceDeg(picked, 0) < 0.5, `expected ~0°, got ${picked}`);
+ });
+
+ it("follows the rectangle round as it turns", () => {
+ // Turning the shape 30° counter-clockwise swings the ridge from 90° to
+ // 60°, so the faces move with it: 150° and 330°, and south-ish wins.
+ const ring = ringFromMetres(rotate(RECT_10x6, 30), STOCKHOLM);
+ assert.ok(Math.abs(ridgeAzimuthDeg(ring) - 60) < 0.5);
+ assert.ok(Math.abs(preferredAzimuthDeg(ring, STOCKHOLM.lat) - 150) < 0.5);
+ });
+
+ it("keeps the ridge a line rather than an arrow", () => {
+ // Drawing the same rectangle the other way round is the same roof.
+ const ring = ringFromMetres(RECT_10x6, STOCKHOLM);
+ const reversed = ringFromMetres([...RECT_10x6].reverse(), STOCKHOLM);
+ assert.ok(Math.abs(ridgeAzimuthDeg(ring) - ridgeAzimuthDeg(reversed)) < 0.5);
+ });
+
+ it("flips to the opposite face", () => {
+ assert.equal(flipAzimuthDeg(180), 0);
+ assert.equal(flipAzimuthDeg(0), 180);
+ assert.equal(flipAzimuthDeg(270), 90);
+ assert.equal(flipAzimuthDeg(350), 170);
+ });
+
+ it("measures the shorter way round the compass", () => {
+ assert.equal(angularDistanceDeg(350, 10), 20);
+ assert.equal(angularDistanceDeg(10, 350), 20);
+ assert.equal(angularDistanceDeg(0, 180), 180);
+ });
+});
+
+describe("tilt turns an outline into panel area", () => {
+ it("leaves a flat roof alone", () => {
+ assert.ok(Math.abs(slopeAreaM2(60, 0) - 60) < 1e-9);
+ });
+
+ it("recovers the area hidden by the slope", () => {
+ // cos 60° = 0.5, so a 60 m² shadow is cast by 120 m² of roof.
+ assert.ok(Math.abs(slopeAreaM2(60, 60) - 120) < 1e-9);
+ // A 35° roof carries ~22 % more panel than its outline suggests.
+ assert.ok(Math.abs(slopeAreaM2(60, 35) - 73.24) < 0.05);
+ });
+
+ it("stays finite at a wall, where there is no outline to trace", () => {
+ assert.ok(Number.isFinite(slopeAreaM2(60, 90)));
+ });
+
+ it("means a steeper roof is a bigger array for the same drawing", () => {
+ const ring = ringFromMetres(RECT_10x6, STOCKHOLM);
+ const flat = arrayFromRing(ring, { latitude: STOCKHOLM.lat, tiltDeg: 0 });
+ const steep = arrayFromRing(ring, { latitude: STOCKHOLM.lat, tiltDeg: 45 });
+ assert.ok(steep.array.rated_w > flat.array.rated_w,
+ `${steep.array.rated_w} should exceed ${flat.array.rated_w}`);
+ });
+});
+
+describe("capacity", () => {
+ it("uses the same basis as the roof model", () => {
+ // 60 m² × 0.70 packing × 200 W/m² = 8400 W.
+ assert.ok(Math.abs(ratedWFromSlopeArea(60) - 8400) < 1e-9);
+ assert.equal(DEFAULT_PACKING_FACTOR, 0.7);
+ assert.equal(DEFAULT_MODULE_W_PER_M2, 200);
+ });
+
+ it("honours an overridden packing factor", () => {
+ assert.ok(Math.abs(ratedWFromSlopeArea(60, 0.5, 200) - 6000) < 1e-9);
+ });
+});
+
+describe("the entry written back to config", () => {
+ it("carries only the four fields weather.pv_arrays defines", () => {
+ const ring = ringFromMetres(RECT_10x6, STOCKHOLM);
+ const out = arrayFromRing(ring, { latitude: STOCKHOLM.lat });
+ assert.deepEqual(
+ Object.keys(out.array).sort(),
+ ["azimuth_deg", "name", "rated_w", "tilt_deg"],
+ );
+ });
+
+ it("describes a south-facing 35° roof from the drawing alone", () => {
+ const ring = ringFromMetres(RECT_10x6, STOCKHOLM);
+ const out = arrayFromRing(ring, { latitude: STOCKHOLM.lat });
+ assert.equal(out.array.azimuth_deg, 180);
+ assert.equal(out.array.tilt_deg, 35);
+ assert.equal(out.array.name, "Roof south");
+ assert.ok(Math.abs(out.planAreaM2 - 60) < 0.6);
+ assert.ok(Math.abs(out.slopeAreaM2 - 73.2) < 0.6);
+ assert.ok(Math.abs(out.array.rated_w - 10250) < 100);
+ assert.deepEqual(out.azimuthCandidates.slice().sort((a, b) => a - b), [0, 180]);
+ });
+
+ it("lets an explicit azimuth override the guess", () => {
+ const ring = ringFromMetres(RECT_10x6, STOCKHOLM);
+ const out = arrayFromRing(ring, { latitude: STOCKHOLM.lat, azimuthDeg: 0 });
+ assert.equal(out.array.azimuth_deg, 0);
+ assert.equal(out.array.name, "Roof north");
+ });
+
+ it("names a flat roof for what it is", () => {
+ assert.equal(compassName(180, 0), "Roof flat");
+ assert.equal(compassName(90, 35), "Roof east");
+ assert.equal(compassName(225, 35), "Roof south-west");
+ });
+});
diff --git a/web/index.html b/web/index.html
index 1e46ac6d0..944951156 100644
--- a/web/index.html
+++ b/web/index.html
@@ -5,10 +5,10 @@
FTW
-
+
@@ -97,7 +97,6 @@ Settings
-
@@ -990,14 +989,14 @@ Price bars (top of the chart)
-
+
-
+
-
+
diff --git a/web/leaflet-vendor.test.mjs b/web/leaflet-vendor.test.mjs
deleted file mode 100644
index ab29ad2db..000000000
--- a/web/leaflet-vendor.test.mjs
+++ /dev/null
@@ -1,36 +0,0 @@
-import assert from 'node:assert/strict';
-import { existsSync, readFileSync } from 'node:fs';
-import { dirname, join } from 'node:path';
-import test from 'node:test';
-import { fileURLToPath } from 'node:url';
-
-const webRoot = dirname(fileURLToPath(import.meta.url));
-const weather = readFileSync(join(webRoot, 'settings', 'tabs', 'weather.js'), 'utf8');
-const vendor = join(webRoot, 'vendor', 'leaflet');
-
-test('weather map loads Leaflet from the vendored copy, not unpkg', () => {
- assert.doesNotMatch(weather, /unpkg\.com/);
- assert.match(weather, /\/vendor\/leaflet\/leaflet\.js/);
- assert.match(weather, /\/vendor\/leaflet\/leaflet\.css/);
-});
-
-test('weather map sends a Referer on OSM tiles so volunteer servers do not 403', () => {
- assert.match(weather, /tile\.openstreetmap\.org/);
- assert.match(weather, /referrerPolicy:\s*"strict-origin-when-cross-origin"/);
- const leaflet = readFileSync(join(vendor, 'leaflet.js'), 'utf8');
- assert.match(leaflet, /typeof this\.options\.referrerPolicy/);
-});
-
-test('vendored Leaflet 1.9.4 files are present', () => {
- for (const rel of [
- 'leaflet.js',
- 'leaflet.css',
- 'LICENSE',
- 'README.md',
- 'images/marker-icon.png',
- 'images/marker-icon-2x.png',
- 'images/marker-shadow.png',
- ]) {
- assert.ok(existsSync(join(vendor, rel)), rel + ' must be vendored');
- }
-});
diff --git a/web/maplibre-vendor.test.mjs b/web/maplibre-vendor.test.mjs
new file mode 100644
index 000000000..761a6e8a7
--- /dev/null
+++ b/web/maplibre-vendor.test.mjs
@@ -0,0 +1,44 @@
+import assert from 'node:assert/strict';
+import { existsSync, readFileSync } from 'node:fs';
+import { dirname, join } from 'node:path';
+import test from 'node:test';
+import { fileURLToPath } from 'node:url';
+
+const webRoot = dirname(fileURLToPath(import.meta.url));
+const weather = readFileSync(join(webRoot, 'settings', 'tabs', 'weather.js'), 'utf8');
+const vendor = join(webRoot, 'vendor', 'maplibre');
+
+test('weather map loads MapLibre from the vendored copy, not a CDN', () => {
+ assert.doesNotMatch(weather, /unpkg\.com/);
+ assert.doesNotMatch(weather, /cdn\.jsdelivr\.net/);
+ assert.match(weather, /\/vendor\/maplibre\//);
+ assert.match(weather, /maplibre-gl\.mjs/);
+ assert.match(weather, /maplibre-gl\.css/);
+});
+
+test('weather map sends a Referer on OSM tiles so volunteer servers do not 403', () => {
+ assert.match(weather, /tile\.openstreetmap\.org/);
+ assert.match(weather, /referrerPolicy:\s*"strict-origin-when-cross-origin"/);
+ // The vendored build must actually honor a per-request referrerPolicy from
+ // transformRequest, or the opt-in above is a no-op.
+ const shared = readFileSync(join(vendor, 'maplibre-gl-shared.mjs'), 'utf8');
+ assert.match(shared, /referrerPolicy/);
+});
+
+test('vendored MapLibre GL JS 6.7.0 files are present', () => {
+ for (const rel of [
+ 'maplibre-gl.mjs',
+ 'maplibre-gl-shared.mjs',
+ 'maplibre-gl-worker.mjs',
+ 'maplibre-gl.css',
+ 'LICENSE.txt',
+ 'README.md',
+ ]) {
+ assert.ok(existsSync(join(vendor, rel)), rel + ' must be vendored');
+ }
+ // v6 is code-split: the entry imports its chunks by relative URL, so the
+ // names above are load-bearing, not just files that happen to exist.
+ const entry = readFileSync(join(vendor, 'maplibre-gl.mjs'), 'utf8');
+ assert.match(entry, /maplibre-gl-shared\.mjs/);
+ assert.match(entry, /maplibre-gl-worker\.mjs/);
+});
diff --git a/web/settings.js b/web/settings.js
index cce548233..b3b29ef65 100644
--- a/web/settings.js
+++ b/web/settings.js
@@ -65,6 +65,24 @@
if (e.target === modal) modal.classList.add("hidden");
});
+ // The help bubble is pure CSS and normally opens downward, but the modal
+ // body is the scroll container that clips it: a badge in the last section
+ // has less room below it than a long help text needs. Flip the bubble
+ // upward when the badge sits in the lower part of the visible modal, so it
+ // grows into space that exists. Measured on hover because scroll position,
+ // not the badge, decides which way is open.
+ modal.addEventListener("mouseover", function (e) {
+ var badge = e.target && e.target.closest && e.target.closest(".help");
+ if (!badge) return;
+ var scroller = badge.closest(".modal-body");
+ if (!scroller) return;
+ var box = scroller.getBoundingClientRect();
+ var at = badge.getBoundingClientRect();
+ badge.classList.toggle(
+ "help-up", at.top + at.height / 2 > box.top + box.height * 0.55
+ );
+ });
+
tabsEl.addEventListener("click", function (e) {
if (e.target.tagName === "BUTTON" && e.target.dataset.tab) {
tabsEl.querySelectorAll("button").forEach(function (b) {
@@ -263,9 +281,17 @@
}
function escHtml(s) {
- var div = document.createElement("div");
- div.textContent = s == null ? "" : String(s);
- return div.innerHTML;
+ // Plain string replaces rather than the textContent/innerHTML trick: that
+ // trick never escapes quotes, and every caller that builds an attribute
+ // (value="...", title="...", data-help="...") has its value cut short at
+ // the first embedded quote — the rest of the text silently becomes junk
+ // attribute names.
+ return String(s == null ? "" : s)
+ .replace(/&/g, "&")
+ .replace(//g, ">")
+ .replace(/"/g, """)
+ .replace(/'/g, "'");
}
function renderTab(tab) {
diff --git a/web/settings/tabs/control.js b/web/settings/tabs/control.js
index 0e2b95877..9f5056e33 100644
--- a/web/settings/tabs/control.js
+++ b/web/settings/tabs/control.js
@@ -1,4 +1,6 @@
-// Settings → Control tab: site + fuse scalars that feed the PI loop.
+// Settings → Control tab: site + fuse scalars that feed the PI loop,
+// plus the forecast inputs (location, weather source, PV arrays, roof
+// geometry) rendered via tabs/weather.js.
(function () {
var S = (window.FTWSettings = window.FTWSettings || { tabs: {} });
S.tabs = S.tabs || {};
@@ -81,7 +83,16 @@
"Headroom below max amps so the inverter's own per-phase limiter doesn't trip first. Defaults to 0.5 A.",
"0.1") +
'' +
- '';
+ '' +
+ // The forecast inputs — location, weather source, PV arrays, roof
+ // geometry. The sections are owned by tabs/weather.js, which no
+ // longer has a tab button of its own.
+ ((S.tabs.weather && S.tabs.weather.render) ? S.tabs.weather.render(ctx) : "");
+ },
+ after: function (ctx) {
+ // Wire the forecast-input sections rendered above (map, PV arrays,
+ // roof derivation) — weather.js owns their behaviour.
+ if (S.tabs.weather && S.tabs.weather.after) S.tabs.weather.after(ctx);
},
};
})();
diff --git a/web/settings/tabs/planner.js b/web/settings/tabs/planner.js
index 579126aa7..c6cc9bda0 100644
--- a/web/settings/tabs/planner.js
+++ b/web/settings/tabs/planner.js
@@ -1,4 +1,6 @@
-// Settings → Planner tab: MPC planner scalars.
+// Settings → Planner tab: MPC planner scalars. The forecast inputs the
+// planner consumes (location, weather source, PV arrays, roof geometry)
+// render on the Control tab via tabs/weather.js.
(function () {
var S = (window.FTWSettings = window.FTWSettings || { tabs: {} });
S.tabs = S.tabs || {};
@@ -71,7 +73,7 @@
kHtml = 'PV forecast safety k is not set in YAML. The Plan card slider owns it, anywhere from 0 to 2 in steps of 0.05.
' +
'';
}
- return '';
+ var engineHtml = '' +
'Engine controls — leave these unless you are debugging.
' +
'' +
- '' +
+ '';
+ return mpcHtml + engineHtml +
'' +
- 'The planner requires working price + weather forecasts. When disabled the system runs in the manual mode set on the Control page.' +
+ 'The planner reads its forecast inputs — location, weather source, PV arrays, roof geometry — from the Control tab, and needs a working price forecast from the Price tab. When disabled the system runs in the manual mode set on the Control page.' +
'
';
},
after: function (ctx) {
@@ -221,6 +224,7 @@
})
.catch(function () {}); // unreachable → line stays hidden
}
+
},
};
diff --git a/web/settings/tabs/weather.js b/web/settings/tabs/weather.js
index 2eb6b174c..a380ef71c 100644
--- a/web/settings/tabs/weather.js
+++ b/web/settings/tabs/weather.js
@@ -1,31 +1,82 @@
-// Settings → Weather tab: forecast provider + location + PV arrays.
-// Owns its own Leaflet loader + PV-array editor + 3D preview loader
-// so the Settings shell stays weather-agnostic.
+// Forecast-input sections of the Settings → Control tab: forecast
+// provider + location + PV arrays. Registered as S.tabs.weather but
+// rendered inside the Control tab (control.js delegates here); there is
+// no Weather tab button. Owns its own MapLibre loader + PV-array editor
+// + 3D preview loader so the Settings shell stays weather-agnostic.
(function () {
var S = (window.FTWSettings = window.FTWSettings || { tabs: {} });
S.tabs = S.tabs || {};
- var leafletLoading = null;
- function loadLeaflet() {
- if (window.L) return Promise.resolve();
- if (leafletLoading) return leafletLoading;
- leafletLoading = new Promise(function (resolve, reject) {
+ // MapLibre GL JS (BSD-3), vendored under /vendor/maplibre and loaded on
+ // demand exactly the way this tab already lazy-loads its other heavy
+ // optional dependency — the picker is ~1 MB and only these sections need
+ // it. Shipping it on the box follows the same policy as /vendor/three and
+ // the Leaflet copy this replaced: no third-party JS from a CDN, and the map
+ // must load when the gateway cannot reach the internet.
+ //
+ // v6 ships ESM only and is code-split (the entry pulls in
+ // maplibre-gl-shared.mjs and spawns maplibre-gl-worker.mjs by relative
+ // URL), so it is loaded with a dynamic import() rather than a