diff --git a/assets/js/elevation_profile_hook.ts b/assets/js/elevation_profile_hook.ts index c31b451d..6f144e7d 100644 --- a/assets/js/elevation_profile_hook.ts +++ b/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 } ] : []), { diff --git a/lib/microwaveprop/canopy.ex b/lib/microwaveprop/canopy.ex new file mode 100644 index 00000000..7784061b --- /dev/null +++ b/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 + _ -> 0 + end + end + + defp default_dir do + Application.get_env(:microwaveprop, :canopy_tiles_dir, "/data/canopy") + end +end diff --git a/lib/microwaveprop/rover/candidate_detail.ex b/lib/microwaveprop/rover/candidate_detail.ex index 06005711..b51c62b4 100644 --- a/lib/microwaveprop/rover/candidate_detail.ex +++ b/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) diff --git a/lib/microwaveprop/rover/path_terrain.ex b/lib/microwaveprop/rover/path_terrain.ex index 652d2c00..9a618261 100644 --- a/lib/microwaveprop/rover/path_terrain.ex +++ b/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). diff --git a/lib/microwaveprop_web/live/path_live.ex b/lib/microwaveprop_web/live/path_live.ex index a420d349..9dce28f6 100644 --- a/lib/microwaveprop_web/live/path_live.ex +++ b/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) diff --git a/lib/microwaveprop_web/live/rover_live.ex b/lib/microwaveprop_web/live/rover_live.ex index 563f4c00..7bc07428 100644 --- a/lib/microwaveprop_web/live/rover_live.ex +++ b/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 = - if has_buildings? do - building_line = - Enum.map_join(points, " ", fn pt -> - x = pt.dist_km / total_km * 240.0 - top = pt.elev + round(Map.get(pt, :building_m, 0.0)) - y = 60.0 - (top - min_e) / span * 50.0 - "#{Float.round(x, 2)},#{Float.round(y, 2)}" - end) + terrain_back = + points + |> Enum.reverse() + |> Enum.map_join(" ", fn pt -> + x = pt.dist_km / total_km * 240.0 + y = 60.0 - (pt.elev - min_e) / span * 50.0 + "#{Float.round(x, 2)},#{Float.round(y, 2)}" + end) - # Close back along the terrain so the building polygon sits on top. - terrain_back = - points - |> Enum.reverse() - |> Enum.map_join(" ", fn pt -> - x = pt.dist_km / total_km * 240.0 - y = 60.0 - (pt.elev - min_e) / span * 50.0 - "#{Float.round(x, 2)},#{Float.round(y, 2)}" - end) + obstacle_fill = fn obstacle_key -> + obstacle_line = + Enum.map_join(points, " ", fn pt -> + x = pt.dist_km / total_km * 240.0 + top = pt.elev + round(Map.get(pt, obstacle_key, 0)) + y = 60.0 - (top - min_e) / span * 50.0 + "#{Float.round(x, 2)},#{Float.round(y, 2)}" + end) - building_line <> " " <> terrain_back - end + obstacle_line <> " " <> terrain_back + 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""" + 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