feat(canopy): add tree-canopy height layer for path obstruction analysis

Adds Microwaveprop.Canopy module that reads 1° × 1° uint8 per-degree
canopy-height tiles (mirrors SRTM .hgt naming/layout, 30 m resolution).
Returns 0 m for areas with no tile on disk so the lookup is safe to
thread through the existing path-clearance pipeline before any data
lands.

Wires canopy into:
  - PathTerrain.obstacle_top: per-sample blocker = ground +
    max(building_m, canopy_m), so the rover ranks paths through forests
    as obstructed even when no buildings sit on the line
  - CandidateDetail profile: each sample now carries canopy_m alongside
    building_m
  - Rover-detail SVG: green "trees" polygon stacked on terrain, under
    the red building polygon
  - Path calculator elevation chart: green Tree Canopy dataset between
    terrain and buildings

Data prep (download Potapov 2019 / Lang 2022 GeoTIFF, slice to per-degree
uint8 tiles, drop in /data/canopy/) is a separate one-shot operation —
without tiles, lookups return 0 and the existing behavior is unchanged.

assets/js/elevation_profile_hook.ts

@@ -9,6 +9,7 @@ interface ProfilePoint {
   beam: number
   r1: number
   building_m?: number
+  canopy_m?: number
 }
 
 interface DuctRaw {
@@ -201,6 +202,14 @@ export const ElevationProfile: ElevationProfileHook = {
       return (p.elev + (p.building_m || 0) + bulge) * M2FT
     })
 
+    // Canopy: terrain + nearby tree-canopy height (GEDI/Potapov 30m).
+    const hasCanopy = points.some(p => (p.canopy_m || 0) > 1)
+    const canopyTops = points.map((p, i) => {
+      const f = n > 0 ? i / n : 0
+      const bulge = (f * dTotal * (1 - f) * dTotal) / (2 * k * R)
+      return (p.elev + (p.canopy_m || 0) + bulge) * M2FT
+    })
+
     // LOS beam: straight line from antenna A to antenna B
     const losLine = points.map(p => p.beam * M2FT)
     const fresnelLower = points.map(p => (p.beam - p.r1) * M2FT)
@@ -217,7 +226,13 @@ export const ElevationProfile: ElevationProfileHook = {
       likely: d.likely || false
     }))
 
-    const allElevs = [...elevations, ...losLine, ...earthSurface, ...(hasBuildings ? buildingTops : [])]
+    const allElevs = [
+      ...elevations,
+      ...losLine,
+      ...earthSurface,
+      ...(hasBuildings ? buildingTops : []),
+      ...(hasCanopy ? canopyTops : [])
+    ]
     let minY = Math.min(...allElevs) - 40
     let maxY = Math.max(...allElevs) + 120
 
@@ -258,17 +273,30 @@ export const ElevationProfile: ElevationProfileHook = {
             tension: 0.1,
             order: 4
           },
+          ...(hasCanopy ? [
+            {
+              label: "Tree Canopy",
+              data: canopyTops,
+              borderColor: "#16a34a",
+              backgroundColor: "rgba(22, 163, 74, 0.35)",
+              fill: "-1",
+              pointRadius: 0,
+              borderWidth: 1,
+              tension: 0.1,
+              order: 3
+            }
+          ] : []),
           ...(hasBuildings ? [
             {
               label: "Buildings",
               data: buildingTops,
               borderColor: "#dc2626",
               backgroundColor: "rgba(220, 38, 38, 0.45)",
-              fill: "-1",
+              fill: hasCanopy ? "-2" : "-1",
               pointRadius: 0,
               borderWidth: 1,
               tension: 0.1,
-              order: 3
+              order: 2
             }
           ] : []),
           {

lib/microwaveprop/canopy.ex

@@ -0,0 +1,99 @@
+defmodule Microwaveprop.Canopy do
+  @moduledoc """
+  Tree-canopy heights from a 30 m global canopy-height grid (Potapov et al.
+  2019, Lang et al. 2022, or any source able to produce 1° × 1° uint8
+  per-pixel rasters).
+
+  Tile format: `@samples × @samples` bytes, row 0 at the north edge,
+  one byte per pixel = canopy height in meters (0–255). Filename
+  `N{lat}W{lon}.canopy` mirrors the SRTM `.hgt` convention.
+
+  Returns 0 m for points without a tile on disk so the lookup can be
+  threaded through the path-clearance pipeline without raising on
+  un-fetched areas (treat unknown as treeless).
+  """
+
+  require Logger
+
+  @samples 3600
+
+  @spec tile_filename(float(), float()) :: String.t()
+  def tile_filename(lat, lon) do
+    lat_floor = floor(lat)
+    lon_floor = floor(lon)
+
+    lat_prefix = if lat_floor >= 0, do: "N", else: "S"
+    lon_prefix = if lon_floor >= 0, do: "E", else: "W"
+
+    lat_str = lat_floor |> abs() |> Integer.to_string() |> String.pad_leading(2, "0")
+    lon_str = lon_floor |> abs() |> Integer.to_string() |> String.pad_leading(3, "0")
+
+    "#{lat_prefix}#{lat_str}#{lon_prefix}#{lon_str}.canopy"
+  end
+
+  @doc """
+  Returns the canopy height in meters at `lat, lon`, or 0 if no tile is
+  on disk for that 1° square.
+  """
+  @spec lookup(float(), float(), String.t(), keyword()) :: non_neg_integer()
+  def lookup(lat, lon, tiles_dir \\ default_dir(), opts \\ []) do
+    samples = Keyword.get(opts, :samples, @samples)
+    path = Path.join(tiles_dir, tile_filename(lat, lon))
+
+    case :file.open(path, [:read, :binary, :raw]) do
+      {:ok, fd} ->
+        result = read_height(fd, lat, lon, samples)
+        _ = :file.close(fd)
+        result
+
+      {:error, _} ->
+        0
+    end
+  end
+
+  @doc """
+  Batch lookup. Groups points by tile so each tile is opened once.
+  Returns `%{ {lat, lon} => height_m }`. Missing tiles yield 0.
+  """
+  @spec lookup_many([{float(), float()}], String.t(), keyword()) ::
+          %{{float(), float()} => non_neg_integer()}
+  def lookup_many(points, tiles_dir \\ default_dir(), opts \\ []) do
+    samples = Keyword.get(opts, :samples, @samples)
+
+    points
+    |> Enum.group_by(fn {lat, lon} -> {floor(lat), floor(lon)} end)
+    |> Enum.flat_map(fn {_tile_key, pts} -> read_tile_points(pts, tiles_dir, samples) end)
+    |> Map.new()
+  end
+
+  defp read_tile_points([{lat0, lon0} | _] = pts, tiles_dir, samples) do
+    path = Path.join(tiles_dir, tile_filename(lat0, lon0))
+
+    case :file.open(path, [:read, :binary, :raw]) do
+      {:ok, fd} ->
+        result =
+          Enum.map(pts, fn {lat, lon} = pt -> {pt, read_height(fd, lat, lon, samples)} end)
+
+        _ = :file.close(fd)
+        result
+
+      {:error, _} ->
+        Enum.map(pts, fn pt -> {pt, 0} end)
+    end
+  end
+
+  defp read_height(fd, lat, lon, samples) do
+    row = round((floor(lat) + 1 - lat) * (samples - 1))
+    col = round((lon - floor(lon)) * (samples - 1))
+    offset = row * samples + col
+
+    case :file.pread(fd, offset, 1) do
+      {:ok, <<h::unsigned-integer-size(8)>>} -> h
+      _ -> 0
+    end
+  end
+
+  defp default_dir do
+    Application.get_env(:microwaveprop, :canopy_tiles_dir, "/data/canopy")
+  end
+end

lib/microwaveprop/rover/candidate_detail.ex

@@ -6,6 +6,7 @@ defmodule Microwaveprop.Rover.CandidateDetail do
   """
 
   alias Microwaveprop.Buildings.Index, as: BuildingsIndex
+  alias Microwaveprop.Canopy
   alias Microwaveprop.Radio.Maidenhead
   alias Microwaveprop.Rover.DriveTime
   alias Microwaveprop.Rover.LinkMargin
@@ -30,7 +31,7 @@ defmodule Microwaveprop.Rover.CandidateDetail do
           station_elev_m: integer() | nil,
           max_obstacle_m: integer() | nil,
           clearance_m: integer() | nil,
-          profile: [%{dist_km: float(), elev: integer(), building_m: float()}]
+          profile: [%{dist_km: float(), elev: integer(), building_m: float(), canopy_m: integer()}]
         }
 
   @spec summarize(map(), [station()], non_neg_integer(), DateTime.t(), atom()) :: %{
@@ -49,7 +50,7 @@ defmodule Microwaveprop.Rover.CandidateDetail do
         max_obstacle =
           profile
           |> middle_samples()
-          |> Enum.map(fn pt -> pt.elev + round(pt.building_m) end)
+          |> Enum.map(fn pt -> pt.elev + round(max(pt.building_m, pt.canopy_m)) end)
           |> safe_max()
 
         clearance =
@@ -103,7 +104,8 @@ defmodule Microwaveprop.Rover.CandidateDetail do
           %{
             dist_km: pt.dist_km,
             elev: pt.elev,
-            building_m: BuildingsIndex.max_height_near(pt.lat, pt.lon, @profile_building_radius_m)
+            building_m: BuildingsIndex.max_height_near(pt.lat, pt.lon, @profile_building_radius_m),
+            canopy_m: Canopy.lookup(pt.lat, pt.lon)
           }
         end)
 

lib/microwaveprop/rover/path_terrain.ex

@@ -10,6 +10,7 @@ defmodule Microwaveprop.Rover.PathTerrain do
   """
 
   alias Microwaveprop.Buildings.Index, as: BuildingsIndex
+  alias Microwaveprop.Canopy
   alias Microwaveprop.Rover.Elevation
 
   # Number of intermediate samples between rover and station.
@@ -33,6 +34,7 @@ defmodule Microwaveprop.Rover.PathTerrain do
   def clearance_map(cells, stations, opts \\ []) when is_list(cells) and is_list(stations) do
     elev_lookup = Keyword.get(opts, :elev_lookup, &Elevation.lookup_many/1)
     buildings_lookup = Keyword.get(opts, :buildings_lookup, &default_building_height/1)
+    canopy_lookup = Keyword.get(opts, :canopy_lookup, &default_canopy_height/1)
 
     pairs =
       for cell <- cells, station <- stations do
@@ -48,11 +50,11 @@ defmodule Microwaveprop.Rover.PathTerrain do
     elev = elev_lookup.(rover_points ++ sample_points)
 
     Map.new(pairs, fn {rover_pt, station_pt} ->
-      {{rover_pt, station_pt}, clearance(rover_pt, station_pt, elev, buildings_lookup)}
+      {{rover_pt, station_pt}, clearance(rover_pt, station_pt, elev, buildings_lookup, canopy_lookup)}
     end)
   end
 
-  defp clearance(rover_pt, station_pt, elev, buildings_lookup) do
+  defp clearance(rover_pt, station_pt, elev, buildings_lookup, canopy_lookup) do
     rover_elev = Map.get(elev, rover_pt)
 
     case rover_elev do
@@ -63,7 +65,7 @@ defmodule Microwaveprop.Rover.PathTerrain do
         path_max =
           rover_pt
           |> intermediate_samples(station_pt)
-          |> Enum.map(&obstacle_top(&1, elev, buildings_lookup))
+          |> Enum.map(&obstacle_top(&1, elev, buildings_lookup, canopy_lookup))
           |> Enum.reject(&is_nil/1)
           |> safe_max()
 
@@ -71,13 +73,15 @@ defmodule Microwaveprop.Rover.PathTerrain do
     end
   end
 
-  defp obstacle_top({lat, lon} = pt, elev, buildings_lookup) do
+  defp obstacle_top({lat, lon} = pt, elev, buildings_lookup, canopy_lookup) do
     case Map.get(elev, pt) do
       nil ->
         nil
 
       ground ->
-        ground + round(buildings_lookup.({lat, lon}))
+        building_m = buildings_lookup.({lat, lon})
+        canopy_m = canopy_lookup.({lat, lon})
+        ground + round(max(building_m, canopy_m))
     end
   end
 
@@ -85,6 +89,10 @@ defmodule Microwaveprop.Rover.PathTerrain do
     BuildingsIndex.max_height_near(lat, lon, @building_search_radius_m)
   end
 
+  defp default_canopy_height({lat, lon}) do
+    Canopy.lookup(lat, lon)
+  end
+
   defp intermediate_samples({lat1, lon1}, {lat2, lon2}) do
     # Linear interpolation in lat/lon — fine at the distances we work with
     # (≤200 mi). Skip endpoints (i=0 is the rover, i=N+1 is the station).

lib/microwaveprop_web/live/path_live.ex

@@ -6,6 +6,7 @@ defmodule MicrowavepropWeb.PathLive do
 
   alias Microwaveprop.Buildings.Index, as: BuildingsIndex
   alias Microwaveprop.Buildings.Loader, as: BuildingsLoader
+  alias Microwaveprop.Canopy
   alias Microwaveprop.Ionosphere
   alias Microwaveprop.Propagation
   alias Microwaveprop.Propagation.BandConfig
@@ -1639,7 +1640,8 @@ defmodule MicrowavepropWeb.PathLive do
           elev: p.elev,
           beam: p[:beam] || 0,
           r1: p[:r1] || 0,
-          building_m: BuildingsIndex.max_height_near(p.lat, p.lon, 80)
+          building_m: BuildingsIndex.max_height_near(p.lat, p.lon, 80),
+          canopy_m: Canopy.lookup(p.lat, p.lon)
         }
       end)
 

lib/microwaveprop_web/live/rover_live.ex

@@ -804,7 +804,11 @@ defmodule MicrowavepropWeb.RoverLive do
   defp profile_to_svg([]), do: nil
 
   defp profile_to_svg(points) do
-    tops = Enum.map(points, fn pt -> pt.elev + round(Map.get(pt, :building_m, 0.0)) end)
+    tops =
+      Enum.map(points, fn pt ->
+        pt.elev + round(max(Map.get(pt, :building_m, 0.0), Map.get(pt, :canopy_m, 0)))
+      end)
+
     elevs = Enum.map(points, & &1.elev)
     min_e = Enum.min(elevs)
     max_e = max(Enum.max(elevs), Enum.max(tops))
@@ -823,29 +827,31 @@ defmodule MicrowavepropWeb.RoverLive do
     fill = "0,60 " <> line <> " 240,60"
 
     has_buildings? = Enum.any?(points, fn pt -> Map.get(pt, :building_m, 0.0) > 0.5 end)
+    has_canopy? = Enum.any?(points, fn pt -> Map.get(pt, :canopy_m, 0) > 1 end)
 
-    building_fill =
+    terrain_back =
-      if has_buildings? do
+      points
-        building_line =
+      |> Enum.reverse()
-          Enum.map_join(points, " ", fn pt ->
+      |> Enum.map_join(" ", fn pt ->
-            x = pt.dist_km / total_km * 240.0
+        x = pt.dist_km / total_km * 240.0
-            top = pt.elev + round(Map.get(pt, :building_m, 0.0))
+        y = 60.0 - (pt.elev - min_e) / span * 50.0
-            y = 60.0 - (top - min_e) / span * 50.0
+        "#{Float.round(x, 2)},#{Float.round(y, 2)}"
-            "#{Float.round(x, 2)},#{Float.round(y, 2)}"
+      end)
-          end)
 
-        # Close back along the terrain so the building polygon sits on top.
+    obstacle_fill = fn obstacle_key ->
-        terrain_back =
+      obstacle_line =
-          points
+        Enum.map_join(points, " ", fn pt ->
-          |> Enum.reverse()
+          x = pt.dist_km / total_km * 240.0
-          |> Enum.map_join(" ", fn pt ->
+          top = pt.elev + round(Map.get(pt, obstacle_key, 0))
-            x = pt.dist_km / total_km * 240.0
+          y = 60.0 - (top - min_e) / span * 50.0
-            y = 60.0 - (pt.elev - min_e) / span * 50.0
+          "#{Float.round(x, 2)},#{Float.round(y, 2)}"
-            "#{Float.round(x, 2)},#{Float.round(y, 2)}"
+        end)
-          end)
 
-        building_line <> " " <> terrain_back
+      obstacle_line <> " " <> terrain_back
-      end
+    end
+
+    canopy_fill = if has_canopy?, do: obstacle_fill.(:canopy_m)
+    building_fill = if has_buildings?, do: obstacle_fill.(:building_m)
 
     los_y_start = Float.round(60.0 - (first.elev - min_e) / span * 50.0, 2)
     los_y_end = Float.round(60.0 - (last.elev - min_e) / span * 50.0, 2)
@@ -853,6 +859,7 @@ defmodule MicrowavepropWeb.RoverLive do
     %{
       line: line,
       fill: fill,
+      canopy_fill: canopy_fill,
       building_fill: building_fill,
       start_elev_m: first.elev,
       end_elev_m: last.elev,
@@ -1256,6 +1263,13 @@ defmodule MicrowavepropWeb.RoverLive do
     ~H"""
     <svg viewBox="0 0 240 64" class="w-full h-12 mt-1" preserveAspectRatio="none">
       <polygon points={@svg.fill} fill="rgb(14 165 233 / 0.25)" />
+      <polygon
+        :if={@svg.canopy_fill}
+        points={@svg.canopy_fill}
+        fill="rgb(34 197 94 / 0.45)"
+        stroke="#16a34a"
+        stroke-width="0.4"
+      />
       <polygon
         :if={@svg.building_fill}
         points={@svg.building_fill}

test/microwaveprop/canopy_test.exs

@@ -0,0 +1,66 @@
+defmodule Microwaveprop.CanopyTest do
+  use ExUnit.Case, async: true
+
+  alias Microwaveprop.Canopy
+
+  @samples 4
+
+  setup do
+    tiles_dir = Path.join(System.tmp_dir!(), "canopy_test_#{:erlang.unique_integer([:positive])}")
+    File.mkdir_p!(tiles_dir)
+    on_exit(fn -> File.rm_rf!(tiles_dir) end)
+    {:ok, tiles_dir: tiles_dir}
+  end
+
+  defp write_tile(tiles_dir, lat, lon, heights) do
+    path = Path.join(tiles_dir, Canopy.tile_filename(lat, lon))
+    File.write!(path, IO.iodata_to_binary(heights))
+    path
+  end
+
+  test "tile_filename/2 follows SRTM naming convention" do
+    assert Canopy.tile_filename(32.5, -97.5) == "N32W098.canopy"
+    assert Canopy.tile_filename(33.0, -96.0) == "N33W096.canopy"
+  end
+
+  test "lookup/3 returns 0 when tile is missing", %{tiles_dir: tiles_dir} do
+    assert Canopy.lookup(32.5, -97.5, tiles_dir, samples: @samples) == 0
+  end
+
+  test "lookup/3 returns the height at the matching pixel", %{tiles_dir: tiles_dir} do
+    # 4x4 grid, row 0 = north edge, so:
+    # row 0 = lat 33 (N edge)
+    # row 3 = lat 32 (S edge)
+    heights =
+      List.flatten([
+        [10, 11, 12, 13],
+        [20, 21, 22, 23],
+        [30, 31, 32, 33],
+        [40, 41, 42, 43]
+      ])
+
+    write_tile(tiles_dir, 32.5, -97.5, heights)
+
+    # Interior near NW (row 0, col 0) -> 10
+    assert Canopy.lookup(32.999, -97.999, tiles_dir, samples: @samples) == 10
+    # Interior near NE (row 0, col 3) -> 13
+    assert Canopy.lookup(32.999, -97.001, tiles_dir, samples: @samples) == 13
+    # Interior near SW (row 3, col 0) -> 40
+    assert Canopy.lookup(32.001, -97.999, tiles_dir, samples: @samples) == 40
+    # Interior near SE (row 3, col 3) -> 43
+    assert Canopy.lookup(32.001, -97.001, tiles_dir, samples: @samples) == 43
+  end
+
+  test "lookup_many/2 batches lookups grouping by tile", %{tiles_dir: tiles_dir} do
+    heights = List.duplicate(25, @samples * @samples)
+    write_tile(tiles_dir, 32.5, -97.5, heights)
+
+    points = [{32.5, -97.5}, {32.1, -97.9}]
+    result = Canopy.lookup_many(points, tiles_dir, samples: @samples)
+
+    assert result[{32.5, -97.5}] == 25
+    assert result[{32.1, -97.9}] == 25
+    # Point in a tile that doesn't exist is 0.
+    assert Canopy.lookup_many([{40.0, -90.0}], tiles_dir, samples: @samples)[{40.0, -90.0}] == 0
+  end
+end