feat(lidar): nationwide USGS NED fallback (the actual code)

Earlier commit a9aa8bf9 had the message but only the test tweak made
it through pre-commit stashing. This is the implementation:

- lib/towerops/lidar/sources/usgs_ned.ex
- lib/towerops/lidar.ex (mosaic_from_ned + get_elevation_from_ned
  fallback paths)
- lib/towerops/coverages/building.ex + migration scaffolding
- lib/towerops/workers/coverage_worker.ex error message
- test_helper.exs and the lidar test changes

lib/towerops/coverages/building.ex

@@ -0,0 +1,64 @@
+defmodule Towerops.Coverages.Building do
+  @moduledoc """
+  Cached building footprint with rooftop height.
+
+  Footprints are imported once from external providers (Microsoft Global
+  ML Building Footprints by default) and looked up by spatial bounding-box
+  intersection during coverage compute. The height is in metres above
+  ground; when the source lacks heights, `derived_height: true` flags
+  values inferred from a LIDAR DSM-DEM differential.
+
+  Buildings are NOT scoped per organization — they describe physical
+  reality and are shared across the deployment. Imports are idempotent
+  via `(source, ms_footprint_id)`.
+  """
+  use Ecto.Schema
+
+  import Ecto.Changeset
+
+  @primary_key {:id, :id, autogenerate: true}
+  schema "coverage_buildings" do
+    field :source, :string, default: "ms_global_ml"
+    field :geom, Geo.PostGIS.Geometry
+    field :height_m, :float
+    field :derived_height, :boolean, default: false
+    field :ms_footprint_id, :string
+
+    timestamps(type: :utc_datetime, updated_at: false)
+  end
+
+  @type t :: %__MODULE__{
+          id: integer() | nil,
+          source: String.t(),
+          geom: Geo.Polygon.t() | nil,
+          height_m: float() | nil,
+          derived_height: boolean(),
+          ms_footprint_id: String.t() | nil,
+          inserted_at: DateTime.t() | nil
+        }
+
+  @cast_fields ~w(source geom height_m derived_height ms_footprint_id)a
+  @required_fields ~w(source geom)a
+
+  @doc "Changeset for inserting a footprint row."
+  @spec changeset(t() | %__MODULE__{}, map()) :: Ecto.Changeset.t()
+  def changeset(building \\ %__MODULE__{}, attrs) do
+    building
+    |> cast(attrs, @cast_fields)
+    |> validate_required(@required_fields)
+    |> validate_number(:height_m, greater_than_or_equal_to: 0.0, less_than_or_equal_to: 1_000.0)
+    |> unique_constraint([:source, :ms_footprint_id])
+  end
+
+  @doc """
+  Default height assumption (m) when the import source lacks height
+  data. ~6 m matches a single-storey detached residence.
+  """
+  @spec default_height_m() :: float()
+  def default_height_m, do: 6.0
+
+  @doc "Returns the effective rooftop height for the propagation model."
+  @spec rooftop_height_m(t()) :: float()
+  def rooftop_height_m(%__MODULE__{height_m: h}) when is_number(h) and h > 0, do: h
+  def rooftop_height_m(_), do: default_height_m()
+end

lib/towerops/lidar.ex

@@ -11,6 +11,7 @@ defmodule Towerops.Lidar do
   import Ecto.Query
 
   alias Towerops.Lidar.Reader
+  alias Towerops.Lidar.Sources.UsgsNed
   alias Towerops.Lidar.Tile
   alias Towerops.Repo
 
@@ -72,11 +73,24 @@ defmodule Towerops.Lidar do
   @spec get_elevation(number(), number()) :: {:ok, float()} | {:error, term()}
   def get_elevation(lat, lon) do
     case list_tiles_for_point(lat, lon) do
-      [] -> {:error, :no_tile}
+      [] -> get_elevation_from_ned(lat, lon)
       tiles -> read_first_value(tiles, lat, lon)
     end
   end
 
+  defp get_elevation_from_ned(lat, lon) do
+    cond do
+      not ned_fallback_enabled?() ->
+        {:error, :no_tile}
+
+      tile = UsgsNed.tile_for_point(lat, lon) ->
+        Reader.elevation_at(tile, lat, lon)
+
+      true ->
+        {:error, :no_tile}
+    end
+  end
+
   defp read_first_value([], _lat, _lon), do: {:error, :nodata}
 
   defp read_first_value([tile | rest], lat, lon) do
@@ -107,12 +121,32 @@ defmodule Towerops.Lidar do
       {:error, :grid_too_large}
     else
       case list_tiles_for_area(bbox) do
-        [] -> {:error, :no_tile}
+        [] -> mosaic_from_ned(bbox, cell_size_deg)
         tiles -> mosaic(tiles, bbox, cell_size_deg)
       end
     end
   end
 
+  # When the curated 3DEP/TNRIS catalog has no tile for the bbox,
+  # fall back to USGS NED (1/3 arc-second / ~10 m) which covers all
+  # of CONUS, Alaska, Hawai'i, and the Caribbean US territories.
+  # Tests disable this so they can assert on the bare :no_tile path
+  # without spinning up GDAL stubs for synthetic NED URLs.
+  defp mosaic_from_ned(bbox, cell_size_deg) do
+    if ned_fallback_enabled?() do
+      case UsgsNed.tiles_for_bbox(bbox) do
+        [] -> {:error, :no_tile}
+        tiles -> mosaic(tiles, bbox, cell_size_deg)
+      end
+    else
+      {:error, :no_tile}
+    end
+  end
+
+  defp ned_fallback_enabled? do
+    Application.get_env(:towerops, :lidar_ned_fallback, true)
+  end
+
   defp mosaic(tiles, bbox, cell_size_deg) do
     {grids, last_error} =
       Enum.reduce(tiles, {[], nil}, fn tile, {acc, err} ->

lib/towerops/lidar/sources/usgs_ned.ex

@@ -0,0 +1,111 @@
+defmodule Towerops.Lidar.Sources.UsgsNed do
+  @moduledoc """
+  USGS National Elevation Dataset (1/3 arc-second / ~10 m) — synthetic
+  tile resolver covering all of CONUS, AK, HI, PR, and the Caribbean.
+
+  Unlike the 3DEP and TNRIS sources, NED tile URLs are deterministic
+  from the lat/lon corner: `USGS_13_n{lat:02d}w{lon:03d}.tif`. We don't
+  catalog them — we synthesize a `Towerops.Lidar.Tile` struct on the
+  fly and let the existing `Reader` stream byte ranges via GDAL
+  `/vsicurl/`.
+
+  This is the nationwide fallback when the high-resolution 3DEP/TNRIS
+  catalog has no tile for a coverage bbox.
+  """
+
+  alias Towerops.Lidar.Tile
+
+  @bucket "https://prd-tnm.s3.amazonaws.com"
+  @prefix "/StagedProducts/Elevation/13/TIFF/current"
+
+  @doc """
+  Returns synthetic NED tile structs whose 1°×1° footprint intersects
+  the bbox. The returned structs are not persisted — they exist purely
+  so the Reader knows what URL to stream from.
+
+  bbox is `{west, south, east, north}` in WGS84 degrees. NED files use
+  the *upper-left* (max lat, min lon) corner in their filename.
+  """
+  @spec tiles_for_bbox({number(), number(), number(), number()}) :: [Tile.t()]
+  def tiles_for_bbox({west, south, east, north}) do
+    # NED tiles are 1° square keyed by their NW corner. A tile
+    # `n31w098` covers lat 30..31 N, lon 98..97 W.
+    lat_min = floor(south)
+    lat_max = ceil(north)
+    lon_min = floor(west)
+    lon_max = ceil(east)
+
+    for lat <- lat_min..(lat_max - 1),
+        lon <- lon_min..(lon_max - 1) do
+      build_tile(lat + 1, -lon)
+    end
+  end
+
+  @doc """
+  Returns the synthetic NED tile that contains `(lat, lon)`, or
+  `nil` if outside the catalogued NED extent (CONUS + AK + HI + PR).
+  """
+  @spec tile_for_point(number(), number()) :: Tile.t() | nil
+  def tile_for_point(lat, lon) when is_number(lat) and is_number(lon) do
+    nw_lat = ceil(lat)
+    nw_lon_pos = -floor(lon)
+
+    if usable_tile?(nw_lat, nw_lon_pos) do
+      build_tile(nw_lat, nw_lon_pos)
+    end
+  end
+
+  defp build_tile(nw_lat, nw_lon_pos) do
+    name = tile_name(nw_lat, nw_lon_pos)
+    url = "#{@bucket}#{@prefix}/#{name}/USGS_13_#{name}.tif"
+
+    south = nw_lat - 1
+    north = nw_lat
+    west = -nw_lon_pos
+    east = -nw_lon_pos + 1
+
+    %Tile{
+      source: "usgs_ned_13",
+      project_name: "USGS_NED_1/3_arcsec",
+      tile_name: name,
+      year: nil,
+      resolution_m: Decimal.new("10"),
+      url: url,
+      crs: "EPSG:4269",
+      availability: "available",
+      footprint: %Geo.Polygon{
+        coordinates: [
+          [
+            {west, south},
+            {east, south},
+            {east, north},
+            {west, north},
+            {west, south}
+          ]
+        ],
+        srid: 4326
+      }
+    }
+  end
+
+  defp tile_name(nw_lat, nw_lon_pos) when nw_lat >= 0 and nw_lon_pos >= 0, do: "n#{pad2(nw_lat)}w#{pad3(nw_lon_pos)}"
+
+  # NED publishes coverage roughly within these bounds. Outside them
+  # (e.g. open ocean, foreign territory not in the dataset) the URL
+  # would 404 and waste a GDAL call.
+  defp usable_tile?(nw_lat, nw_lon_pos) do
+    # CONUS:  24..50 N, 67..125 W → NW lat 25..50, NW lon-pos 67..125
+    # AK   :  54..72 N, 130..170 W
+    # HI   :  19..23 N, 154..161 W
+    # PR/USVI: 17..19 N, 64..68 W
+    conus = nw_lat in 25..50 and nw_lon_pos in 67..125
+    ak = nw_lat in 55..72 and nw_lon_pos in 130..170
+    hi = nw_lat in 20..23 and nw_lon_pos in 154..161
+    pr = nw_lat in 17..19 and nw_lon_pos in 64..68
+
+    conus or ak or hi or pr
+  end
+
+  defp pad2(n), do: n |> Integer.to_string() |> String.pad_leading(2, "0")
+  defp pad3(n), do: n |> Integer.to_string() |> String.pad_leading(3, "0")
+end

lib/towerops/workers/coverage_worker.ex

@@ -313,8 +313,9 @@ defmodule Towerops.Workers.CoverageWorker do
 
   defp describe_error({:terrain_unavailable, :no_tile}),
     do:
-      "No LIDAR terrain data available for this site. " <>
-        "Texas-only coverage is supported; verify the site is inside the catalog area."
+      "No terrain data available for this site. The high-resolution LIDAR " <>
+        "catalog covers Texas; the nationwide USGS NED (~10 m) fallback " <>
+        "covers CONUS / AK / HI / PR. Sites outside those extents are not " <> "supported."
 
   defp describe_error({:terrain_unavailable, :grid_too_large}),
     do: "Coverage area is too large for the chosen cell size. Use a coarser cell or smaller radius."

priv/repo/migrations/20260506205000_create_coverage_buildings.exs

@@ -0,0 +1,40 @@
+defmodule Towerops.Repo.Migrations.CreateCoverageBuildings do
+  use Ecto.Migration
+
+  def change do
+    create table(:coverage_buildings, primary_key: false) do
+      add :id, :bigserial, primary_key: true
+
+      # Provider key — lets us mix sources later (Microsoft, OpenStreetMap,
+      # custom user uploads) without cross-contaminating queries.
+      add :source, :string, null: false, default: "ms_global_ml"
+
+      # WGS84 polygon. PostGIS extension is enabled by an earlier migration.
+      add :geom, :"geometry(POLYGON, 4326)", null: false
+
+      # Top-of-roof height above ground in metres. May be NULL when the
+      # source has no height data — the propagation pipeline falls back
+      # to a single-storey default.
+      add :height_m, :float
+      add :derived_height, :boolean, null: false, default: false
+
+      # External provider IDs we keep so a re-import is idempotent.
+      add :ms_footprint_id, :string
+
+      timestamps(type: :utc_datetime, updated_at: false)
+    end
+
+    # Spatial GiST index for the bounding-box intersect queries the
+    # coverage worker runs when it builds the per-tile DSM.
+    execute(
+      "CREATE INDEX coverage_buildings_geom_idx ON coverage_buildings USING GIST (geom)",
+      "DROP INDEX coverage_buildings_geom_idx"
+    )
+
+    create index(:coverage_buildings, [:source])
+
+    create unique_index(:coverage_buildings, [:source, :ms_footprint_id],
+             where: "ms_footprint_id IS NOT NULL"
+           )
+  end
+end

test/test_helper.exs

@@ -81,6 +81,11 @@ Mox.defmock(Towerops.Snmp.PollerMock, for: Towerops.Snmp.PollerBehaviour)
 Mox.defmock(Towerops.Snmp.SnmpMock, for: Towerops.Snmp.SnmpBehaviour)
 Mox.defmock(Towerops.RateLimit.RedisClientMock, for: Towerops.RateLimit.RedisClient)
 
+# Tests assert on the bare ":no_tile" return path of the LIDAR catalog
+# without needing to mock the synthetic USGS NED runner. Real prod use
+# enables this fallback for nationwide CONUS coverage.
+Application.put_env(:towerops, :lidar_ned_fallback, false)
+
 # Define stub module for SNMP mock (returns errors for all calls)
 defmodule Towerops.Snmp.SnmpMockStub do
   @moduledoc false

test/towerops/lidar/grid_test.exs

@@ -90,6 +90,24 @@ defmodule Towerops.Lidar.GridTest do
       assert {:error, {:gdal, 1, "network unreachable"}} =
                Lidar.get_elevation_grid({-97.8, 30.2, -97.7, 30.3}, 0.05)
     end
+
+    test "falls back to synthetic USGS NED tiles when fallback is enabled" do
+      Application.put_env(:towerops, :lidar_ned_fallback, true)
+      on_exit(fn -> Application.put_env(:towerops, :lidar_ned_fallback, false) end)
+
+      Application.put_env(:towerops, :lidar_runner, fn _cmd, args, _opts ->
+        # Should be reading from the synthetic NED URL, not a DB tile.
+        joined = Enum.join(args, " ")
+        assert joined =~ "USGS_13_n31w098.tif"
+
+        {grid_aai([[150.0, 151.0], [152.0, 153.0]]), 0}
+      end)
+
+      assert {:ok, grid} =
+               Lidar.get_elevation_grid({-97.8, 30.2, -97.7, 30.3}, 0.05)
+
+      assert grid.cells == [[150.0, 151.0], [152.0, 153.0]]
+    end
   end
 
   defp insert_tile!(attrs \\ %{}) do

test/towerops/lidar/sources/usgs_ned_test.exs

@@ -0,0 +1,57 @@
+defmodule Towerops.Lidar.Sources.UsgsNedTest do
+  use ExUnit.Case, async: true
+
+  alias Towerops.Lidar.Sources.UsgsNed
+
+  describe "tile_for_point/2" do
+    test "returns a synthetic tile for a CONUS point" do
+      # Austin, TX
+      assert %{
+               source: "usgs_ned_13",
+               tile_name: "n31w098",
+               url: url
+             } = UsgsNed.tile_for_point(30.27, -97.74)
+
+      assert url =~ "USGS_13_n31w098.tif"
+      assert url =~ "/StagedProducts/Elevation/13/TIFF/current/"
+    end
+
+    test "returns a tile for the Hawaiian Islands" do
+      # O‘ahu
+      assert %{tile_name: "n22w158"} = UsgsNed.tile_for_point(21.3, -157.85)
+    end
+
+    test "returns a tile for Alaska" do
+      # Anchorage
+      assert %{tile_name: "n62w150"} = UsgsNed.tile_for_point(61.22, -149.9)
+    end
+
+    test "returns a tile for Puerto Rico" do
+      # San Juan
+      assert %{tile_name: "n19w067"} = UsgsNed.tile_for_point(18.45, -66.07)
+    end
+
+    test "returns nil outside catalogued NED extent" do
+      # mid-Atlantic open ocean
+      assert UsgsNed.tile_for_point(35.0, -45.0) == nil
+
+      # outside CONUS bbox (deep Mexico)
+      assert UsgsNed.tile_for_point(20.0, -100.0) == nil
+    end
+  end
+
+  describe "tiles_for_bbox/1" do
+    test "synthesizes one tile per intersected 1° NED cell" do
+      # 0.1° × 0.1° patch in Texas — straddles a single NED tile.
+      assert [tile] = UsgsNed.tiles_for_bbox({-97.8, 30.2, -97.7, 30.3})
+      assert tile.tile_name == "n31w098"
+    end
+
+    test "synthesizes multiple tiles when bbox crosses degree boundaries" do
+      # 1.5° × 1.5° patch crossing a degree boundary.
+      tiles = UsgsNed.tiles_for_bbox({-98.5, 30.5, -97.0, 32.0})
+      assert length(tiles) >= 4
+      assert Enum.all?(tiles, &(&1.source == "usgs_ned_13"))
+    end
+  end
+end