feat(rover): add canopy clutter penalty + Canopy.BulkFetch downloader

Rover ranking now applies a per-cell tree-canopy clutter penalty (1 dB / 8 m, capped at 4 dB) so cells in dense forests get down-ranked even when their per-station path clearance survives. PathTerrain already accounts for trees ON the link path; this surfaces "I'm in a forest, every direction is foliage." Canopy.BulkFetch downloads Lang et al. 2020 10 m global canopy-height 3-degree COGs and slices each into 1° × 1° uint8 .canopy tiles via gdal_translate, the runtime image now ships gdal-bin. Run on prod via: bin/microwaveprop rpc 'Microwaveprop.Canopy.BulkFetch.run(32.8, -97.0, 200)' Tile cache lives under /data/canopy/. Already-sliced tiles are skipped.
2026-04-26 12:57:49 -05:00 · 2026-04-26 12:57:49 -05:00 · 7105ac61df
commit 7105ac61df
parent 1a85a5f540
3 changed files with 222 additions and 5 deletions
--- a/2
+++ b/2
@ -171,7 +171,7 @@ RUN --mount=type=cache,target=/var/cache/apt,sharing=locked \
    --mount=type=cache,target=/var/lib/apt,sharing=locked \
    rm -f /etc/apt/apt.conf.d/docker-clean \
    && apt-get update \
-    && apt-get install -y --no-install-recommends cdo
+    && apt-get install -y --no-install-recommends cdo gdal-bin

 RUN sed -i '/en_US.UTF-8/s/^# //g' /etc/locale.gen && locale-gen

--- a/lib/microwaveprop/canopy/bulk_fetch.ex
+++ b/lib/microwaveprop/canopy/bulk_fetch.ex
@ -0,0 +1,183 @@
+defmodule Microwaveprop.Canopy.BulkFetch do
+  @moduledoc """
+  One-shot prefetch of Lang et al. 2020 10 m global canopy-height COGs
+  (https://share.phys.ethz.ch/~pf/nlangdata/ETH_GlobalCanopyHeight_10m_2020_version1/3deg_cogs/),
+  sliced into per-degree uint8 tiles compatible with `Microwaveprop.Canopy`.
+
+  Run on prod (mix tasks don't work there) via:
+
+      bin/microwaveprop rpc 'Microwaveprop.Canopy.BulkFetch.run(32.8, -97.0, 200)'
+
+  Source COGs are 3° × 3° mosaics, ~150 MB each. Each is sliced into
+  the 1° × 1° uint8 .canopy tiles the runtime expects in
+  `/data/canopy/`. Requires `gdal_translate` in the runtime image.
+  """
+
+  alias Microwaveprop.Canopy
+
+  require Logger
+
+  @cog_base_url "https://share.phys.ethz.ch/~pf/nlangdata/ETH_GlobalCanopyHeight_10m_2020_version1/3deg_cogs"
+  @samples 3600
+
+  @doc """
+  Downloads + slices canopy tiles covering a `radius_mi` circle around
+  `(lat, lon)`. Skips per-degree tiles already on disk. Returns
+  `%{downloaded, sliced, cached, failed}` counters.
+  """
+  @spec run(float(), float(), pos_integer(), keyword()) :: map()
+  def run(lat, lon, radius_mi, opts \\ []) do
+    cache_dir = Keyword.get(opts, :cache_dir, default_cache_dir())
+    staging_dir = Keyword.get(opts, :staging_dir, default_staging_dir())
+    File.mkdir_p!(cache_dir)
+    File.mkdir_p!(staging_dir)
+
+    lat_deg = radius_mi / 69.0
+    lon_deg = radius_mi / (69.0 * max(:math.cos(lat * :math.pi() / 180.0), 0.1))
+
+    bbox = %{
+      south: lat - lat_deg,
+      north: lat + lat_deg,
+      west: lon - lon_deg,
+      east: lon + lon_deg
+    }
+
+    degrees = degree_tiles(bbox)
+    cog_tiles = Enum.uniq_by(degrees, &cog_tile_key/1)
+
+    Logger.info(
+      "canopy bulk fetch: lat=#{lat} lon=#{lon} radius_mi=#{radius_mi} " <>
+        "degree_tiles=#{length(degrees)} cog_tiles=#{length(cog_tiles)}"
+    )
+
+    state = %{downloaded: 0, sliced: 0, cached: 0, failed: 0}
+
+    state =
+      Enum.reduce(cog_tiles, state, fn deg, acc ->
+        process_cog_tile(deg, staging_dir, cache_dir, degrees, acc)
+      end)
+
+    Logger.info(
+      "canopy bulk fetch DONE: downloaded=#{state.downloaded} sliced=#{state.sliced} " <>
+        "cached=#{state.cached} failed=#{state.failed}"
+    )
+
+    Map.put(state, :cache_dir, cache_dir)
+  end
+
+  defp process_cog_tile({cog_lat, cog_lon}, staging_dir, cache_dir, all_degrees, acc) do
+    cog_path = Path.join(staging_dir, cog_filename(cog_lat, cog_lon))
+
+    case ensure_cog(cog_path, cog_lat, cog_lon) do
+      {:ok, :downloaded} ->
+        slice_cog(cog_path, cog_lat, cog_lon, cache_dir, all_degrees, %{acc | downloaded: acc.downloaded + 1})
+
+      {:ok, :cached} ->
+        slice_cog(cog_path, cog_lat, cog_lon, cache_dir, all_degrees, acc)
+
+      {:error, reason} ->
+        Logger.warning("canopy: skipping COG #{cog_filename(cog_lat, cog_lon)}: #{inspect(reason)}")
+        %{acc | failed: acc.failed + 1}
+    end
+  end
+
+  defp ensure_cog(path, cog_lat, cog_lon) do
+    if File.exists?(path) and File.stat!(path).size > 1_000_000 do
+      {:ok, :cached}
+    else
+      url = "#{@cog_base_url}/#{cog_filename(cog_lat, cog_lon)}"
+      Logger.info("canopy: downloading COG #{cog_filename(cog_lat, cog_lon)}")
+
+      case Req.get(url, into: File.stream!(path), retry: false, receive_timeout: 600_000) do
+        {:ok, %{status: 200}} -> {:ok, :downloaded}
+        {:ok, %{status: status}} -> {:error, {:http, status}}
+        {:error, reason} -> {:error, reason}
+      end
+    end
+  end
+
+  defp slice_cog(cog_path, cog_lat, cog_lon, cache_dir, all_degrees, acc) do
+    degrees_in_cog = Enum.filter(all_degrees, fn d -> cog_tile_key(d) == {cog_lat, cog_lon} end)
+    Enum.reduce(degrees_in_cog, acc, &slice_degree(&1, cog_path, cache_dir, &2))
+  end
+
+  defp slice_degree({lat, lon}, cog_path, cache_dir, acc) do
+    out_path = Path.join(cache_dir, Canopy.tile_filename(lat, lon))
+
+    cond do
+      File.exists?(out_path) -> %{acc | cached: acc.cached + 1}
+      match?(:ok, slice_one(cog_path, lat, lon, out_path)) -> %{acc | sliced: acc.sliced + 1}
+      true -> %{acc | failed: acc.failed + 1}
+    end
+  end
+
+  defp slice_one(cog_path, lat, lon, out_path) do
+    # gdal_translate -of ENVI -ot Byte -outsize 3600 3600 \
+    #   -projwin <ulx> <uly> <lrx> <lry> input.tif /tmp/N32W097.canopy
+    # ENVI puts the raw bytes in `out_path` plus sidecar `.hdr` files
+    # which we delete. Width before height = `-outsize 3600 3600`.
+    ulx = lon * 1.0
+    uly = (lat + 1) * 1.0
+    lrx = (lon + 1) * 1.0
+    lry = lat * 1.0
+
+    args = [
+      "-q",
+      "-of",
+      "ENVI",
+      "-ot",
+      "Byte",
+      "-outsize",
+      Integer.to_string(@samples),
+      Integer.to_string(@samples),
+      "-projwin",
+      Float.to_string(ulx),
+      Float.to_string(uly),
+      Float.to_string(lrx),
+      Float.to_string(lry),
+      cog_path,
+      out_path
+    ]
+
+    case System.cmd("gdal_translate", args, stderr_to_stdout: true) do
+      {_, 0} ->
+        # Discard ENVI sidecar metadata.
+        Enum.each([".hdr", ".aux.xml"], fn ext -> _ = File.rm(out_path <> ext) end)
+        :ok
+
+      {output, status} ->
+        Logger.warning("canopy: gdal_translate failed (status #{status}): #{output}")
+        _ = File.rm(out_path)
+        {:error, {:gdal, status}}
+    end
+  end
+
+  defp degree_tiles(bbox) do
+    south = floor(bbox.south)
+    north = floor(bbox.north)
+    west = floor(bbox.west)
+    east = floor(bbox.east)
+
+    for lat <- south..north, lon <- west..east, do: {lat, lon}
+  end
+
+  # Lang COGs are bucketed in 3° × 3° steps anchored at multiples of 3.
+  defp cog_tile_key({lat, lon}) do
+    lat_anchor = lat - Integer.mod(lat, 3)
+    lon_anchor = lon - Integer.mod(lon, 3)
+    {lat_anchor, lon_anchor}
+  end
+
+  defp cog_filename(lat, lon) do
+    lat_prefix = if lat >= 0, do: "N", else: "S"
+    lon_prefix = if lon >= 0, do: "E", else: "W"
+    lat_str = lat |> abs() |> Integer.to_string() |> String.pad_leading(2, "0")
+    lon_str = lon |> abs() |> Integer.to_string() |> String.pad_leading(3, "0")
+
+    "ETH_GlobalCanopyHeight_10m_2020_#{lat_prefix}#{lat_str}#{lon_prefix}#{lon_str}_Map.tif"
+  end
+
+  defp default_cache_dir, do: Application.get_env(:microwaveprop, :canopy_tiles_dir, "/data/canopy")
+
+  defp default_staging_dir, do: Application.get_env(:microwaveprop, :canopy_staging_dir, "/data/canopy/staging")
+end
--- a/lib/microwaveprop/rover/compute.ex
+++ b/lib/microwaveprop/rover/compute.ex
@ -9,6 +9,7 @@ defmodule Microwaveprop.Rover.Compute do

  alias Microwaveprop.Buildings.Index, as: BuildingsIndex
  alias Microwaveprop.Buildings.Loader, as: BuildingsLoader
+  alias Microwaveprop.Canopy
  alias Microwaveprop.Propagation
  alias Microwaveprop.Radio.Maidenhead
  alias Microwaveprop.Rover.Aggregator
@ -56,6 +57,15 @@ defmodule Microwaveprop.Rover.Compute do
  @building_clutter_m_per_db 5.0
  @building_clutter_cap_db 6.0

+  # Per-cell tree-canopy clutter penalty. Same shape as the building
+  # clutter penalty but scaled gentler (1 dB / 8 m, capped at 4 dB) —
+  # forests are common away from roads so the penalty is a tiebreaker,
+  # not a hard exclusion. Path-clearance already drops cells whose
+  # link to a station is blocked by trees; this penalty surfaces the
+  # "I'm in a forest, every direction is foliage" signal.
+  @canopy_clutter_m_per_db 8.0
+  @canopy_clutter_cap_db 4.0
+
  @tier_excellent 10.0
  @tier_good 3.0
  @tier_marginal 0.0
@ -89,6 +99,9 @@ defmodule Microwaveprop.Rover.Compute do
    buildings_clutter_lookup =
      Keyword.get(deps, :buildings_clutter_lookup, &default_building_clutter/1)

+    canopy_clutter_lookup =
+      Keyword.get(deps, :canopy_clutter_lookup, &default_canopy_clutter/1)
+
    hilltop_snap = Keyword.get(deps, :hilltop_snap, &Hilltop.snap/1)
    progress = Keyword.get(deps, :progress, fn _, _, _ -> :ok end)

@ -107,7 +120,7 @@ defmodule Microwaveprop.Rover.Compute do
    bbox = bbox_around(home, radius_km)

    road_enabled? = Application.get_env(:microwaveprop, :rover_road_proximity_enabled, true)
-    total_steps = if road_enabled?, do: 8, else: 7
+    total_steps = if road_enabled?, do: 9, else: 8
    counter = :counters.new(1, [:atomics])

    step = fn label ->
@ -147,6 +160,8 @@ defmodule Microwaveprop.Rover.Compute do

    step.("Scanning building clutter")
    clutter_map = time_step("building_clutter", fn -> buildings_clutter_lookup.(in_radius) end)
+    step.("Scanning tree canopy")
+    canopy_map = time_step("canopy_clutter", fn -> canopy_clutter_lookup.(in_radius) end)
    step.("Scoring cells")

    home_elev = home.elev_m || 0
@ -156,7 +171,8 @@ defmodule Microwaveprop.Rover.Compute do
      clearance: clearance_map,
      prominence: prominence_map,
      road: road_map,
-      clutter: clutter_map
+      clutter: clutter_map,
+      canopy: canopy_map
    }

    cells =
@ -213,13 +229,15 @@ defmodule Microwaveprop.Rover.Compute do
      clearance: clearance_map,
      prominence: prominence_map,
      road: road_map,
-      clutter: clutter_map
+      clutter: clutter_map,
+      canopy: canopy_map
    } = cell_maps

    elev_m = Map.get(elev_map, {cell.lat, cell.lon})
    prominence_m = Map.get(prominence_map, {cell.lat, cell.lon})
    road_km = Map.get(road_map, {cell.lat, cell.lon})
    building_height_m = Map.get(clutter_map, {cell.lat, cell.lon})
+    canopy_height_m = Map.get(canopy_map, {cell.lat, cell.lon})
    base_margin = LinkMargin.link_margin_from_score(cell.score, mode)

    margins =
@ -239,7 +257,8 @@ defmodule Microwaveprop.Rover.Compute do

        score =
          agg_db + prominence_db(prominence_m) - drive_penalty(dist_km) -
-            road_penalty_db(road_km) - building_clutter_db(building_height_m)
+            road_penalty_db(road_km) - building_clutter_db(building_height_m) -
+            canopy_clutter_db(canopy_height_m)

        %{
          lat: cell.lat,
@ -248,6 +267,7 @@ defmodule Microwaveprop.Rover.Compute do
          prominence_m: prominence_m,
          road_km: road_km,
          building_height_m: building_height_m,
+          canopy_height_m: canopy_height_m,
          score: score,
          distance_km: dist_km,
          tier_color: tier_color(score)
@ -291,6 +311,20 @@ defmodule Microwaveprop.Rover.Compute do
    end)
  end

+  defp canopy_clutter_db(nil), do: 0.0
+  defp canopy_clutter_db(h) when is_number(h) and h <= 0, do: 0.0
+
+  defp canopy_clutter_db(h) when is_number(h) do
+    db = h / @canopy_clutter_m_per_db
+    min(db, @canopy_clutter_cap_db)
+  end
+
+  defp default_canopy_clutter(cells) do
+    cells
+    |> Enum.map(fn cell -> {cell.lat, cell.lon} end)
+    |> Canopy.lookup_many()
+  end
+
  defp terrain_db(nil), do: 0.0

  defp terrain_db(clearance_m) do