feat(web-standalone): antenna height from the tallest building within 50 m of a repeater

AGENTS.md

@@ -30,6 +30,7 @@ web-standalone/         Self-contained map tool (no DB)
   changelog.json        "What's new" entries, newest first; every user-visible change adds one
   server.py             stdlib-only HTTP server + proxy for upstream feeds (mc-radar, meshcore.io, PDOK, mc-spamdetector.nl)
   spamdetector.py       parses a spam-detector attack page into Step 1 clues (pure logic, tested)
+  bag3d.py              3D BAG: RD conversion, footprints, the tallest roof within reach of a node (pure logic, tested)
   tests/                pytest for the pure logic; fixtures are trimmed copies of real upstream pages
   tools/                accuracy.mjs (estimator error against known targets), its fixture, build-fixture.py
   Dockerfile            Alpine, non-root, stdlib-only
@@ -179,7 +180,7 @@ new hardcoded hex values.
 
 ### Don't hammer upstream feeds per-node
 
-The proxied upstream APIs (mc-radar, map.meshcore.io, PDOK, mc-spamdetector.nl) are third-party
+The proxied upstream APIs (mc-radar, map.meshcore.io, PDOK, mc-spamdetector.nl, 3D BAG) are third-party
 and rate-limitable; mc-radar answers 429 after a few hundred proven-link fetches in a day.
 Batch/cache lookups (as the existing code does, including the in-memory cache in `server.py`
 for the spam detector) rather than firing one request per node/prefix.

README.md

@@ -58,6 +58,7 @@ control in front if you do; the proxy has no auth of its own.
 | map.meshcore.io | Global node feed (fallback). |
 | PDOK AHN | Terrain/surface elevation (NL). |
 | mc-spamdetector.nl | Incident list, and an incident's entry hops as Step 1 clues. |
+| 3D BAG (api.3dbag.nl) | The tallest building within 50 m of an observer, for its antenna height in the line-of-sight check (NL). |
 
 Loaded by the page directly, not through the proxy: OpenFreeMap vector tiles (the Liberty basemap, no key) and AWS Terrain Tiles (relief in the 3D view, for looking at; line of sight is measured on AHN).
 

web-standalone/.dockerignore

@@ -1,5 +1,6 @@
 *
 !server.py
 !spamdetector.py
+!bag3d.py
 !index.html
 !changelog.json

web-standalone/Dockerfile

@@ -7,7 +7,7 @@ LABEL org.opencontainers.image.source="https://github.com/khagele/meshcore-trian
       org.opencontainers.image.description="MeshCore triangulator (web-standalone): static UI + stdlib proxy to mc-radar/meshcore/PDOK"
 
 WORKDIR /app
-COPY server.py spamdetector.py index.html changelog.json ./
+COPY server.py spamdetector.py bag3d.py index.html changelog.json ./
 
 # Inside the container we must bind all interfaces; restrict exposure when
 # publishing the port instead (e.g. -p 127.0.0.1:8000:8000).

web-standalone/bag3d.py

@@ -0,0 +1,196 @@
+"""3D BAG lookups for the triangulator (#74).
+
+A repeater on a roof has its antenna at roof height plus a mast, not at the
+12 m the estimator assumed for every node. 3D BAG (TU Delft, CC BY 4.0)
+publishes per-building footprints with roof and ground heights. This module
+turns the responses of its OGC API Features endpoint into "the tallest
+building within reach of this point, how high its roof is above the ground,
+and how far away it is". Advertised positions are metres off and a mast can
+stand next to the building it serves, so the tallest building within
+50 m is the antenna's structure rather than the footprint under the point.
+
+Pure logic, stdlib only, so it is unit-testable (tests/test_bag3d.py).
+server.py fetches the payload and serves the result under /proxy/3dbag/.
+
+The API takes bounding boxes in EPSG:7415 (RD New + NAP), so lat/lon is
+converted with the usual polynomial approximation, good to about 25 cm.
+"""
+import math
+
+# Schreutelkops-Adema approximation, WGS84 -> RD New.
+_RD_X0, _RD_Y0 = 155000.0, 463000.0
+_LAT0, _LON0 = 52.15517440, 5.38720621
+
+_RD_X_TERMS = [(0, 1, 190094.945), (1, 1, -11832.228), (2, 1, -114.221), (0, 3, -32.391),
+               (1, 0, -0.705), (3, 1, -2.340), (1, 3, -0.608), (0, 2, -0.008), (2, 3, 0.148)]
+_RD_Y_TERMS = [(1, 0, 309056.544), (0, 2, 3638.893), (2, 0, 73.077), (1, 2, -157.984),
+               (3, 0, 59.788), (0, 1, 0.433), (2, 2, -6.439), (1, 1, -0.032), (0, 4, 0.092),
+               (1, 4, -0.054)]
+_WGS_LAT_TERMS = [(0, 1, 3235.65389), (2, 0, -32.58297), (0, 2, -0.24750), (2, 1, -0.84978),
+                  (0, 3, -0.06550), (2, 2, -0.01709), (1, 0, -0.00738), (4, 0, 0.00530),
+                  (2, 3, -0.00039), (4, 1, 0.00033), (1, 1, -0.00012)]
+_WGS_LON_TERMS = [(1, 0, 5260.52916), (1, 1, 105.94684), (1, 2, 2.45656), (3, 0, -0.81885),
+                  (1, 3, 0.05594), (3, 1, -0.05607), (0, 1, 0.01199), (3, 2, -0.00256),
+                  (1, 4, 0.00128), (0, 2, 0.00022), (2, 0, -0.00022), (5, 0, 0.00026)]
+
+
+def wgs84_to_rd(lat, lon):
+    """RD New x, y in metres for a WGS84 position."""
+    dlat = 0.36 * (lat - _LAT0)
+    dlon = 0.36 * (lon - _LON0)
+    x = _RD_X0 + sum(c * dlat ** p * dlon ** q for p, q, c in _RD_X_TERMS)
+    y = _RD_Y0 + sum(c * dlat ** p * dlon ** q for p, q, c in _RD_Y_TERMS)
+    return x, y
+
+
+def rd_to_wgs84(x, y):
+    """WGS84 lat, lon for an RD New position. Used by the tests to build a
+    lat/lon inside a known footprint; the app only goes the other way."""
+    dx = (x - _RD_X0) * 1e-5
+    dy = (y - _RD_Y0) * 1e-5
+    lat = _LAT0 + sum(c * dx ** p * dy ** q for p, q, c in _WGS_LAT_TERMS) / 3600.0
+    lon = _LON0 + sum(c * dx ** p * dy ** q for p, q, c in _WGS_LON_TERMS) / 3600.0
+    return lat, lon
+
+
+def bbox_around(lat, lon, radius_m):
+    """An RD bbox string for the API, radius_m either side of the point."""
+    x, y = wgs84_to_rd(lat, lon)
+    return "%.0f,%.0f,%.0f,%.0f" % (x - radius_m, y - radius_m, x + radius_m, y + radius_m)
+
+
+def _transform(payload, feature):
+    meta = payload.get("metadata") or {}
+    transform = meta.get("transform") or feature.get("transform") or {}
+    scale = transform.get("scale") or [1, 1, 1]
+    translate = transform.get("translate") or [0, 0, 0]
+    return scale, translate
+
+
+def buildings(payload):
+    """Every Building in a response: id, roof height above ground, the raw
+    heights, and its footprint rings in RD metres (outer ring first per
+    surface, inner rings are holes)."""
+    out = []
+    for feature in payload.get("features") or []:
+        scale, translate = _transform(payload, feature)
+        vertices = feature.get("vertices") or []
+
+        def xy(index):
+            v = vertices[index]
+            return (v[0] * scale[0] + translate[0], v[1] * scale[1] + translate[1])
+
+        for object_id, obj in (feature.get("CityObjects") or {}).items():
+            if obj.get("type") != "Building":
+                continue
+            attrs = obj.get("attributes") or {}
+            roof = attrs.get("b3_h_dak_max")
+            ground = attrs.get("b3_h_maaiveld")
+            surfaces = []
+            for geometry in obj.get("geometry") or []:
+                if str(geometry.get("lod")) != "0":
+                    continue
+                for surface in geometry.get("boundaries") or []:
+                    rings = [[xy(i) for i in ring] for ring in surface if ring]
+                    if rings:
+                        surfaces.append(rings)
+            out.append({
+                "id": object_id,
+                "roof_m": (roof - ground) if isinstance(roof, (int, float)) and isinstance(ground, (int, float)) else None,
+                "h_dak_max": roof,
+                "h_maaiveld": ground,
+                "h_nok": attrs.get("b3_h_nok"),
+                "surfaces": surfaces,
+            })
+    return out
+
+
+def point_in_ring(x, y, ring):
+    """Ray casting; a point on the edge counts as inside."""
+    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):
+            cross = (x2 - x1) * (y - y1) / (y2 - y1) + x1
+            if x <= cross:
+                inside = not inside
+    return inside
+
+
+def point_in_surface(x, y, rings):
+    if not rings or not point_in_ring(x, y, rings[0]):
+        return False
+    return not any(point_in_ring(x, y, hole) for hole in rings[1:])
+
+
+def building_at(payload, lat, lon):
+    """The building whose footprint contains the point, or None."""
+    x, y = wgs84_to_rd(lat, lon)
+    for building in buildings(payload):
+        if any(point_in_surface(x, y, rings) for rings in building["surfaces"]):
+            return building
+    return None
+
+
+def _segment_distance(px, py, ax, ay, bx, by):
+    dx, dy = bx - ax, by - ay
+    if dx == 0 and dy == 0:
+        return math.hypot(px - ax, py - ay)
+    t = max(0.0, min(1.0, ((px - ax) * dx + (py - ay) * dy) / (dx * dx + dy * dy)))
+    return math.hypot(px - (ax + t * dx), py - (ay + t * dy))
+
+
+def footprint_distance(x, y, surfaces):
+    """Metres from the point to the nearest footprint edge, 0 inside."""
+    best = None
+    for rings in surfaces:
+        if point_in_surface(x, y, rings):
+            return 0.0
+        outer = rings[0] if rings else []
+        n = len(outer)
+        for i in range(n):
+            ax, ay = outer[i]
+            bx, by = outer[(i + 1) % n]
+            d = _segment_distance(x, y, ax, ay, bx, by)
+            if best is None or d < best:
+                best = d
+    return best
+
+
+def tallest_near(payloads, lat, lon, radius_m):
+    """The tallest building whose footprint comes within radius_m of the
+    point, over one or more response pages. Ties go to the nearer one."""
+    if isinstance(payloads, dict):
+        payloads = [payloads]
+    x, y = wgs84_to_rd(lat, lon)
+    best = None
+    total = 0
+    for payload in payloads:
+        for building in buildings(payload):
+            total += 1
+            if building["roof_m"] is None or not building["surfaces"]:
+                continue
+            distance = footprint_distance(x, y, building["surfaces"])
+            if distance is None or distance > radius_m:
+                continue
+            key = (round(building["roof_m"], 1), -distance)
+            if best is None or key > best[0]:
+                best = (key, building, distance)
+    return best, total
+
+
+def roof_for_point(payloads, lat, lon, radius_m=50.0):
+    """What the page asks for: the roof height above ground of the tallest
+    building within radius_m, how far its footprint is, or nulls."""
+    best, total = tallest_near(payloads, lat, lon, radius_m)
+    if not best:
+        return {"roof_m": None, "building": None, "distance_m": None, "buildings": total}
+    _, building, distance = best
+    return {
+        "roof_m": round(building["roof_m"], 1),
+        "building": building["id"],
+        "distance_m": round(distance, 1),
+        "buildings": total,
+    }

web-standalone/changelog.json

@@ -1,4 +1,10 @@
 [
+  {
+    "id": "2026-09-07-antenna-height-from-the-roof",
+    "date": "2026-09-07",
+    "title": "A repeater on a roof gets its real antenna height",
+    "body": "The line-of-sight check assumed every observer's antenna at 12 m. When a cluster is locked the tool now looks up, in 3D BAG, the tallest building within 50 m of each observer; those get its roof plus a 3 m mast, the others keep 12 m. The cluster status says how many, and an observer's popup shows its antenna height and how far that roof is. Netherlands only, like the terrain check."
+  },
   {
     "id": "2026-09-07-fresh-visit",
     "date": "2026-09-07",

web-standalone/index.html

@@ -1737,8 +1737,50 @@ <h3>Unknown Node Estimate</h3>
     const panelCluster = document.getElementById("panel-cluster");
     const panelTerrain = document.getElementById("panel-terrain");
 
+    // Assumed when no building stands within ROOF_SEARCH_M of the node. Half
+    // of the fixture's observers stand within 12 m of one (32 of 60 sampled
+    // against 3D BAG), roof median 10.2 m above ground, p90 25 m. Advertised
+    // positions are metres off and a mast can stand next to the building it
+    // serves, so the tallest building within 50 m counts, and its roof plus
+    // a mast replaces this constant (#74).
     const REPEATER_ANTENNA_HEIGHT_M = 12;
+    const ROOFTOP_MAST_M = 3;
+    const ROOF_SEARCH_M = 50;
     const TARGET_ANTENNA_HEIGHT_M = 1.5;
+
+    // Antenna height above the ground for an observer: its roof plus a mast
+    // when 3D BAG put it on a building, the flat constant otherwise.
+    function observerAntennaHeightM(node) {
+      return Number.isFinite(node.roofHeightM) ? node.roofHeightM + ROOFTOP_MAST_M : REPEATER_ANTENNA_HEIGHT_M;
+    }
+
+    // One 3D BAG lookup per node, through server.py, which caches by
+    // position. Sets node.roofHeightM (metres above ground) and
+    // node.roofDistanceM when a building stands within ROOF_SEARCH_M; leaves
+    // them unset otherwise or when the lookup fails. Counts are returned for
+    // the status line.
+    async function fetchRoofHeights(nodes) {
+      let onRoof = 0;
+      let answered = 0;
+      await Promise.all(nodes.map(async (node) => {
+        if (!Number.isFinite(node.lat) || !Number.isFinite(node.lon)) return;
+        try {
+          const response = await fetch(`/proxy/3dbag/roof?lat=${node.lat.toFixed(5)}&lon=${node.lon.toFixed(5)}`);
+          if (!response.ok) return;
+          const data = await response.json();
+          answered += 1;
+          if (Number.isFinite(data.roof_m) && data.roof_m > 0) {
+            node.roofHeightM = data.roof_m;
+            node.roofDistanceM = Number.isFinite(data.distance_m) ? data.distance_m : null;
+            node.roofBuilding = data.building;
+            onRoof += 1;
+          }
+        } catch (error) {
+          // A failed lookup leaves the constant in place; the status line says how many answered.
+        }
+      }));
+      return { onRoof, answered, asked: nodes.length };
+    }
     // Swept against the accuracy fixture (#53, #69): 2 km to 20 km moves median
     // error 1.9 -> 2.0 km and p90 9.0 -> 9.1 km. This constant is inert. Error
     // is set by observer geometry, not by assumed coverage size.
@@ -2905,6 +2947,10 @@ <h3>Unknown Node Estimate</h3>
         if (moved) {
           node.moved = true;
           delete node.provenRadiusKm;
+          // The roof belonged to the position it left.
+          delete node.roofHeightM;
+          delete node.roofDistanceM;
+          delete node.roofBuilding;
         } else {
           delete node.moved;
         }
@@ -3361,7 +3407,7 @@ <h3>Unknown Node Estimate</h3>
         fetchAhnSurfaceElevation(to.lat, to.lon),
         Promise.all(interior.map((point) => fetchAhnLayerValue(point.lat, point.lon, "dsm_05m", surfaceCache)))
       ]);
-      const startHeight = fromGround + REPEATER_ANTENNA_HEIGHT_M;
+      const startHeight = fromGround + observerAntennaHeightM(from);
       const endHeight = targetGround + TARGET_ANTENNA_HEIGHT_M;
       const pathM = distanceKm * 1000;
       let worstLossDb = 0;
@@ -4193,6 +4239,13 @@ <h3>Unknown Node Estimate</h3>
       if (runId !== discoveryRunId) {
         return;
       }
+      // Roof heights for the same nodes the link budget covers (#74). The
+      // cluster copies are what the pools are built from below, so the value
+      // travels with the node from here on.
+      const roofs = await fetchRoofHeights(fetchTargets);
+      if (runId !== discoveryRunId) {
+        return;
+      }
       const supportNodes = Array.from(supportNodeMap.values()).slice(0, 40);
       const observerPoolMap = new Map();
       clusterNodes.forEach((node) => {
@@ -4285,7 +4338,11 @@ <h3>Unknown Node Estimate</h3>
         radiusNotes.push(`${observersWithoutReceiveQuorum} observer${observersWithoutReceiveQuorum === 1 ? "" : "s"} had under ${RECEIVE_QUORUM} inbound links, so used all directions`);
       }
       const implausibleNote = radiusNotes.length ? ` Coverage radii: ${radiusNotes.join("; ")}.` : "";
-      clusterStatus.textContent = `Cluster near ${centroid.lat.toFixed(3)}, ${centroid.lon.toFixed(3)}: ${uniquePrefixCount(dedupedNodes)} observers active.${hullNote}${dedupNote}${linkCapNote}${implausibleNote}`;
+      // Which observers got a roof height, and which stay on the constant.
+      const roofNote = roofs.answered
+        ? ` Antenna heights: ${roofs.onRoof} of ${roofs.answered} observers have a building within ${ROOF_SEARCH_M} m per 3D BAG and use its roof, the rest assume ${REPEATER_ANTENNA_HEIGHT_M} m.`
+        : ` Antenna heights: 3D BAG did not answer, all observers assume ${REPEATER_ANTENNA_HEIGHT_M} m.`;
+      clusterStatus.textContent = `Cluster near ${centroid.lat.toFixed(3)}, ${centroid.lon.toFixed(3)}: ${uniquePrefixCount(dedupedNodes)} observers active.${hullNote}${dedupNote}${linkCapNote}${implausibleNote}${roofNote}`;
       syncClusterEditorFromLockedRegion();
       setDiscoveryCollapsed(true);
       renderCandidateMap([]);
@@ -4677,6 +4734,9 @@ <h3>Unknown Node Estimate</h3>
               lat: n.lat, lon: n.lon,
               weight: Number.isFinite(n.weight) ? n.weight : 1,
               provenRadiusKm: Number.isFinite(n.provenRadiusKm) ? n.provenRadiusKm : null,
+              antennaM: observerAntennaHeightM(n),
+              onRoof: Number.isFinite(n.roofHeightM) ? 1 : 0,
+              roofDistanceM: Number.isFinite(n.roofDistanceM) ? n.roofDistanceM : null,
               moved: n.moved ? 1 : 0,
               movedKm: n.moved && positionOverrides.get(nodeId(n))
                 ? haversineKm(positionOverrides.get(nodeId(n)).upstream, n) : null
@@ -5285,10 +5345,17 @@ <h3>Unknown Node Estimate</h3>
           ? "" : `<br>Proven-link radius: ${Number(p.provenRadiusKm).toFixed(1)} km`;
         const moved = p.moved === 1 || p.moved === "1";
         const movedLine = moved ? `<br>Moved ${Number(p.movedKm).toFixed(1)} km from its advertised position; range from link count.` : "";
+        const onRoof = p.onRoof === 1 || p.onRoof === "1";
+        const roofWhere = onRoof
+          ? (Number(p.roofDistanceM) > 0 ? `tallest roof within ${Number(p.roofDistanceM).toFixed(0)} m, per 3D BAG, + mast` : "its roof per 3D BAG + mast")
+          : "assumed";
+        const antennaLine = Number.isFinite(Number(p.antennaM))
+          ? `<br>Antenna ${Number(p.antennaM).toFixed(0)} m (${roofWhere})`
+          : "";
         const popup = show3dPopup(e.lngLat, `
           <div>
             <p class="popup-title">${p.title}</p>
-            <p class="popup-copy">Short ID: ${p.shortId}<br>Lat/Lon: ${Number(p.lat).toFixed(5)}, ${Number(p.lon).toFixed(5)}<br>Weight: ${String(Number(p.weight).toFixed(2)).replace(/\.00$/, "")}${radius}${movedLine}</p>
+            <p class="popup-copy">Short ID: ${p.shortId}<br>Lat/Lon: ${Number(p.lat).toFixed(5)}, ${Number(p.lon).toFixed(5)}<br>Weight: ${String(Number(p.weight).toFixed(2)).replace(/\.00$/, "")}${radius}${antennaLine}${movedLine}</p>
             ${moved ? `<button type="button" class="link-remove-btn" data-putback="${p.nodeKey}">Put it back</button>` : ""}
           </div>
         `);