diff --git a/Dockerfile b/Dockerfile index 96395d67..b2c4938e 100644 --- a/Dockerfile +++ b/Dockerfile @@ -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 diff --git a/lib/microwaveprop/canopy/bulk_fetch.ex b/lib/microwaveprop/canopy/bulk_fetch.ex new file mode 100644 index 00000000..a9e37f4a --- /dev/null +++ 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 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 diff --git a/lib/microwaveprop/rover/compute.ex b/lib/microwaveprop/rover/compute.ex index 66cc9058..c420bb0d 100644 --- 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