feat(coverage): tree-canopy clutter via LANDFIRE EVH

Closes the last clutter source. NLOS heatmaps now see trees as
obstacles in addition to buildings:

- New Towerops.Coverages.TreeCanopy module fetches an AAIGrid of
  canopy heights for the coverage bbox via the same /vsicurl/
  + Lidar.Reader pipeline used for terrain. Source defaults to
  the public LANDFIRE 2022 EVH COG (~30 m, all CONUS) and is
  override-able via :tree_canopy_url.
- EVH cells are integer codes; decode_evh maps the canopy range
  (101..199 → height-100 m), the herbaceous range (11..29 →
  (code-10) * 0.1 m), and zeroes everything else.
- Profile.sample/5 takes a :canopy grid alongside :clutter and
  uses max(building_height, canopy_height) so a tree poking above
  a one-storey roof still shadows.
- CoverageWorker fetches the canopy grid once per job (best
  effort — failures log + continue with nil) and threads it
  through the env map. LOS keeps clutter empty + canopy nil
  (height-above-clutter view); NLOS gets both.
- Tests: tree_canopy_test for the disabled / override paths;
  test_helper disables the LANDFIRE fetch by default so the
  rest of the suite doesn't stream remote bytes.
This commit is contained in:
Graham McIntire 2026-05-06 16:54:59 -05:00
parent ce9bd744c7
commit 6e4b07026f
5 changed files with 190 additions and 26 deletions

View file

@ -59,29 +59,38 @@ defmodule Towerops.Coverages.Profile do
@doc """
Same as `sample/4` but accepts `:clutter` (a list of building polygons
the caller has prefetched for this path). Each polygon must look like
`%{height_m: float(), coords: [{lon, lat}, ...]}` closed ring,
WGS84. Sample points that fall inside a polygon get the building's
rooftop height added to terrain so the diffraction model "sees" the
obstacle.
the caller has prefetched for this path) and `:canopy` (an AAIGrid of
tree-canopy heights in metres, e.g. from
`Towerops.Coverages.TreeCanopy.fetch_grid/2`).
This is the cnHeat NLOS-mode input. Pass `clutter: []` (or omit) for
bare-terrain LOS sampling.
At each sample point the DSM is `terrain + max(building, canopy)`
the tallest of the two clutter sources wins so a tree poking above a
one-storey building still shadows. Pass empty / nil for either to
disable that clutter input.
"""
@spec sample(map(), {number(), number()}, {number(), number()}, pos_integer(), keyword()) ::
[float()]
def sample(grid, {from_lat, from_lon}, _to, 1, opts) do
[dsm_height(grid, from_lat, from_lon, Keyword.get(opts, :clutter, []))]
[
dsm_height(
grid,
from_lat,
from_lon,
Keyword.get(opts, :clutter, []),
Keyword.get(opts, :canopy)
)
]
end
def sample(grid, {from_lat, from_lon}, {to_lat, to_lon}, n, opts) when n >= 2 do
clutter = Keyword.get(opts, :clutter, [])
canopy = Keyword.get(opts, :canopy)
Enum.map(0..(n - 1), fn i ->
t = i / (n - 1)
lat = from_lat + (to_lat - from_lat) * t
lon = from_lon + (to_lon - from_lon) * t
dsm_height(grid, lat, lon, clutter)
dsm_height(grid, lat, lon, clutter, canopy)
end)
end
@ -112,14 +121,23 @@ defmodule Towerops.Coverages.Profile do
end
end
# Returns the DSM height at (lat, lon) — terrain plus the rooftop
# of any building polygon that contains the point. Tree canopy
# could plug in here once a height-keyed source is wired up.
defp dsm_height(grid, lat, lon, []), do: elevation_or_zero(grid, lat, lon)
defp dsm_height(grid, lat, lon, clutter) do
# Returns the DSM height at (lat, lon) — terrain plus the tallest
# of the rooftop / tree-canopy clutter at that point. Either
# clutter source may be empty / nil to disable.
defp dsm_height(grid, lat, lon, clutter, canopy) do
base = elevation_or_zero(grid, lat, lon)
base + clutter_at(clutter, lat, lon)
building = clutter_at(clutter, lat, lon)
canopy_h = canopy_at(canopy, lat, lon)
base + max(building, canopy_h)
end
defp canopy_at(nil, _lat, _lon), do: 0.0
defp canopy_at(grid, lat, lon) when is_map(grid) do
case elevation_at(grid, lat, lon) do
h when is_number(h) and h > 0 -> h
_ -> 0.0
end
end
defp clutter_at([], _lat, _lon), do: 0.0

View file

@ -0,0 +1,88 @@
defmodule Towerops.Coverages.TreeCanopy do
@moduledoc """
Tree canopy height source for coverage compute.
We use LANDFIRE Existing Vegetation Height (EVH) a 30 m nationwide
raster published as a Cloud-Optimized GeoTIFF by USFS. It encodes
vegetation heights in metres in the upper code range (101-199
height = code - 100), with low codes reserved for non-tree classes.
Default URL points to the public LANDFIRE 2022 release. Overridable
via `config :towerops, :tree_canopy_url, "https://..."` so deployments
can point at a local mirror or a more recent release.
Like LIDAR DEMs, we never mirror raster bytes `gdal_translate` reads
exactly the bytes needed for a bbox via `/vsicurl/` HTTP byte ranges.
"""
alias Towerops.Lidar.Reader
alias Towerops.Lidar.Tile
@default_url "https://landfire.gov/data/LF2022/LF20_EVH_220.tif"
@doc """
Returns an AAIGrid of canopy heights (metres) covering the bbox at
the requested cell size. Returns `{:error, _}` if the canopy
source is unreachable or the bbox falls outside CONUS.
Heights are derived from the LANDFIRE EVH code:
* 101..199 tree canopy, height = code - 100 metres
* 11..29 herbaceous / shrub, height = (code - 10) * 0.1 m
* else 0 (non-vegetation, water, bare ground, etc.)
Callers receive a grid whose cells are already height-in-metres so
downstream sampling treats it identically to terrain.
"""
@spec fetch_grid({number(), number(), number(), number()}, number()) ::
{:ok, map()} | {:error, term()}
def fetch_grid(bbox, cell_size_deg) do
if enabled?() do
do_fetch(bbox, cell_size_deg)
else
{:error, :disabled}
end
end
@doc "Returns the URL configured for the canopy raster."
@spec source_url() :: String.t()
def source_url do
Application.get_env(:towerops, :tree_canopy_url, @default_url)
end
defp enabled? do
Application.get_env(:towerops, :tree_canopy_enabled, true)
end
defp do_fetch(bbox, cell_size_deg) do
tile = %Tile{
source: "landfire_evh",
project_name: "LANDFIRE_EVH",
tile_name: "evh",
url: source_url(),
crs: "EPSG:5070",
availability: "available"
}
case Reader.elevation_grid(tile, bbox, cell_size_deg) do
{:ok, grid} -> {:ok, decode_evh(grid)}
{:error, reason} -> {:error, reason}
end
end
# LANDFIRE EVH cells are integer codes; rewrite them as height-in-m.
defp decode_evh(grid) do
decoded =
Enum.map(grid.cells, fn row ->
Enum.map(row, fn
v when is_number(v) -> evh_to_height(v)
other -> other
end)
end)
%{grid | cells: decoded, nodata_value: -9999.0}
end
defp evh_to_height(code) when code >= 101 and code <= 199, do: (code - 100) * 1.0
defp evh_to_height(code) when code >= 11 and code <= 29, do: (code - 10) * 0.1
defp evh_to_height(_), do: 0.0
end

View file

@ -36,6 +36,7 @@ defmodule Towerops.Workers.CoverageWorker do
alias Towerops.Coverages.Profile
alias Towerops.Coverages.Propagation
alias Towerops.Coverages.Raster
alias Towerops.Coverages.TreeCanopy
alias Towerops.Repo
require Logger
@ -81,27 +82,46 @@ defmodule Towerops.Workers.CoverageWorker do
_ = update_progress(coverage, "computing", 5),
{:ok, grid} <- fetch_terrain(bbox, coverage.cell_size_m, lat),
clutter = Buildings.for_bbox(bbox),
canopy = fetch_canopy(bbox, coverage.cell_size_m, lat),
_ = update_progress(coverage, "computing", 15),
{:ok, tiers} <- compute_all_tiers(coverage, antenna, lat, lon, bbox, grid, clutter) do
{:ok, tiers} <-
compute_all_tiers(coverage, antenna, lat, lon, bbox, grid, clutter, canopy) do
finalize(coverage, tiers)
else
{:error, reason} -> fail(coverage, reason)
end
end
# Best-effort canopy fetch — non-fatal if disabled / out-of-area.
# Only NLOS sampling consumes it, so a nil here just means
# "compute as if no trees are present" for that mode.
defp fetch_canopy(bbox, cell_size_m, lat) do
cell_size_deg = cell_size_m / (111_000.0 * :math.cos(lat * :math.pi() / 180.0))
case TreeCanopy.fetch_grid(bbox, cell_size_deg) do
{:ok, grid} ->
grid
{:error, reason} ->
Logger.info("CoverageWorker: tree canopy unavailable (#{inspect(reason)}); skipping")
nil
end
end
# Runs compute_pixels twice per SM-height tier (LOS / NLOS) and
# writes per-tier raster + PNG sets. LOS treats clutter as ground
# baseline (cnHeat "Height above clutter"); NLOS samples the
# prefetched building list into the path profile so rooftop
# diffraction shadows show up ("Height above ground").
defp compute_all_tiers(coverage, antenna, lat, lon, bbox, grid, clutter) do
defp compute_all_tiers(coverage, antenna, lat, lon, bbox, grid, clutter, canopy) do
indexed = Enum.with_index(@sm_height_tiers_m)
initial = {%{los: %{}, nlos: %{}}, nil}
env = %{antenna: antenna, lat: lat, lon: lon, bbox: bbox, grid: grid, clutter: clutter, canopy: canopy}
{results, error} =
Enum.reduce_while(indexed, initial, fn {height_m, idx}, {acc, _last} ->
case compute_tier(coverage, antenna, lat, lon, bbox, grid, clutter, height_m) do
case compute_tier(coverage, env, height_m) do
{:ok, %{los: los_paths, nlos: nlos_paths}} ->
update_progress(
coverage,
@ -133,16 +153,16 @@ defmodule Towerops.Workers.CoverageWorker do
end
end
defp compute_tier(coverage, antenna, lat, lon, bbox, grid, clutter, height_m) do
defp compute_tier(coverage, env, height_m) do
cov = %{coverage | receiver_height_m: height_m}
los_ctx = %{antenna: antenna, grid: grid, bbox: bbox, clutter: []}
nlos_ctx = %{antenna: antenna, grid: grid, bbox: bbox, clutter: clutter}
los_ctx = Map.merge(env, %{clutter: [], canopy: nil})
nlos_ctx = Map.merge(env, %{clutter: env.clutter, canopy: env.canopy})
with {:ok, los_out} <- compute_pixels(cov, los_ctx, lat, lon),
with {:ok, los_out} <- compute_pixels(cov, los_ctx, env.lat, env.lon),
{:ok, los_paths} <-
Raster.write(coverage, los_out.pixels, los_out.dims, tier_suffix(height_m, :los)),
{:ok, nlos_out} <- compute_pixels(cov, nlos_ctx, lat, lon),
{:ok, nlos_out} <- compute_pixels(cov, nlos_ctx, env.lat, env.lon),
{:ok, nlos_paths} <-
Raster.write(coverage, nlos_out.pixels, nlos_out.dims, tier_suffix(height_m, :nlos)) do
{:ok, %{los: los_paths, nlos: nlos_paths}}
@ -264,10 +284,11 @@ defmodule Towerops.Workers.CoverageWorker do
distance_m
) do
%{antenna: antenna, grid: grid, clutter: clutter, eirp: eirp} = ctx
canopy = Map.get(ctx, :canopy)
# 1. Path profile (terrain + building clutter from TX → RX).
# 1. Path profile (terrain + building/tree clutter from TX → RX).
samples = max(3, min(@profile_samples, ceil(distance_m / coverage.cell_size_m) + 2))
profile = Profile.sample(grid, tx_latlon, rx_latlon, samples, clutter: clutter)
profile = Profile.sample(grid, tx_latlon, rx_latlon, samples, clutter: clutter, canopy: canopy)
# Replace TX endpoint with the ground elevation under the antenna (so the
# path loss model sees the correct antenna height above terrain).

View file

@ -86,6 +86,11 @@ Mox.defmock(Towerops.RateLimit.RedisClientMock, for: Towerops.RateLimit.RedisCli
# enables this fallback for nationwide CONUS coverage.
Application.put_env(:towerops, :lidar_ned_fallback, false)
# Same idea for the LANDFIRE tree-canopy COG: tests disable it so the
# worker doesn't try to stream LANDFIRE bytes during compute. Tests
# that exercise the canopy code path opt back in explicitly.
Application.put_env(:towerops, :tree_canopy_enabled, false)
# Define stub module for SNMP mock (returns errors for all calls)
defmodule Towerops.Snmp.SnmpMockStub do
@moduledoc false

View file

@ -0,0 +1,32 @@
defmodule Towerops.Coverages.TreeCanopyTest do
use ExUnit.Case, async: false
alias Towerops.Coverages.TreeCanopy
setup do
on_exit(fn ->
Application.delete_env(:towerops, :tree_canopy_enabled)
Application.delete_env(:towerops, :tree_canopy_url)
end)
end
describe "fetch_grid/2" do
test "returns :disabled when the canopy source is turned off" do
Application.put_env(:towerops, :tree_canopy_enabled, false)
assert {:error, :disabled} =
TreeCanopy.fetch_grid({-97.8, 30.2, -97.7, 30.3}, 0.001)
end
end
describe "source_url/0" do
test "defaults to the public LANDFIRE EVH COG" do
assert TreeCanopy.source_url() =~ "landfire.gov"
end
test "respects the runtime override" do
Application.put_env(:towerops, :tree_canopy_url, "https://example.com/canopy.tif")
assert TreeCanopy.source_url() == "https://example.com/canopy.tif"
end
end
end