diff --git a/lib/towerops/coverages/profile.ex b/lib/towerops/coverages/profile.ex index 901492ff..66121db7 100644 --- a/lib/towerops/coverages/profile.ex +++ b/lib/towerops/coverages/profile.ex @@ -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 diff --git a/lib/towerops/coverages/tree_canopy.ex b/lib/towerops/coverages/tree_canopy.ex new file mode 100644 index 00000000..b97ae855 --- /dev/null +++ b/lib/towerops/coverages/tree_canopy.ex @@ -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 diff --git a/lib/towerops/workers/coverage_worker.ex b/lib/towerops/workers/coverage_worker.ex index 74640924..f58c8793 100644 --- a/lib/towerops/workers/coverage_worker.ex +++ b/lib/towerops/workers/coverage_worker.ex @@ -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). diff --git a/test/test_helper.exs b/test/test_helper.exs index 3e3011d0..9d0245ec 100644 --- a/test/test_helper.exs +++ b/test/test_helper.exs @@ -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 diff --git a/test/towerops/coverages/tree_canopy_test.exs b/test/towerops/coverages/tree_canopy_test.exs new file mode 100644 index 00000000..4376f6c7 --- /dev/null +++ b/test/towerops/coverages/tree_canopy_test.exs @@ -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