From f77c493fa34e4917286cca267288cd6a425ed2e5 Mon Sep 17 00:00:00 2001 From: Graham McIntire Date: Wed, 6 May 2026 16:21:50 -0500 Subject: [PATCH] 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 | 64 ++++++++++ lib/towerops/lidar.ex | 38 +++++- lib/towerops/lidar/sources/usgs_ned.ex | 111 ++++++++++++++++++ lib/towerops/workers/coverage_worker.ex | 5 +- ...260506205000_create_coverage_buildings.exs | 40 +++++++ test/test_helper.exs | 5 + test/towerops/lidar/grid_test.exs | 18 +++ test/towerops/lidar/sources/usgs_ned_test.exs | 57 +++++++++ 8 files changed, 334 insertions(+), 4 deletions(-) create mode 100644 lib/towerops/coverages/building.ex create mode 100644 lib/towerops/lidar/sources/usgs_ned.ex create mode 100644 priv/repo/migrations/20260506205000_create_coverage_buildings.exs create mode 100644 test/towerops/lidar/sources/usgs_ned_test.exs diff --git a/lib/towerops/coverages/building.ex b/lib/towerops/coverages/building.ex new file mode 100644 index 00000000..f5504d8f --- /dev/null +++ b/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 diff --git a/lib/towerops/lidar.ex b/lib/towerops/lidar.ex index c5c613a8..2725eeb1 100644 --- a/lib/towerops/lidar.ex +++ b/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} -> diff --git a/lib/towerops/lidar/sources/usgs_ned.ex b/lib/towerops/lidar/sources/usgs_ned.ex new file mode 100644 index 00000000..d4990346 --- /dev/null +++ b/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 diff --git a/lib/towerops/workers/coverage_worker.ex b/lib/towerops/workers/coverage_worker.ex index 64f400fa..c9ba4ef3 100644 --- a/lib/towerops/workers/coverage_worker.ex +++ b/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." diff --git a/priv/repo/migrations/20260506205000_create_coverage_buildings.exs b/priv/repo/migrations/20260506205000_create_coverage_buildings.exs new file mode 100644 index 00000000..e8d56cf9 --- /dev/null +++ b/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 diff --git a/test/test_helper.exs b/test/test_helper.exs index c91dc941..3e3011d0 100644 --- a/test/test_helper.exs +++ b/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 diff --git a/test/towerops/lidar/grid_test.exs b/test/towerops/lidar/grid_test.exs index 1a124f52..ae9a4cbf 100644 --- a/test/towerops/lidar/grid_test.exs +++ b/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 diff --git a/test/towerops/lidar/sources/usgs_ned_test.exs b/test/towerops/lidar/sources/usgs_ned_test.exs new file mode 100644 index 00000000..97c25a87 --- /dev/null +++ b/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