feat(sim): add land-cover clutter inputs

DISCRETE_EVENT_SIM.md

@@ -9,7 +9,7 @@ To start one simulation with the default configurations, run:
 
 ```python3 loraMesh.py [nr_nodes]```
 
-If no argument is given, you first have to place the nodes on a plot. After you place a node, you can change its [role](https://meshtastic.org/docs/settings/config/device#role), hopLimit, height (elevation) and antenna gain. These settings will automatically save when you place a new node or when you start the simulation.
+If no argument is given, you first have to place the nodes on a plot. After you place a node, you can change its [role](https://meshtastic.org/docs/settings/config/device#role), hopLimit, antenna height above local ground, and antenna gain. These settings will automatically save when you place a new node or when you start the simulation.
 
 ![](/img/configNode.png)
 
@@ -23,6 +23,71 @@ Short deterministic smoke runs can also override the configured duration and mes
 
 ```python3 loraMesh.py 2 --no-gui --simtime-seconds 5 --period-seconds 0.5```
 
+The same headless path can import positioned real-mesh nodes. `--from-map`
+reads a Meshtastic map `/api/v1/nodes` JSON endpoint; the public default is
+`https://meshtastic.liamcottle.net/api/v1/nodes`, but you can pass another
+compatible endpoint URL. These map endpoints usually return a broad node list,
+so pass a local area-of-interest bounding box. `--map-bbox` uses the common
+`min_lat,min_lon,max_lat,max_lon` order that most GIS tools call
+`south,west,north,east`; you can copy those four numbers from OpenStreetMap's
+Export panel, geojson.io's bbox readout, QGIS, or any other tool that shows the
+extent of the map view or selected polygon. Keep the box tight enough for the
+local scenario you want to simulate:
+
+```python3 loraMesh.py --from-map https://meshtastic.liamcottle.net/api/v1/nodes --map-bbox 41.50,41.50,41.82,41.86 --map-limit 50 --no-gui```
+
+You can also import positioned nodes from the NodeDB cached by a local
+Meshtastic device. This uses the Python client `interface.nodesByNum` data that
+backs `meshtastic --nodes`, not the pretty-printed table. Use TCP for a network
+device, or omit `--nodedb-host` to use Meshtastic serial auto-detection. For a
+quick local-device run, pass the device address and cap the imported node count:
+
+```python3 loraMesh.py --from-nodedb --nodedb-host 192.168.1.23 --map-limit 50 --no-gui```
+
+NodeDB often contains old or far-away positions. Add `--map-bbox` when you want
+to restrict the run to one local area:
+
+```python3 loraMesh.py --from-nodedb --nodedb-host 192.168.1.23 --map-bbox 41.50,41.50,41.82,41.86 --map-limit 50 --no-gui```
+
+Imported nodes use the same `HM` antenna height and `hopLimit` defaults as
+generated and file-backed scenarios. Change those config values when the
+position source does not carry the simulation value you want.
+
+Terrain obstruction can be added to map, NodeDB, or origin-backed scenario inputs
+without creating a custom terrain file. `--terrain-srtm` downloads missing SRTM
+HGT tiles from Mapzen Terrain Tiles on AWS into a local cache and feeds the
+terrain grid directly into terrain-aware node geometry:
+
+```python3 loraMesh.py --from-nodedb --nodedb-host 192.168.1.23 --map-limit 50 --terrain-srtm --no-gui```
+
+With an explicit `--map-bbox`, SRTM samples that whole requested rectangle. When
+the terrain bbox is derived from imported or file-backed nodes, Meshtasticator
+keeps the download smaller: it loads tiles around the selected nodes and along
+flat-link candidate paths, instead of downloading every tile in a large
+edge-to-edge rectangle. When publishing screenshots, reports, or derived
+datasets from this terrain source, attribute the terrain data to
+[Mapzen Terrain Tiles on AWS](https://registry.opendata.aws/terrain-tiles/),
+SRTM/NASA, and their underlying open elevation sources:
+
+```python3 loraMesh.py --from-map https://meshtastic.liamcottle.net/api/v1/nodes --map-bbox 41.50,41.50,41.82,41.86 --map-limit 50 --terrain-srtm --no-gui```
+
+Map payload `altitude` values are absolute GPS/MSL altitude, not antenna height,
+so map import keeps using `HM` as the fallback antenna height above local
+ground. When `--terrain-srtm` is enabled, each map node is checked
+against its own SRTM ground sample: plausible positive map altitudes are used as
+absolute node altitude, while missing, below-ground, or implausibly high values
+fall back to `SRTM ground + antenna height` for 3D distance calculations.
+
+Land-cover clutter is a separate optional CSV grid. Use it for broad urban,
+open, water, or forest excess-loss inputs without pretending Meshtasticator is a
+building-level ray tracer:
+
+```python3 loraMesh.py --from-file nodeConfig.yaml --terrain-srtm --clutter-grid clutter.csv --no-gui```
+
+`tools/osm_to_clutter_csv.py` can build a coarse clutter grid from public
+OpenStreetMap building, landuse, natural, and water polygons. The simulator
+never fetches OpenStreetMap data implicitly.
+
 If you placed the nodes yourself, after a simulation the number of nodes, their coordinates and configuration are automatically saved and you can rerun the scenario with:
 
  ```python3 loraMesh.py --from-file```

batchSim.py

@@ -232,7 +232,7 @@ def __init__(self, conf, x, y):
                 x, y = coords[nodeId]
 
                 # We create a NodeConfig object so that MeshNode will use that
-                nodeConfig = NodeConfig(nodeId, Point(x, y, routerTypeConf.HM), routerTypeConf.PERIOD, antenna_gain=routerTypeConf.GL, hop_limit=routerTypeConf.hopLimit)
+                nodeConfig = NodeConfig(nodeId, Point(x, y, routerTypeConf.HM), routerTypeConf.PERIOD, routerTypeConf.PTX, routerTypeConf.FREQ, antenna_gain=routerTypeConf.GL, hop_limit=routerTypeConf.hopLimit)
                 node_configs.append(nodeConfig)
 
             if SHOW_GRAPH:

lib/clutter.py

@@ -0,0 +1,252 @@
+"""Optional land-cover clutter loss for radio links.
+
+Terrain handles hills, curvature, and Fresnel obstruction. It does not know
+whether a lowland path crosses apartment blocks, a beach/coastal opening, or a
+mountain-side vantage point looking down into the city. This module adds that
+separate, data-driven clutter term from a small raster CSV.
+"""
+
+import bisect
+import csv
+import math
+from pathlib import Path
+
+from lib.csv_validation import finite_float, finite_lat_lon
+from lib.terrain import latlon_to_xy, terrain_ground_elevation
+
+
+class ClutterGrid:
+    """Nearest-cell lookup for small land-cover rasters."""
+
+    def __init__(self, samples):
+        self.samples = samples
+        self.xs = sorted({x for x, _, _ in samples})
+        self.ys = sorted({y for _, y, _ in samples})
+        self.by_xy = {(x, y): clutter_class for x, y, clutter_class in samples}
+        self.is_regular = len(self.xs) * len(self.ys) == len(samples)
+
+    @classmethod
+    def from_csv(cls, path, origin_lat=None, origin_lon=None):
+        samples = []
+        with open(path, newline="", encoding="utf-8") as fh:
+            reader = csv.DictReader(fh)
+            for row_number, row in enumerate(reader, start=2):
+                has_lat_lon = "lat" in row and "lon" in row
+                if has_lat_lon and origin_lat is not None and origin_lon is not None:
+                    lat, lon = finite_lat_lon(row, "clutter", row_number)
+                    x, y = latlon_to_xy(lat, lon, origin_lat, origin_lon)
+                elif "x_m" in row and "y_m" in row:
+                    x = finite_float(row, "x_m", "clutter", row_number)
+                    y = finite_float(row, "y_m", "clutter", row_number)
+                elif has_lat_lon:
+                    raise ValueError("lat/lon clutter CSV requires GEO_ORIGIN_LAT and GEO_ORIGIN_LON")
+                else:
+                    raise ValueError("clutter CSV needs x_m/y_m or lat/lon columns")
+
+                clutter_class = row.get("clutter_class")
+                if clutter_class is None or not clutter_class.strip():
+                    raise ValueError("clutter CSV needs clutter_class column")
+                samples.append((x, y, clutter_class.strip().lower()))
+
+        if not samples:
+            raise ValueError(f"clutter CSV has no samples: {path}")
+        return cls(samples)
+
+    @staticmethod
+    def _nearest_axis_value(values, value):
+        index = bisect.bisect_left(values, value)
+        if index <= 0:
+            return values[0]
+        if index >= len(values):
+            return values[-1]
+
+        before = values[index - 1]
+        after = values[index]
+        return before if abs(value - before) <= abs(after - value) else after
+
+    def class_at(self, x, y):
+        if self.is_regular:
+            nearest_x = self._nearest_axis_value(self.xs, x)
+            nearest_y = self._nearest_axis_value(self.ys, y)
+            return self.by_xy[(nearest_x, nearest_y)]
+
+        _, _, clutter_class = min(
+            self.samples,
+            key=lambda sample: math.hypot(x - sample[0], y - sample[1]),
+        )
+        return clutter_class
+
+
+def _clutter_grid(conf):
+    if not conf.CLUTTER_ENABLED or not conf.CLUTTER_GRID_FILE:
+        return None
+
+    # Lat/lon CSVs are projected into scenario-local meters, so the projection
+    # origin is part of the loaded grid identity, not just metadata.
+    cache_identity = (conf.CLUTTER_GRID_FILE, conf.GEO_ORIGIN_LAT, conf.GEO_ORIGIN_LON)
+    cached_identity = getattr(conf, "_clutter_grid_identity", None)
+    if getattr(conf, "_clutter_grid", None) is not None and cached_identity == cache_identity:
+        return conf._clutter_grid
+
+    path = Path(conf.CLUTTER_GRID_FILE)
+    conf._clutter_grid = ClutterGrid.from_csv(path, conf.GEO_ORIGIN_LAT, conf.GEO_ORIGIN_LON)
+    conf._clutter_grid_identity = cache_identity
+    return conf._clutter_grid
+
+
+def _class_loss_db_per_km(conf, clutter_class):
+    if clutter_class in {"urban", "building"}:
+        return conf.CLUTTER_URBAN_LOSS_DB_PER_KM
+    if clutter_class in {"suburban", "residential"}:
+        return conf.CLUTTER_SUBURBAN_LOSS_DB_PER_KM
+    if clutter_class in {"forest", "wood"}:
+        return conf.CLUTTER_FOREST_LOSS_DB_PER_KM
+    if clutter_class in {"water", "coastal_water"}:
+        return conf.CLUTTER_WATER_LOSS_DB_PER_KM
+    return conf.CLUTTER_OPEN_LOSS_DB_PER_KM
+
+
+def clutter_path_features(conf, tx_point, rx_point):
+    """Return coarse land-cover fractions along a radio path.
+
+    The radio calibration model needs reusable features, not node-pair lookup
+    tables. These fractions let a fitted model learn patterns such as "urban
+    lowland paths behave differently from open coastal paths" and apply that
+    lesson to new generated node pairs.
+    """
+    grid = _clutter_grid(conf)
+    if grid is None:
+        return {
+            "urban_fraction": 0.0,
+            "open_fraction": 0.0,
+            "water_fraction": 0.0,
+            "forest_fraction": 0.0,
+            "endpoint_urban_count": 0.0,
+        }
+
+    cache = getattr(conf, "_clutter_feature_cache", None)
+    if cache is None:
+        cache = {}
+        conf._clutter_feature_cache = cache
+
+    cache_key = (
+        round(tx_point.x, 2),
+        round(tx_point.y, 2),
+        round(tx_point.z, 2),
+        round(rx_point.x, 2),
+        round(rx_point.y, 2),
+        round(rx_point.z, 2),
+        conf.CLUTTER_GRID_FILE,
+        conf.GEO_ORIGIN_LAT,
+        conf.GEO_ORIGIN_LON,
+        conf.CLUTTER_PROFILE_SAMPLES,
+    )
+    if cache_key in cache:
+        return cache[cache_key]
+
+    samples = max(1, conf.CLUTTER_PROFILE_SAMPLES)
+    class_counts = {}
+    for index in range(samples):
+        fraction = (index + 0.5) / samples
+        x = tx_point.x + (rx_point.x - tx_point.x) * fraction
+        y = tx_point.y + (rx_point.y - tx_point.y) * fraction
+        clutter_class = grid.class_at(x, y)
+        class_counts[clutter_class] = class_counts.get(clutter_class, 0) + 1
+
+    endpoint_urban_count = 0
+    for point in (tx_point, rx_point):
+        endpoint_class = grid.class_at(point.x, point.y)
+        if endpoint_class in {"urban", "building", "suburban", "residential"}:
+            endpoint_urban_count += 1
+
+    features = {
+        "urban_fraction": class_counts.get("urban", 0) / samples,
+        "open_fraction": class_counts.get("open", 0) / samples,
+        "water_fraction": (class_counts.get("water", 0) + class_counts.get("coastal_water", 0)) / samples,
+        "forest_fraction": (class_counts.get("forest", 0) + class_counts.get("wood", 0)) / samples,
+        "endpoint_urban_count": float(endpoint_urban_count),
+    }
+    cache[cache_key] = features
+    return features
+
+
+def _is_high_vantage(conf, point):
+    ground = terrain_ground_elevation(conf, point)
+    return ground is not None and ground >= conf.CLUTTER_HIGH_VANTAGE_ELEVATION_M
+
+
+def clutter_obstruction_loss(conf, tx_point, rx_point):
+    """Estimate extra land-cover clutter loss in dB for a TX/RX path."""
+    grid = _clutter_grid(conf)
+    if grid is None:
+        return 0.0
+
+    cache = getattr(conf, "_clutter_loss_cache", None)
+    if cache is None:
+        cache = {}
+        conf._clutter_loss_cache = cache
+
+    cache_key = (
+        round(tx_point.x, 2),
+        round(tx_point.y, 2),
+        round(tx_point.z, 2),
+        round(rx_point.x, 2),
+        round(rx_point.y, 2),
+        round(rx_point.z, 2),
+        conf.CLUTTER_GRID_FILE,
+        conf.GEO_ORIGIN_LAT,
+        conf.GEO_ORIGIN_LON,
+        conf.CLUTTER_PROFILE_SAMPLES,
+        conf.CLUTTER_URBAN_LOSS_DB_PER_KM,
+        conf.CLUTTER_SUBURBAN_LOSS_DB_PER_KM,
+        conf.CLUTTER_FOREST_LOSS_DB_PER_KM,
+        conf.CLUTTER_OPEN_LOSS_DB_PER_KM,
+        conf.CLUTTER_WATER_LOSS_DB_PER_KM,
+        conf.CLUTTER_URBAN_ENDPOINT_LOSS_DB,
+        conf.CLUTTER_HIGH_VANTAGE_ELEVATION_M,
+        conf.CLUTTER_HIGH_VANTAGE_LOSS_FACTOR,
+        conf.CLUTTER_COASTAL_PATH_LOSS_FACTOR,
+        conf.CLUTTER_COASTAL_SAMPLE_FRACTION,
+        conf.CLUTTER_MAX_LOSS_DB,
+    )
+    if cache_key in cache:
+        return cache[cache_key]
+
+    horizontal_distance = math.hypot(rx_point.x - tx_point.x, rx_point.y - tx_point.y)
+    if horizontal_distance <= 0:
+        return 0.0
+
+    samples = max(1, conf.CLUTTER_PROFILE_SAMPLES)
+    class_counts = {}
+    path_loss_rate = 0.0
+    for index in range(samples):
+        fraction = (index + 0.5) / samples
+        x = tx_point.x + (rx_point.x - tx_point.x) * fraction
+        y = tx_point.y + (rx_point.y - tx_point.y) * fraction
+        clutter_class = grid.class_at(x, y)
+        class_counts[clutter_class] = class_counts.get(clutter_class, 0) + 1
+        path_loss_rate += _class_loss_db_per_km(conf, clutter_class)
+
+    path_loss = (path_loss_rate / samples) * (horizontal_distance / 1000.0)
+
+    # Coastal or sea-adjacent paths are often real line-of-sight corridors. Do
+    # not let a few nearby urban cells make them look like street-canyon links.
+    water_samples = class_counts.get("water", 0) + class_counts.get("coastal_water", 0)
+    open_samples = class_counts.get("open", 0) + class_counts.get("beach", 0)
+    if (water_samples + open_samples) / samples >= conf.CLUTTER_COASTAL_SAMPLE_FRACTION:
+        path_loss *= conf.CLUTTER_COASTAL_PATH_LOSS_FACTOR
+
+    tx_high = _is_high_vantage(conf, tx_point)
+    rx_high = _is_high_vantage(conf, rx_point)
+    if tx_high or rx_high:
+        path_loss *= conf.CLUTTER_HIGH_VANTAGE_LOSS_FACTOR
+
+    endpoint_loss = 0.0
+    for point, high_vantage in ((tx_point, tx_high), (rx_point, rx_high)):
+        endpoint_class = grid.class_at(point.x, point.y)
+        if endpoint_class in {"urban", "building", "suburban", "residential"} and not high_vantage:
+            endpoint_loss += conf.CLUTTER_URBAN_ENDPOINT_LOSS_DB
+
+    loss = min(path_loss + endpoint_loss, conf.CLUTTER_MAX_LOSS_DB)
+    cache[cache_key] = loss
+    return loss

lib/common.py

@@ -1,11 +1,15 @@
 import random
-import os
 
 import numpy as np
 
 from lib import phy
 from lib.point import Point
 
+
+def node_antenna_height(node):
+    """Return antenna height above ground, falling back to legacy Point.z."""
+    return getattr(node, "antennaHeight", getattr(node, "antenna_height", node.position.z))
+
 def find_random_position(conf, node_configs) -> (float, float):
     """Given a simulation config and list of existing node configs/nodes, find a
     randomly chosen position for the next node such that it is within the
@@ -79,7 +83,7 @@ def setup_asymmetric_links(conf, nodes):
                 nodeA = nodes[a]
                 nodeB = nodes[b]
                 distAB = nodeA.position.euclidean_distance(nodeB.position)
-                pathLossAB = phy.estimate_path_loss(conf, distAB, conf.FREQ, nodeA.position.z, nodeB.position.z)
+                pathLossAB = phy.estimate_path_loss(conf, distAB, conf.FREQ, node_antenna_height(nodeA), node_antenna_height(nodeB))
 
                 offsetAB = conf.LINK_OFFSET[(a, b)]
                 offsetBA = conf.LINK_OFFSET[(b, a)]