feat(buildings): integrate MS footprints into path clearance + map render

Phase 2/3/5 of the building-blockage support:

* Parser streams csv.gz tiles into compact records (centroid, radius,
  height) — keeps RAM ~32 B per polygon, drops -1.0-height entries.
* Index buckets records into 0.01° (~1 km) grid cells in a public ETS
  table for concurrent reads; max_height_near/3 and records_near/3
  scan only the 9 nearest buckets.
* Loader lazily parses any cached quadkey before a Calculate so the
  scorer is terrain-only when tiles aren't on disk yet.
* Rover.PathTerrain now treats each path-sample obstacle as
  terrain_elev + max_local_building_height, so blocked LOS through
  recent construction actually reduces clearance.
* Selecting a candidate pushes a candidate_buildings event with the
  buildings near each rover→station path; the hook renders them as
  height-coloured circles (yellow <10 m, orange 10-30 m, red >=30 m)
  with a tooltip showing height in meters.
This commit is contained in:
Graham McIntire 2026-04-26 10:44:45 -05:00
parent a150e02cd0
commit 0f15e99fb9
No known key found for this signature in database
GPG key ID: F4ABF488E6029E59
10 changed files with 467 additions and 5 deletions

View file

@ -28,6 +28,7 @@ interface RoverMapHook extends ViewHook {
qualityLayer: L.GridLayer | null
selectedCandidateMarker: L.CircleMarker | null
candidatePathLayer: L.LayerGroup
candidateBuildingsLayer: L.LayerGroup
gridLayer: L.LayerGroup
gridVisible: boolean
cellLookup: Map<string, Cell>
@ -90,6 +91,7 @@ export const RoverMap: Partial<RoverMapHook> = {
this.qualityLayer = null
this.selectedCandidateMarker = null
this.candidatePathLayer = L.layerGroup()
this.candidateBuildingsLayer = L.layerGroup()
this.gridLayer = L.layerGroup()
this.gridVisible = true
this.cellLookup = new Map()
@ -183,6 +185,7 @@ export const RoverMap: Partial<RoverMapHook> = {
map.fitBounds(this.driveCircle.getBounds(), { padding: [20, 20] })
this.candidatePathLayer.addTo(map)
this.candidateBuildingsLayer.addTo(map)
this.gridLayer.addTo(map)
updateGridOverlay(map, this.gridLayer)
@ -294,6 +297,30 @@ export const RoverMap: Partial<RoverMapHook> = {
).addTo(this.candidatePathLayer)
}
})
this.handleEvent("candidate_buildings", (
{ buildings }: {
buildings: { lat: number; lon: number; radius_m: number; height_m: number }[]
}
) => {
this.candidateBuildingsLayer.clearLayers()
if (!buildings || buildings.length === 0) return
for (const b of buildings) {
// Color by height: tall = red, medium = orange, low = yellow.
const color = b.height_m >= 30 ? "#dc2626" : b.height_m >= 10 ? "#f97316" : "#facc15"
L.circle([b.lat, b.lon], {
radius: Math.max(b.radius_m, 8),
color,
weight: 1,
fillColor: color,
fillOpacity: 0.45,
interactive: true
})
.bindTooltip(`${Math.round(b.height_m)} m`, { direction: "top", offset: [0, -4] })
.addTo(this.candidateBuildingsLayer)
}
})
},
reconnected(this: RoverMapHook) {

View file

@ -24,6 +24,7 @@ defmodule Microwaveprop.Application do
{PartitionSupervisor,
child_spec: Task.Supervisor, name: Microwaveprop.TaskSupervisor, partitions: System.schedulers_online()},
Microwaveprop.Cache,
Microwaveprop.Buildings.Index,
Microwaveprop.Propagation.ScoreCache,
Microwaveprop.Weather.GridCache,
{Microwaveprop.Weather.IemRateLimiter,

View file

@ -0,0 +1,164 @@
defmodule Microwaveprop.Buildings.Index do
@moduledoc """
In-memory spatial index over Microsoft building footprint records.
Buckets buildings into a ~1.1 km grid so that a max-height query in
a small radius scans only the 9 nearest buckets.
Records (`Microwaveprop.Buildings.Parser.record`) are bucketed by
centroid; queries widen by the building's `max_radius_m` so a long
warehouse straddling two buckets is still found by points in either.
Loading is one-shot per tile (idempotent via the `:loaded_tiles`
bucket). Queries hit ETS directly with no GenServer involvement.
"""
use GenServer
alias Microwaveprop.Buildings.Parser
@table :buildings_spatial_index
@loaded_tiles_table :buildings_loaded_tiles
# 0.01° ≈ 1.1 km at our latitude — strikes a balance between
# bucket scan size and number of empty buckets allocated.
@bucket_step 0.01
@spec start_link(keyword()) :: GenServer.on_start()
def start_link(opts), do: GenServer.start_link(__MODULE__, opts, name: __MODULE__)
@impl true
def init(_opts) do
ensure_tables()
{:ok, %{}}
end
@doc """
Returns the list of building records within `radius_m` of `(lat, lon)`.
Same scan path as `max_height_near/3` useful for surfacing the
individual obstacles to a UI layer.
"""
@spec records_near(float(), float(), number()) :: [Parser.record()]
def records_near(lat, lon, radius_m) do
ensure_tables()
radius_deg = radius_m / 111_000.0
cos_lat = :math.cos(lat * :math.pi() / 180.0)
{min_bx, max_bx} =
bucket_range(lon - radius_deg / max(cos_lat, 0.1), lon + radius_deg / max(cos_lat, 0.1))
{min_by, max_by} = bucket_range(lat - radius_deg, lat + radius_deg)
for by <- min_by..max_by,
bx <- min_bx..max_bx,
{_, list} <- :ets.lookup(@table, {bx, by}),
rec <- list,
within?(rec, lat, lon, cos_lat, radius_m) do
rec
end
end
defp within?(rec, lat, lon, cos_lat, radius_m) do
dlat_m = (rec.centroid_lat - lat) * 111_000.0
dlon_m = (rec.centroid_lon - lon) * 111_000.0 * cos_lat
dist_m = :math.sqrt(dlat_m * dlat_m + dlon_m * dlon_m)
dist_m - rec.max_radius_m <= radius_m
end
@doc """
Returns the maximum building `height_m` within `radius_m` of
`(lat, lon)`. `0.0` if no buildings indexed within the radius.
"""
@spec max_height_near(float(), float(), number()) :: float()
def max_height_near(lat, lon, radius_m) do
ensure_tables()
radius_deg = radius_m / 111_000.0
cos_lat = :math.cos(lat * :math.pi() / 180.0)
{min_bx, max_bx} = bucket_range(lon - radius_deg / max(cos_lat, 0.1), lon + radius_deg / max(cos_lat, 0.1))
{min_by, max_by} = bucket_range(lat - radius_deg, lat + radius_deg)
for by <- min_by..max_by,
bx <- min_bx..max_bx,
{_, list} <- :ets.lookup(@table, {bx, by}),
reduce: 0.0 do
acc -> max(acc, scan_bucket(list, lat, lon, cos_lat, radius_m))
end
end
defp scan_bucket(list, lat, lon, cos_lat, radius_m) do
Enum.reduce(list, 0.0, fn rec, acc ->
dlat_m = (rec.centroid_lat - lat) * 111_000.0
dlon_m = (rec.centroid_lon - lon) * 111_000.0 * cos_lat
dist_m = :math.sqrt(dlat_m * dlat_m + dlon_m * dlon_m)
if dist_m - rec.max_radius_m <= radius_m,
do: max(acc, rec.height_m),
else: acc
end)
end
@doc """
Insert a list of `Parser.record/0` maps into the spatial index.
Idempotent at the tile level (callers gate via `loaded?/1`); within
a tile, inserting the same record twice will produce a duplicate
entry but will not affect query correctness.
"""
@spec insert_records([Parser.record()]) :: :ok
def insert_records(records) when is_list(records) do
ensure_tables()
records
|> Enum.group_by(fn r -> {bucket_x(r.centroid_lon), bucket_y(r.centroid_lat)} end)
|> Enum.each(fn {key, recs} ->
existing =
case :ets.lookup(@table, key) do
[{_, list}] -> list
[] -> []
end
:ets.insert(@table, {key, recs ++ existing})
end)
:ok
end
@doc "Mark a quadkey as loaded so callers can skip re-parsing its tile."
@spec mark_loaded(String.t()) :: :ok
def mark_loaded(quadkey) do
ensure_tables()
:ets.insert(@loaded_tiles_table, {quadkey, true})
:ok
end
@spec loaded?(String.t()) :: boolean()
def loaded?(quadkey) do
ensure_tables()
:ets.member(@loaded_tiles_table, quadkey)
end
@doc "Wipe all index state (test helper / dev convenience)."
@spec clear() :: :ok
def clear do
ensure_tables()
:ets.delete_all_objects(@table)
:ets.delete_all_objects(@loaded_tiles_table)
:ok
end
defp ensure_tables do
if :ets.whereis(@table) == :undefined do
:ets.new(@table, [:set, :named_table, :public, read_concurrency: true])
end
if :ets.whereis(@loaded_tiles_table) == :undefined do
:ets.new(@loaded_tiles_table, [:set, :named_table, :public, read_concurrency: true])
end
end
defp bucket_range(min, max) do
{bucket_index(min), bucket_index(max)}
end
defp bucket_x(lon), do: bucket_index(lon)
defp bucket_y(lat), do: bucket_index(lat)
defp bucket_index(v), do: trunc(Float.floor(v / @bucket_step))
end

View file

@ -0,0 +1,60 @@
defmodule Microwaveprop.Buildings.Loader do
@moduledoc """
Lazily loads cached Microsoft footprint csv.gz tiles into the
spatial `Microwaveprop.Buildings.Index`. Idempotent: a quadkey
loaded once stays loaded for the lifetime of the BEAM.
Tiles that haven't been pre-fetched are skipped with a debug log —
the rover scorer falls back to terrain-only path analysis when the
index is missing for a region.
"""
alias Microwaveprop.Buildings.Index
alias Microwaveprop.Buildings.MsFootprints
alias Microwaveprop.Buildings.Parser
require Logger
@zoom 9
@doc """
Ensures every quadkey covering `bbox` has its csv.gz tile parsed
into the index. Returns the list of quadkeys that were *attempted*
(whether already loaded, freshly loaded, or missing on disk).
"""
@spec ensure_loaded_for_bbox(map()) :: [String.t()]
def ensure_loaded_for_bbox(bbox) do
quadkeys = MsFootprints.quadkeys_for_bbox(bbox, @zoom)
Enum.each(quadkeys, &ensure_loaded/1)
quadkeys
end
@doc """
Load a single quadkey. No-op if already loaded; logs and skips if
the csv.gz file is missing on disk.
"""
@spec ensure_loaded(String.t()) :: :ok
def ensure_loaded(quadkey) do
if Index.loaded?(quadkey) do
:ok
else
do_load(quadkey)
end
end
defp do_load(quadkey) do
path = Path.join(MsFootprints.cache_dir(), "#{quadkey}.csv.gz")
if File.exists?(path) do
records = path |> Parser.parse_tile() |> Enum.to_list()
Index.insert_records(records)
Index.mark_loaded(quadkey)
Logger.info("buildings index: loaded #{quadkey} (#{length(records)} polygons)")
else
Logger.debug("buildings index: tile #{quadkey} not on disk; skipping")
Index.mark_loaded(quadkey)
end
:ok
end
end

View file

@ -0,0 +1,79 @@
defmodule Microwaveprop.Buildings.Parser do
@moduledoc """
Streaming parser for Microsoft Building Footprint csv.gz tiles.
Each line is a GeoJSON `Feature` whose `geometry.coordinates` is a
Polygon (or MultiPolygon) and whose `properties.height` is meters
above ground. We collapse each polygon to a centroid + bounding
radius that's all the rover path-clearance scorer needs, and it
keeps the in-memory footprint to ~32 bytes per building.
Polygons with `height = -1.0` (no estimate from the ML model) are
dropped including them would just add noise.
"""
@type record :: %{
centroid_lat: float(),
centroid_lon: float(),
max_radius_m: float(),
height_m: float()
}
@doc """
Returns a `Stream` of `t:record/0` for every polygon in `path` whose
height is known. Lazy; suitable for tiles with millions of rows.
"""
@spec parse_tile(String.t()) :: Enumerable.t()
def parse_tile(path) do
path
|> File.stream!([:read, :compressed])
|> Stream.flat_map(&parse_line/1)
end
defp parse_line(line) do
case Jason.decode(String.trim(line)) do
{:ok, %{"properties" => %{"height" => h}, "geometry" => geom}} when h > 0.0 ->
case ring_from_geometry(geom) do
[] -> []
ring -> [build_record(ring, h)]
end
_ ->
[]
end
end
defp ring_from_geometry(%{"type" => "Polygon", "coordinates" => [outer | _]}), do: outer
defp ring_from_geometry(%{"type" => "MultiPolygon", "coordinates" => [[outer | _] | _]}), do: outer
defp ring_from_geometry(_), do: []
defp build_record(ring, height) do
{sum_lat, sum_lon, n, min_lat, min_lon, max_lat, max_lon} =
Enum.reduce(ring, {0.0, 0.0, 0, 90.0, 180.0, -90.0, -180.0}, fn
[lon, lat], {sl, slon, n, mnla, mnlo, mxla, mxlo} ->
{
sl + lat,
slon + lon,
n + 1,
min(mnla, lat),
min(mnlo, lon),
max(mxla, lat),
max(mxlo, lon)
}
end)
centroid_lat = sum_lat / n
centroid_lon = sum_lon / n
# Approximate max radius: half-diagonal of the bbox in meters.
dlat_m = (max_lat - min_lat) * 111_000.0 / 2.0
dlon_m = (max_lon - min_lon) * 111_000.0 * :math.cos(centroid_lat * :math.pi() / 180.0) / 2.0
%{
centroid_lat: centroid_lat,
centroid_lon: centroid_lon,
max_radius_m: :math.sqrt(dlat_m * dlat_m + dlon_m * dlon_m),
height_m: height * 1.0
}
end
end

View file

@ -7,6 +7,7 @@ defmodule Microwaveprop.Rover.Compute do
plus the top 5 candidate parking spots.
"""
alias Microwaveprop.Buildings.Loader, as: BuildingsLoader
alias Microwaveprop.Propagation
alias Microwaveprop.Radio.Maidenhead
alias Microwaveprop.Rover.Aggregator
@ -103,6 +104,7 @@ defmodule Microwaveprop.Rover.Compute do
points = Enum.map(in_radius, &{&1.lat, &1.lon})
elev_map = time_step("elev_lookup", fn -> elev_lookup.(points) end)
_ = time_step("buildings_load", fn -> BuildingsLoader.ensure_loaded_for_bbox(bbox) end)
clearance_map = time_step("clearance", fn -> clearance_lookup.(in_radius, selected_stations) end)
prominence_map = time_step("prominence", fn -> prominence_lookup.(in_radius) end)

View file

@ -9,11 +9,18 @@ defmodule Microwaveprop.Rover.PathTerrain do
clearance means the path is blocked by an intervening ridge.
"""
alias Microwaveprop.Buildings.Index, as: BuildingsIndex
alias Microwaveprop.Rover.Elevation
# Number of intermediate samples between rover and station.
@sample_count 12
# Search radius around each path sample point when looking for
# buildings that obstruct the line of sight. ~150 m catches the
# typical width of a downtown high-rise plus a margin for building
# centroid imprecision.
@building_search_radius_m 150
@type latlon :: {float(), float()}
@doc """
@ -25,6 +32,7 @@ defmodule Microwaveprop.Rover.PathTerrain do
@spec clearance_map([map()], [map()], keyword()) :: %{{latlon(), latlon()} => integer() | nil}
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)
pairs =
for cell <- cells, station <- stations do
@ -40,11 +48,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)}
{{rover_pt, station_pt}, clearance(rover_pt, station_pt, elev, buildings_lookup)}
end)
end
defp clearance(rover_pt, station_pt, elev) do
defp clearance(rover_pt, station_pt, elev, buildings_lookup) do
rover_elev = Map.get(elev, rover_pt)
case rover_elev do
@ -55,7 +63,7 @@ defmodule Microwaveprop.Rover.PathTerrain do
path_max =
rover_pt
|> intermediate_samples(station_pt)
|> Enum.map(&Map.get(elev, &1))
|> Enum.map(&obstacle_top(&1, elev, buildings_lookup))
|> Enum.reject(&is_nil/1)
|> safe_max()
@ -63,6 +71,20 @@ defmodule Microwaveprop.Rover.PathTerrain do
end
end
defp obstacle_top({lat, lon} = pt, elev, buildings_lookup) do
case Map.get(elev, pt) do
nil ->
nil
ground ->
ground + round(buildings_lookup.({lat, lon}))
end
end
defp default_building_height({lat, lon}) do
BuildingsIndex.max_height_near(lat, lon, @building_search_radius_m)
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).

View file

@ -9,6 +9,7 @@ defmodule MicrowavepropWeb.RoverLive do
use MicrowavepropWeb, :live_view
alias Microwaveprop.Accounts.User
alias Microwaveprop.Buildings.Index, as: BuildingsIndex
alias Microwaveprop.Geocoder
alias Microwaveprop.Propagation
alias Microwaveprop.Radio.CallsignClient
@ -287,6 +288,7 @@ defmodule MicrowavepropWeb.RoverLive do
candidate ->
detail = build_candidate_detail(candidate, socket)
paths = candidate_path_payload(detail, socket.assigns.fixed_stations)
{:noreply,
socket
@ -294,7 +296,10 @@ defmodule MicrowavepropWeb.RoverLive do
|> push_event("focus_cell", %{lat: candidate.lat, lon: candidate.lon, zoom: 11})
|> push_event("candidate_paths", %{
candidate: %{lat: candidate.lat, lon: candidate.lon},
paths: candidate_path_payload(detail, socket.assigns.fixed_stations)
paths: paths
})
|> push_event("candidate_buildings", %{
buildings: buildings_along_paths(candidate, paths)
})}
end
end
@ -303,7 +308,8 @@ defmodule MicrowavepropWeb.RoverLive do
{:noreply,
socket
|> assign(selected_candidate: nil)
|> push_event("candidate_paths", %{candidate: nil, paths: []})}
|> push_event("candidate_paths", %{candidate: nil, paths: []})
|> push_event("candidate_buildings", %{buildings: []})}
end
def handle_event("rover_cell_detail", %{"lat" => lat, "lon" => lon}, socket) do
@ -733,6 +739,34 @@ defmodule MicrowavepropWeb.RoverLive do
Map.put(link, :profile_svg, profile_to_svg(link.profile))
end
defp buildings_along_paths(candidate, paths) do
paths
|> Enum.flat_map(fn p ->
sample_along({candidate.lat, candidate.lon}, {p.station_lat, p.station_lon})
end)
|> Enum.flat_map(fn {lat, lon} ->
BuildingsIndex.records_near(lat, lon, 150)
end)
|> Enum.uniq_by(fn r -> {r.centroid_lat, r.centroid_lon} end)
|> Enum.map(fn r ->
%{
lat: r.centroid_lat,
lon: r.centroid_lon,
radius_m: r.max_radius_m,
height_m: r.height_m
}
end)
end
defp sample_along({lat1, lon1}, {lat2, lon2}) do
n = 24
for i <- 1..n do
f = i / (n + 1)
{lat1 + f * (lat2 - lat1), lon1 + f * (lon2 - lon1)}
end
end
defp candidate_path_payload(detail, fixed_stations) do
by_callsign = Map.new(fixed_stations, fn s -> {s.callsign, s} end)

View file

@ -0,0 +1,33 @@
defmodule Microwaveprop.Buildings.IndexTest do
use ExUnit.Case, async: false
alias Microwaveprop.Buildings.Index
setup do
Index.clear()
:ok
end
test "max_height_near/3 returns 0 when no buildings have been indexed" do
assert Index.max_height_near(32.8, -97.0, 100) == 0.0
end
test "insert_records/1 then max_height_near/3 returns the tallest nearby height" do
Index.insert_records([
%{centroid_lat: 32.8000, centroid_lon: -97.0000, max_radius_m: 10.0, height_m: 8.0},
%{centroid_lat: 32.8001, centroid_lon: -97.0001, max_radius_m: 10.0, height_m: 25.0},
# Far away — should not be considered.
%{centroid_lat: 33.0, centroid_lon: -97.5, max_radius_m: 10.0, height_m: 100.0}
])
assert Index.max_height_near(32.8000, -97.0000, 50) == 25.0
end
test "max_height_near/3 returns 0 when buildings are outside the radius" do
Index.insert_records([
%{centroid_lat: 32.8500, centroid_lon: -97.0500, max_radius_m: 10.0, height_m: 50.0}
])
assert Index.max_height_near(32.8000, -97.0000, 100) == 0.0
end
end

View file

@ -0,0 +1,40 @@
defmodule Microwaveprop.Buildings.ParserTest do
use ExUnit.Case, async: true
alias Microwaveprop.Buildings.Parser
@sample_geojsonl """
{"type": "Feature", "properties": {"height": 5.5, "confidence": -1.0}, "geometry": {"type": "Polygon", "coordinates": [[[-97.736, 32.221], [-97.736, 32.222], [-97.735, 32.222], [-97.735, 32.221], [-97.736, 32.221]]]}}
{"type": "Feature", "properties": {"height": 12.0, "confidence": 0.9}, "geometry": {"type": "Polygon", "coordinates": [[[-97.0, 32.5], [-97.0, 32.501], [-96.999, 32.501], [-96.999, 32.5], [-97.0, 32.5]]]}}
{"type": "Feature", "properties": {"height": -1.0, "confidence": -1.0}, "geometry": {"type": "Polygon", "coordinates": [[[-97.5, 32.0], [-97.5, 32.001], [-97.499, 32.001], [-97.499, 32.0], [-97.5, 32.0]]]}}
"""
setup do
path = Path.join(System.tmp_dir!(), "buildings_parser_#{:erlang.unique_integer([:positive])}.csv.gz")
gz = :zlib.gzip(@sample_geojsonl)
File.write!(path, gz)
on_exit(fn -> File.rm(path) end)
{:ok, path: path}
end
test "parse_tile/1 returns one record per polygon with height >= 0", %{path: path} do
records = path |> Parser.parse_tile() |> Enum.to_list()
assert length(records) == 2
assert Enum.all?(records, fn r -> r.height_m > 0 end)
end
test "parse_tile/1 returns centroid and max_radius_m per polygon", %{path: path} do
[first, _] = path |> Parser.parse_tile() |> Enum.to_list()
assert_in_delta first.centroid_lat, 32.2215, 0.001
assert_in_delta first.centroid_lon, -97.7355, 0.001
assert first.max_radius_m > 0
assert first.height_m == 5.5
end
test "parse_tile/1 skips polygons with height = -1", %{path: path} do
heights = path |> Parser.parse_tile() |> Enum.map(& &1.height_m)
refute -1.0 in heights
end
end