feat(rainscatter): classify QSO propagation mechanism from common-volume radar

Adds a per-contact enrichment pipeline that determines whether a QSO was
most likely carried by rain scatter, tropospheric ducting, or ordinary
troposcatter — using IEM n0q composite reflectivity sampled inside the
lens-shaped intersection of 400 km-radius disks around each endpoint.

Pieces:
  * Microwaveprop.Propagation.CommonVolume — lens geometry (haversine,
    in-CV test, bbox, area).
  * contact_common_volume_radar table (1:1 per contact) storing
    aggregate dBZ stats inside the CV + radar_status column on contacts.
  * Microwaveprop.Workers.CommonVolumeRadarWorker — Oban :radar queue,
    fetches the n0q frame at QSO time, iterates pixels inside the CV
    bbox, aggregates rain/heavy/core-pixel counts, max/mean dBZ, and
    coverage percentage.
  * Microwaveprop.Propagation.RainScatterClassifier — rule-based mapper
    from (band, distance, radar stats, duct flags) to one of
    :likely_rainscatter | :rainscatter_possible | :tropo_duct |
    :troposcatter | :unknown.
  * ContactWeatherEnqueueWorker learns a :radar enrichment type and
    enqueues the CV worker on contact submission; pre-2014 contacts
    (outside IEM n0q coverage) are pinned to :unavailable.
  * `mix radar_backfill` bulk-enqueues historical contacts with
    --year / --limit / --dry-run.
  * Contact detail page renders a mechanism badge with supporting
    stats (common-volume area, max dBZ, heavy-rain pixel count,
    coverage %).
This commit is contained in:
Graham McIntire 2026-04-17 15:57:54 -05:00
parent c09446a517
commit 33f5d4edbe
No known key found for this signature in database
GPG key ID: F4ABF488E6029E59
19 changed files with 1107 additions and 3 deletions

View file

@ -68,6 +68,7 @@ config :microwaveprop, Oban,
propagation: 2,
admin: 1,
nexrad: 2,
radar: 2,
ionosphere: 1,
space_weather: 1,
contact_import: 4

View file

@ -83,6 +83,7 @@ config :microwaveprop, Oban,
commercial: 2,
iemre: 10,
narr: 6,
radar: 2,
backfill_enqueue: 1,
exports: 1,
contact_import: 4

View file

@ -179,6 +179,7 @@ if config_env() == :prod do
backfill_enqueue: 1,
admin: 1,
nexrad: 2,
radar: 2,
ionosphere: 1,
space_weather: 1,
exports: 1,

View file

@ -0,0 +1,97 @@
defmodule Microwaveprop.Propagation.CommonVolume do
@moduledoc """
Geometry for the "common volume" of a microwave QSO the lens-shaped
intersection of two circles of radius R around the endpoints.
Rain-scatter requires both stations' beams to intersect a precipitating
cell, so the cell has to live inside both 400 km-radius neighborhoods.
This module provides the geometric primitives a classifier needs:
* `in_common_volume?/4` is a point inside both disks?
* `bounding_box/3` lat/lon rectangle that contains the lens
* `area_km2/3` area of the intersection, for coverage normalization
Distances are spherical-great-circle (haversine) at mean Earth radius.
"""
@earth_radius_km 6371.0
@km_per_deg_lat 111.32
@type latlon :: {float(), float()}
@type bbox :: %{min_lat: float(), max_lat: float(), min_lon: float(), max_lon: float()}
@doc "Is `point` within `radius_km` of both endpoints?"
@spec in_common_volume?(latlon(), latlon(), latlon(), float()) :: boolean()
def in_common_volume?(pos1, pos2, point, radius_km) do
haversine_km(pos1, point) <= radius_km and haversine_km(pos2, point) <= radius_km
end
@doc """
Lat/lon rectangle containing the common volume. Returns `:empty` when
the two circles don't overlap.
"""
@spec bounding_box(latlon(), latlon(), float()) :: bbox() | :empty
def bounding_box({lat1, lon1} = pos1, {lat2, lon2} = pos2, radius_km) do
if haversine_km(pos1, pos2) > 2 * radius_km do
:empty
else
lat_pad = radius_km / @km_per_deg_lat
# Approximate longitude degrees-per-km at the higher-latitude endpoint
# (narrower at higher lat, so more conservative).
ref_lat = max(abs(lat1), abs(lat2))
km_per_deg_lon = @km_per_deg_lat * max(:math.cos(ref_lat * :math.pi() / 180.0), 0.01)
lon_pad = radius_km / km_per_deg_lon
%{
min_lat: max(lat1 - lat_pad, lat2 - lat_pad),
max_lat: min(lat1 + lat_pad, lat2 + lat_pad),
min_lon: max(lon1 - lon_pad, lon2 - lon_pad),
max_lon: min(lon1 + lon_pad, lon2 + lon_pad)
}
end
end
@doc """
Area (km²) of the lens-shaped intersection of two circles of radius
`radius_km` centered at `pos1` and `pos2`. Returns `0.0` when the
circles don't overlap. Uses the standard two-circle lens formula
with the great-circle distance between centers.
"""
@spec area_km2(latlon(), latlon(), float()) :: float()
def area_km2(pos1, pos2, radius_km) do
d = haversine_km(pos1, pos2)
cond do
d >= 2 * radius_km ->
0.0
d <= 0.0 ->
:math.pi() * radius_km * radius_km
true ->
r = radius_km
# Lens area = 2 * [r² * arccos(d/2r) - (d/4) * sqrt(4r² - d²)]
part_arc = r * r * :math.acos(d / (2 * r))
part_tri = d / 4 * :math.sqrt(4 * r * r - d * d)
2 * (part_arc - part_tri)
end
end
@doc "Great-circle distance (km) between two lat/lon points."
@spec haversine_km(latlon(), latlon()) :: float()
def haversine_km({lat1, lon1}, {lat2, lon2}) do
rlat1 = lat1 * :math.pi() / 180.0
rlat2 = lat2 * :math.pi() / 180.0
dlat = (lat2 - lat1) * :math.pi() / 180.0
dlon = (lon2 - lon1) * :math.pi() / 180.0
a =
:math.sin(dlat / 2) * :math.sin(dlat / 2) +
:math.cos(rlat1) * :math.cos(rlat2) *
:math.sin(dlon / 2) * :math.sin(dlon / 2)
c = 2.0 * :math.atan2(:math.sqrt(a), :math.sqrt(1.0 - a))
@earth_radius_km * c
end
end

View file

@ -0,0 +1,109 @@
defmodule Microwaveprop.Propagation.RainScatterClassifier do
@moduledoc """
Rule-based classifier that labels a contact with its most likely
propagation mechanism: rain scatter vs tropospheric ducting vs plain
troposcatter vs unknown.
Takes a small input map so callers can assemble it from whatever data
they have (schema rows at submit time, aggregated SQL in backfill,
etc.) without this module depending on Ecto.
## Inputs
%{
band_mhz: integer(),
distance_km: float(),
# common-volume radar stats (nil if no radar sweep covers the CV)
radar: nil | %{
max_dbz: float() | nil,
heavy_rain_pixel_count: non_neg_integer(),
coverage_pct: float() | nil
},
# true if HRRR detected a trapping layer at *either* endpoint
duct_either_endpoint: boolean()
}
## Outputs
* `:likely_rainscatter` 5-11 GHz, path 800 km, heavy rain ( 35 dBZ
with several pixels) inside the common volume, no duct at either end.
* `:rainscatter_possible` light rain (25-35 dBZ) in the CV, no duct;
could be rainscatter, could be weak tropo.
* `:tropo_duct` HRRR ducting signature at either endpoint dominates.
* `:troposcatter` no duct, no meaningful rain default tropospheric
forward scatter, which carries most "clear-air" microwave contacts.
* `:unknown` not enough information (no radar sweep, poor coverage).
These labels are coarse on purpose: with post-hoc HRRR + n0q data you
can't prove the mechanism, but you *can* separate the cohorts enough
to calibrate band-weights and to show the likely mode on a QSO detail
page.
"""
# Rainscatter is feasible roughly 5-11 GHz. Below, scattering cross-section
# is too small; above, absorption kills the signal before it scatters.
@rs_band_min_mhz 5_000
@rs_band_max_mhz 11_000
# Geometry limit: common volume stops existing beyond ~2R = 800 km.
@rs_max_distance_km 800.0
# Reflectivity thresholds (dBZ) on the n0q composite.
@light_rain_dbz 25.0
@heavy_rain_dbz 35.0
# How many heavy pixels count as a "real" thunderstorm cell in the CV.
@heavy_pixel_floor 3
# Minimum radar coverage pct to trust the reading. Below this the CV
# scan didn't have enough valid pixels to say either way.
@min_coverage_pct 25.0
@type result ::
:likely_rainscatter
| :rainscatter_possible
| :tropo_duct
| :troposcatter
| :unknown
@spec classify(map()) :: result()
def classify(%{duct_either_endpoint: true}), do: :tropo_duct
def classify(%{band_mhz: band, distance_km: dist} = input) do
if rainscatter_geometry?(band, dist) do
classify_from_radar(input)
else
fallback_non_rainscatter(input)
end
end
defp rainscatter_geometry?(band_mhz, distance_km) do
band_mhz >= @rs_band_min_mhz and band_mhz <= @rs_band_max_mhz and
distance_km <= @rs_max_distance_km
end
defp classify_from_radar(%{radar: nil}), do: :unknown
defp classify_from_radar(%{radar: radar}) do
coverage = radar[:coverage_pct] || 0.0
max_dbz = radar[:max_dbz] || 0.0
heavy = radar[:heavy_rain_pixel_count] || 0
cond do
coverage < @min_coverage_pct ->
:unknown
max_dbz >= @heavy_rain_dbz and heavy >= @heavy_pixel_floor ->
:likely_rainscatter
max_dbz >= @light_rain_dbz ->
:rainscatter_possible
true ->
:troposcatter
end
end
defp fallback_non_rainscatter(%{radar: nil}), do: :unknown
defp fallback_non_rainscatter(_), do: :troposcatter
end

View file

@ -300,6 +300,10 @@ defmodule Microwaveprop.Radio do
@spec unprocessed_iemre_contacts(non_neg_integer()) :: [Contact.t()]
def unprocessed_iemre_contacts(limit \\ 500), do: contacts_needing_enrichment(:iemre_status, limit)
@spec unprocessed_radar_contacts(non_neg_integer()) :: [Contact.t()]
def unprocessed_radar_contacts(limit \\ 500),
do: contacts_needing_enrichment(:radar_status, limit, &where(&1, [q], not is_nil(q.pos2)))
@spec set_enrichment_status!([Ecto.UUID.t()], atom(), atom()) :: {non_neg_integer(), nil | [any()]}
def set_enrichment_status!(ids, field, status) do
Contact

View file

@ -38,6 +38,10 @@ defmodule Microwaveprop.Radio.Contact do
values: [:pending, :queued, :processing, :complete, :failed, :unavailable],
default: :pending
field :radar_status, Ecto.Enum,
values: [:pending, :queued, :processing, :complete, :failed, :unavailable],
default: :pending
field :user_submitted, :boolean, default: false
field :submitter_email, :string
field :flagged_invalid, :boolean, default: false

View file

@ -0,0 +1,47 @@
defmodule Microwaveprop.Radio.ContactCommonVolumeRadar do
@moduledoc """
Aggregated n0q composite reflectivity statistics sampled inside the
common volume (lens-shaped intersection of two 400 km circles around
the endpoints) for a single contact, at the QSO time.
One row per contact. Feeds `Microwaveprop.Propagation.RainScatterClassifier`.
"""
use Ecto.Schema
import Ecto.Changeset
alias Microwaveprop.Radio.Contact
@primary_key {:id, :binary_id, autogenerate: true}
@foreign_key_type :binary_id
schema "contact_common_volume_radar" do
belongs_to :contact, Contact
field :observed_at, :utc_datetime
field :common_volume_km2, :float
field :pixel_count, :integer, default: 0
field :rain_pixel_count, :integer, default: 0
field :heavy_rain_pixel_count, :integer, default: 0
field :core_pixel_count, :integer, default: 0
field :max_dbz, :float
field :mean_dbz, :float
field :coverage_pct, :float
timestamps(type: :utc_datetime)
end
@type t :: %__MODULE__{}
@required ~w(contact_id observed_at)a
@optional ~w(common_volume_km2 pixel_count rain_pixel_count heavy_rain_pixel_count
core_pixel_count max_dbz mean_dbz coverage_pct)a
@spec changeset(t() | Ecto.Changeset.t(), map()) :: Ecto.Changeset.t()
def changeset(row, attrs) do
row
|> cast(attrs, @required ++ @optional)
|> validate_required(@required)
|> unique_constraint(:contact_id)
end
end

View file

@ -92,6 +92,19 @@ defmodule Microwaveprop.Weather.NexradClient do
end
end
@doc """
Fetch and decode the n0q PNG frame nearest the given timestamp.
Returns `{:ok, pixels, width}` a flat binary of palette-index bytes
or `{:error, reason}`. Does NOT go through `NexradCache`; callers that
want caching (e.g. real-time point queries) should call `fetch_rain_cells/4`.
"""
@spec fetch_decoded_frame(DateTime.t()) ::
{:ok, binary(), non_neg_integer()} | {:error, term()}
def fetch_decoded_frame(%DateTime{} = dt) do
rounded = round_to_5min(dt)
fetch_and_decode_frame(rounded)
end
defp fetch_and_decode_frame(rounded) do
url = frame_url(rounded)
req_opts = Application.get_env(:microwaveprop, :nexrad_req_options, [])

View file

@ -0,0 +1,256 @@
defmodule Microwaveprop.Workers.CommonVolumeRadarWorker do
@moduledoc """
Per-contact radar enrichment.
For each contact, fetches the IEM n0q composite-reflectivity frame
closest to the QSO timestamp and aggregates the pixels that fall
inside the lens-shaped common volume between the two endpoints (the
intersection of 400 km-radius disks). The aggregated statistics feed
the rain-scatter vs tropo classifier.
Jobs are unique on `contact_id` so the backfill and submission-time
enqueue paths collapse to a single job per contact.
"""
use Oban.Worker,
queue: :radar,
max_attempts: 3,
unique: [
period: :infinity,
states: [:available, :scheduled, :executing, :retryable],
keys: [:contact_id]
]
alias Microwaveprop.Propagation.CommonVolume
alias Microwaveprop.Radio.Contact
alias Microwaveprop.Radio.ContactCommonVolumeRadar
alias Microwaveprop.Repo
alias Microwaveprop.Weather.NexradClient
require Logger
# Radius of each station's effective rain-scatter neighborhood. Beyond
# this, the bistatic geometry stops being plausible given typical rain
# tops at ~15 km AGL.
@radius_km 400.0
# Pixel-to-dBZ thresholds (see NexradClient.pixel_to_dbz/1).
@rain_dbz 25.0
@heavy_rain_dbz 35.0
@core_dbz 40.0
# Sample every Nth pixel inside the common-volume bounding box. n0q is
# 0.005°/px ≈ 0.5 km/px, so step=10 samples at ~5 km resolution —
# enough to catch thunderstorm cores (typically 10+ km wide).
@pixel_step 10
@impl Oban.Worker
def perform(%Oban.Job{args: %{"contact_id" => contact_id}}) do
case Repo.get(Contact, contact_id) do
nil ->
:ok
%Contact{pos1: p1, pos2: p2} = contact when is_map(p1) and is_map(p2) ->
process(contact)
contact ->
mark_unavailable(contact)
:ok
end
end
defp process(%Contact{} = contact) do
pos1 = to_latlon(contact.pos1)
pos2 = to_latlon(contact.pos2)
case CommonVolume.bounding_box(pos1, pos2, @radius_km) do
:empty ->
mark_unavailable(contact)
:ok
_bbox ->
fetch_and_store(contact, pos1, pos2)
end
end
defp fetch_and_store(%Contact{} = contact, pos1, pos2) do
qso_at = contact.qso_timestamp
rounded = NexradClient.round_to_5min(DateTime.from_naive!(qso_at, "Etc/UTC"))
case NexradClient.fetch_decoded_frame(rounded) do
{:ok, pixels, width} ->
stats = aggregate_stats(pixels, width, pos1, pos2)
upsert_row!(contact, rounded, stats)
mark_status(contact, :complete)
:ok
{:error, reason} ->
Logger.info("CommonVolumeRadarWorker: no frame for #{contact.id}: #{inspect(reason)}")
mark_unavailable(contact)
:ok
end
end
@doc """
Aggregate n0q pixel statistics over the common volume of the two
endpoints. Pure function accepts raw pixel buffer + width so it
can be unit-tested without hitting the network.
`step` controls the pixel sampling stride (default 10 5 km/sample at
the n0q resolution). Tests pass `step: 1` to hit every pixel of a
small synthetic buffer.
"""
@spec aggregate_stats(
binary(),
non_neg_integer(),
CommonVolume.latlon(),
CommonVolume.latlon(),
keyword()
) ::
%{
pixel_count: non_neg_integer(),
rain_pixel_count: non_neg_integer(),
heavy_rain_pixel_count: non_neg_integer(),
core_pixel_count: non_neg_integer(),
max_dbz: float() | nil,
mean_dbz: float() | nil,
common_volume_km2: float(),
coverage_pct: float() | nil
}
def aggregate_stats(pixels, width, pos1, pos2, opts \\ []) do
step = Keyword.get(opts, :step, @pixel_step)
case CommonVolume.bounding_box(pos1, pos2, @radius_km) do
:empty ->
empty_stats(CommonVolume.area_km2(pos1, pos2, @radius_km))
%{min_lat: min_lat, max_lat: max_lat, min_lon: min_lon, max_lon: max_lon} ->
height = div(byte_size(pixels), width)
{x_min, y_max} = NexradClient.latlon_to_pixel(min_lat, min_lon)
{x_max, y_min} = NexradClient.latlon_to_pixel(max_lat, max_lon)
x_min = max(x_min, 0)
x_max = min(x_max, width - 1)
y_min = max(y_min, 0)
y_max = min(y_max, height - 1)
pixels
|> collect_cv_pixels(
width,
x_min..x_max,
y_min..y_max,
pos1,
pos2,
step
)
|> finalize_stats(CommonVolume.area_km2(pos1, pos2, @radius_km))
end
end
defp collect_cv_pixels(pixels, width, x_range, y_range, pos1, pos2, step) do
# We accumulate three counters (in_cv, echo_count, cumulative dBZ) plus
# max and threshold counts. Using a reduce keeps GC low on the ~10^5
# pixel scan per contact.
for y <- y_range.first..y_range.last//step,
x <- x_range.first..x_range.last//step,
reduce: init_acc() do
acc -> visit_pixel(acc, pixels, width, x, y, pos1, pos2)
end
end
defp visit_pixel(acc, pixels, width, x, y, pos1, pos2) do
lat = 50.0 - y * 0.005
lon = -126.0 + x * 0.005
if CommonVolume.in_common_volume?(pos1, pos2, {lat, lon}, @radius_km) do
offset = y * width + x
<<_::binary-size(offset), pixel_val::8, _::binary>> = pixels
update_acc(acc, pixel_val)
else
acc
end
end
defp init_acc do
%{
in_cv: 0,
echo_count: 0,
dbz_sum: 0.0,
max_dbz: nil,
rain: 0,
heavy: 0,
core: 0
}
end
defp update_acc(acc, 0), do: %{acc | in_cv: acc.in_cv + 1}
defp update_acc(acc, pixel_val) do
dbz = NexradClient.pixel_to_dbz(pixel_val)
%{
acc
| in_cv: acc.in_cv + 1,
echo_count: acc.echo_count + 1,
dbz_sum: acc.dbz_sum + dbz,
max_dbz: max(acc.max_dbz || dbz, dbz),
rain: acc.rain + if(dbz >= @rain_dbz, do: 1, else: 0),
heavy: acc.heavy + if(dbz >= @heavy_rain_dbz, do: 1, else: 0),
core: acc.core + if(dbz >= @core_dbz, do: 1, else: 0)
}
end
defp finalize_stats(%{in_cv: 0}, cv_area) do
empty_stats(cv_area)
end
defp finalize_stats(acc, cv_area) do
mean = if acc.echo_count > 0, do: acc.dbz_sum / acc.echo_count
%{
pixel_count: acc.in_cv,
rain_pixel_count: acc.rain,
heavy_rain_pixel_count: acc.heavy,
core_pixel_count: acc.core,
max_dbz: acc.max_dbz,
mean_dbz: mean,
common_volume_km2: cv_area,
coverage_pct: 100.0 * acc.echo_count / max(acc.in_cv, 1)
}
end
defp empty_stats(cv_area) do
%{
pixel_count: 0,
rain_pixel_count: 0,
heavy_rain_pixel_count: 0,
core_pixel_count: 0,
max_dbz: nil,
mean_dbz: nil,
common_volume_km2: cv_area,
coverage_pct: nil
}
end
defp upsert_row!(%Contact{id: contact_id}, observed_at, stats) do
%ContactCommonVolumeRadar{}
|> ContactCommonVolumeRadar.changeset(Map.merge(stats, %{contact_id: contact_id, observed_at: observed_at}))
|> Repo.insert(
on_conflict: {:replace_all_except, [:id, :contact_id, :inserted_at]},
conflict_target: :contact_id
)
end
defp mark_unavailable(%Contact{} = contact), do: mark_status(contact, :unavailable)
defp mark_status(%Contact{} = contact, status) do
contact
|> Ecto.Changeset.change(%{radar_status: status})
|> Repo.update!()
end
defp to_latlon(%{"lat" => lat, "lon" => lon}) when is_number(lat) and is_number(lon) do
{lat / 1.0, lon / 1.0}
end
end

View file

@ -6,6 +6,7 @@ defmodule Microwaveprop.Workers.ContactWeatherEnqueueWorker do
alias Microwaveprop.Weather
alias Microwaveprop.Weather.HrrrClient
alias Microwaveprop.Weather.NarrClient
alias Microwaveprop.Workers.CommonVolumeRadarWorker
alias Microwaveprop.Workers.HrrrFetchWorker
alias Microwaveprop.Workers.IemreFetchWorker
alias Microwaveprop.Workers.NarrFetchWorker
@ -22,12 +23,14 @@ defmodule Microwaveprop.Workers.ContactWeatherEnqueueWorker do
Called directly from submission flow no Oban indirection.
Optionally pass a list of types to limit which enrichments run.
"""
def enqueue_for_contact(contact, types \\ [:hrrr, :weather, :terrain, :iemre]) do
def enqueue_for_contact(contact, types \\ [:hrrr, :weather, :terrain, :iemre, :radar]) do
contact = Radio.ensure_positions!(contact)
jobs_by_type = build_jobs_by_type(contact, types)
bulk_jobs =
jobs_by_type[:weather] ++ jobs_by_type[:hrrr] ++ jobs_by_type[:terrain] ++ jobs_by_type[:iemre]
jobs_by_type[:weather] ++
jobs_by_type[:hrrr] ++
jobs_by_type[:terrain] ++ jobs_by_type[:iemre] ++ jobs_by_type[:radar]
if bulk_jobs != [] do
insert_all_chunked(bulk_jobs)
@ -48,7 +51,8 @@ defmodule Microwaveprop.Workers.ContactWeatherEnqueueWorker do
hrrr: if(:hrrr in types, do: build_hrrr_jobs([contact]), else: []),
terrain: if(:terrain in types, do: build_terrain_jobs([contact]), else: []),
iemre: if(:iemre in types, do: build_iemre_jobs([contact]), else: []),
narr: if(:narr in types, do: build_narr_jobs([contact]), else: [])
narr: if(:narr in types, do: build_narr_jobs([contact]), else: []),
radar: if(:radar in types, do: build_radar_jobs([contact]), else: [])
}
end
@ -63,9 +67,24 @@ defmodule Microwaveprop.Workers.ContactWeatherEnqueueWorker do
if contact.pos1 && contact.pos2 do
if :terrain in types, do: mark_status!(ids, :terrain_status, jobs_by_type[:terrain])
if :radar in types, do: mark_radar_status!(contact, ids, jobs_by_type[:radar])
end
end
# IEM n0q composite reflectivity archive only has meaningful CONUS
# coverage post-2014 (pre that, the archive is spotty mosaics). For
# pre-2014 contacts we can't get radar, so pin to :unavailable rather
# than leaving it :pending forever.
defp mark_radar_status!(contact, ids, []) do
if NarrClient.in_coverage?(contact.qso_timestamp) do
Radio.set_enrichment_status!(ids, :radar_status, :unavailable)
else
Radio.set_enrichment_status!(ids, :radar_status, :queued)
end
end
defp mark_radar_status!(_contact, ids, [_ | _]), do: Radio.set_enrichment_status!(ids, :radar_status, :queued)
# Pre-2014 contacts have no HRRR data — NarrClient.in_coverage?/1 is the
# authoritative "HRRR can never serve this" check. Mark :unavailable so
# BackfillEnqueueWorker's :narr filter picks them up and the reconcile
@ -158,6 +177,17 @@ defmodule Microwaveprop.Workers.ContactWeatherEnqueueWorker do
|> Enum.uniq_by(fn changeset -> changeset.changes.args end)
end
def build_radar_jobs(contacts) do
contacts
|> Enum.filter(fn contact ->
contact.pos1 && contact.pos2 && not NarrClient.in_coverage?(contact.qso_timestamp)
end)
|> Enum.map(fn contact ->
CommonVolumeRadarWorker.new(%{"contact_id" => contact.id})
end)
|> Enum.uniq_by(fn changeset -> changeset.changes.args end)
end
def build_weather_jobs(contacts) do
contacts
|> Enum.flat_map(&jobs_for_contact/1)

View file

@ -49,6 +49,8 @@ defmodule MicrowavepropWeb.ContactLive.Show do
narr: nil,
narr_path: [],
iemre: nil,
radar: nil,
mechanism: nil,
terrain: nil,
elevation_profile: nil,
propagation_analysis: nil,
@ -95,6 +97,44 @@ defmodule MicrowavepropWeb.ContactLive.Show do
|> start_async(:native_profile, fn -> load_native_profile(contact) end)
|> start_async(:terrain, fn -> Terrain.get_terrain_profile(contact.id) end)
|> start_async(:iemre, fn -> load_iemre(contact) end)
|> start_async(:radar, fn -> load_radar(contact) end)
end
defp load_radar(contact) do
import Ecto.Query
alias Microwaveprop.Propagation.RainScatterClassifier
alias Microwaveprop.Radio.ContactCommonVolumeRadar
alias Microwaveprop.Repo
radar = Repo.get_by(ContactCommonVolumeRadar, contact_id: contact.id)
hrrr_path = Weather.hrrr_profiles_for_path(contact)
duct_either? = Enum.any?(hrrr_path, &(&1 && &1.ducting_detected))
band_mhz = Decimal.to_integer(Decimal.round(contact.band, 0))
distance_km =
case contact.distance_km do
nil -> 0.0
d -> Decimal.to_float(d)
end
mechanism =
RainScatterClassifier.classify(%{
band_mhz: band_mhz,
distance_km: distance_km,
radar:
if(radar,
do: %{
max_dbz: radar.max_dbz,
heavy_rain_pixel_count: radar.heavy_rain_pixel_count,
coverage_pct: radar.coverage_pct
}
),
duct_either_endpoint: duct_either?
})
%{row: radar, mechanism: mechanism}
end
# The native HRRR profile gives the skew-T a full-atmosphere trace (50
@ -356,6 +396,10 @@ defmodule MicrowavepropWeb.ContactLive.Show do
{:noreply, assign(socket, :iemre, iemre)}
end
def handle_async(:radar, {:ok, %{row: row, mechanism: mechanism}}, socket) do
{:noreply, assign(socket, radar: row, mechanism: mechanism)}
end
def handle_async(:native_profile, {:ok, native_profile}, socket) do
{:noreply, assign(socket, :native_profile, native_profile)}
end
@ -967,6 +1011,33 @@ defmodule MicrowavepropWeb.ContactLive.Show do
<div class="divider" />
<h2 class="text-base font-semibold mb-2">Propagation Mechanism</h2>
<div class="mb-6 border border-base-300 rounded-lg p-4 text-sm">
<div class="flex items-center gap-3 flex-wrap">
<span class={["badge badge-sm", mechanism_badge_class(@mechanism)]}>
{mechanism_label(@mechanism)}
</span>
<span class="opacity-70">{mechanism_explainer(@mechanism)}</span>
</div>
<%= if @radar do %>
<div class="mt-3 flex flex-wrap gap-4 text-xs opacity-80">
<span>
Common volume: {format_number(@radar.common_volume_km2, 0)} km²
</span>
<%= if @radar.max_dbz do %>
<span>Max reflectivity: {format_number(@radar.max_dbz, 1)} dBZ</span>
<% end %>
<span>Heavy-rain pixels: {@radar.heavy_rain_pixel_count}</span>
<%= if @radar.coverage_pct do %>
<span>Radar coverage: {format_number(@radar.coverage_pct, 0)}%</span>
<% end %>
<span class="opacity-60">
Frame: {Calendar.strftime(@radar.observed_at, "%Y-%m-%d %H:%M UTC")}
</span>
</div>
<% end %>
</div>
<h2 class="text-base font-semibold mb-2" id="terrain-heading">Terrain Profile</h2>
<%= if @terrain do %>
<div class="mb-6 border border-base-300 rounded-lg" id="terrain-profile">
@ -1450,6 +1521,45 @@ defmodule MicrowavepropWeb.ContactLive.Show do
defp format_number(n) when is_float(n), do: :erlang.float_to_binary(n, decimals: 1)
defp format_number(n), do: to_string(n)
defp format_number(nil, _), do: ""
defp format_number(n, decimals) when is_float(n) and is_integer(decimals) do
:erlang.float_to_binary(n, decimals: decimals)
end
defp format_number(n, _), do: to_string(n)
defp mechanism_label(:likely_rainscatter), do: "Likely Rain Scatter"
defp mechanism_label(:rainscatter_possible), do: "Rain Scatter Possible"
defp mechanism_label(:tropo_duct), do: "Tropospheric Duct"
defp mechanism_label(:troposcatter), do: "Troposcatter"
defp mechanism_label(:unknown), do: "Unclassified"
defp mechanism_label(_), do: "Analyzing…"
defp mechanism_badge_class(:likely_rainscatter), do: "badge-info"
defp mechanism_badge_class(:rainscatter_possible), do: "badge-info badge-outline"
defp mechanism_badge_class(:tropo_duct), do: "badge-warning"
defp mechanism_badge_class(:troposcatter), do: "badge-ghost"
defp mechanism_badge_class(:unknown), do: "badge-ghost"
defp mechanism_badge_class(_), do: "badge-ghost"
defp mechanism_explainer(:likely_rainscatter),
do:
"Heavy precipitation detected inside the common volume between the two stations — the most likely carrier at this frequency and distance."
defp mechanism_explainer(:rainscatter_possible),
do: "Light-to-moderate precipitation in the common volume, no ducting signature — scatter is plausible but marginal."
defp mechanism_explainer(:tropo_duct),
do: "HRRR profile shows a trapping layer (refractivity gradient below 157 N/km) at one of the endpoints."
defp mechanism_explainer(:troposcatter),
do: "No duct, no meaningful rain inside the common volume — default tropospheric forward scatter."
defp mechanism_explainer(:unknown), do: "Not enough radar/profile coverage to distinguish mechanisms."
defp mechanism_explainer(_), do: ""
defp wind_speed_from_components(u, v) when is_number(u) and is_number(v) do
mps = :math.sqrt(u * u + v * v)
Float.round(mps * 1.944, 1)

View file

@ -0,0 +1,101 @@
defmodule Mix.Tasks.RadarBackfill do
@shortdoc "Enqueue CommonVolumeRadarWorker for contacts with radar_status :pending"
@moduledoc """
Bulk-enqueue per-contact common-volume radar enrichment for historical
contacts.
Eligible contacts: post-2014-10-02 (outside NARR coverage), both pos1
and pos2 present, `radar_status = :pending`. The worker's own unique
constraint on `contact_id` makes re-runs idempotent.
Options:
--limit N Maximum contacts to enqueue per run (default: all)
--year N Only enqueue contacts from this year (optional filter)
--dry-run Print counts; don't enqueue
Example:
mix radar_backfill --year 2024 --limit 5000
"""
use Mix.Task
import Ecto.Query
alias Microwaveprop.Radio
alias Microwaveprop.Radio.Contact
alias Microwaveprop.Repo
alias Microwaveprop.Workers.CommonVolumeRadarWorker
@narr_cutoff ~U[2014-10-02 00:00:00Z]
@chunk 400
@impl Mix.Task
def run(argv) do
Mix.Task.run("app.start")
{opts, _, _} =
OptionParser.parse(argv,
switches: [limit: :integer, year: :integer, dry_run: :boolean]
)
limit = Keyword.get(opts, :limit)
year = Keyword.get(opts, :year)
dry_run? = Keyword.get(opts, :dry_run, false)
total = count_eligible(year)
Mix.shell().info("Eligible contacts: #{total}#{if year, do: " (year #{year})", else: ""}")
if dry_run? or total == 0 do
:ok
else
ids = fetch_ids(year, limit)
Mix.shell().info("Enqueueing #{length(ids)} CommonVolumeRadarWorker jobs...")
ids
|> Enum.chunk_every(@chunk)
|> Enum.with_index(1)
|> Enum.each(fn {batch, i} ->
jobs = Enum.map(batch, &CommonVolumeRadarWorker.new(%{"contact_id" => &1}))
Oban.insert_all(jobs)
Radio.set_enrichment_status!(batch, :radar_status, :queued)
Mix.shell().info(" batch #{i}: #{length(batch)} enqueued")
end)
Mix.shell().info("Done.")
end
end
defp count_eligible(year) do
Contact
|> eligible_query(year)
|> Repo.aggregate(:count, :id)
end
defp fetch_ids(year, limit) do
query =
Contact
|> eligible_query(year)
|> order_by(asc: :qso_timestamp)
|> select([c], c.id)
query = if limit, do: limit(query, ^limit), else: query
Repo.all(query)
end
defp eligible_query(query, year) do
query
|> where([c], not is_nil(c.pos1) and not is_nil(c.pos2))
|> where([c], c.qso_timestamp >= ^@narr_cutoff)
|> where([c], c.radar_status == :pending)
|> maybe_filter_year(year)
end
defp maybe_filter_year(query, nil), do: query
defp maybe_filter_year(query, year) do
start_dt = DateTime.new!(Date.new!(year, 1, 1), ~T[00:00:00], "Etc/UTC")
end_dt = DateTime.new!(Date.new!(year, 12, 31), ~T[23:59:59], "Etc/UTC")
where(query, [c], c.qso_timestamp >= ^start_dt and c.qso_timestamp <= ^end_dt)
end
end

View file

@ -0,0 +1,36 @@
defmodule Microwaveprop.Repo.Migrations.CreateContactCommonVolumeRadar do
use Ecto.Migration
def change do
create table(:contact_common_volume_radar, primary_key: false) do
add :id, :binary_id, primary_key: true
add :contact_id, references(:contacts, type: :binary_id, on_delete: :delete_all),
null: false
add :observed_at, :utc_datetime, null: false
add :common_volume_km2, :float
add :pixel_count, :integer, null: false, default: 0
add :rain_pixel_count, :integer, null: false, default: 0
add :heavy_rain_pixel_count, :integer, null: false, default: 0
add :core_pixel_count, :integer, null: false, default: 0
add :max_dbz, :float
add :mean_dbz, :float
add :coverage_pct, :float
timestamps(type: :utc_datetime)
end
create unique_index(:contact_common_volume_radar, [:contact_id])
create index(:contact_common_volume_radar, [:observed_at])
alter table(:contacts) do
add :radar_status, :string, null: false, default: "pending"
end
create index(:contacts, [:radar_status],
where: "radar_status::text <> 'complete'::text",
name: :contacts_radar_status_pending_index
)
end
end

View file

@ -0,0 +1,102 @@
defmodule Microwaveprop.Propagation.CommonVolumeTest do
use ExUnit.Case, async: true
alias Microwaveprop.Propagation.CommonVolume
describe "in_common_volume?/4" do
test "point at midpoint of two nearby stations is in the common volume" do
pos1 = {32.9, -97.0}
pos2 = {33.1, -96.8}
midpoint = {33.0, -96.9}
assert CommonVolume.in_common_volume?(pos1, pos2, midpoint, 400.0)
end
test "point far from both stations is not in the common volume" do
pos1 = {32.9, -97.0}
pos2 = {33.1, -96.8}
far = {50.0, -50.0}
refute CommonVolume.in_common_volume?(pos1, pos2, far, 400.0)
end
test "point inside circle-1 but outside circle-2 is not in the common volume" do
pos1 = {33.0, -97.0}
# ~900 km east so radius=400 circles don't overlap there
pos2 = {33.0, -87.0}
near_pos1 = {33.0, -96.5}
assert CommonVolume.in_common_volume?(pos1, pos1, near_pos1, 400.0)
refute CommonVolume.in_common_volume?(pos1, pos2, near_pos1, 400.0)
end
test "endpoints themselves are in the common volume when stations are < 2R apart" do
pos1 = {32.9, -97.0}
pos2 = {33.1, -96.8}
assert CommonVolume.in_common_volume?(pos1, pos2, pos1, 400.0)
assert CommonVolume.in_common_volume?(pos1, pos2, pos2, 400.0)
end
test "nothing is in the common volume when stations are > 2R apart" do
pos1 = {33.0, -97.0}
# ~2200 km east
pos2 = {33.0, -73.0}
refute CommonVolume.in_common_volume?(pos1, pos2, pos1, 400.0)
end
end
describe "bounding_box/3" do
test "returns min/max lat/lon covering the entire common volume" do
pos1 = {32.9, -97.0}
pos2 = {33.1, -96.8}
bbox = CommonVolume.bounding_box(pos1, pos2, 400.0)
# Sanity: midpoint is inside bbox
assert bbox.min_lat <= 33.0 and bbox.max_lat >= 33.0
assert bbox.min_lon <= -96.9 and bbox.max_lon >= -96.9
# Bbox shouldn't extend past either individual circle bbox
r_deg_lat = 400.0 / 111.32
assert bbox.min_lat >= (pos1 |> elem(0) |> min(elem(pos2, 0))) - r_deg_lat - 0.01
assert bbox.max_lat <= (pos1 |> elem(0) |> max(elem(pos2, 0))) + r_deg_lat + 0.01
end
test "returns an empty bbox when the circles don't overlap" do
pos1 = {33.0, -97.0}
pos2 = {33.0, -73.0}
bbox = CommonVolume.bounding_box(pos1, pos2, 400.0)
assert bbox == :empty
end
end
describe "area_km2/3" do
test "two coincident circles -> area of full circle" do
pos = {33.0, -97.0}
r = 400.0
assert_in_delta CommonVolume.area_km2(pos, pos, r), :math.pi() * r * r, 1.0
end
test "circles just touching -> small area compared to full circle" do
# ~800 km apart, radius 400 -> barely overlapping
pos1 = {33.0, -97.0}
pos2 = {33.0, -88.48}
area = CommonVolume.area_km2(pos1, pos2, 400.0)
full = :math.pi() * 400.0 * 400.0
assert area >= 0.0
assert area / full < 0.01
end
test "non-overlapping circles -> zero area" do
pos1 = {33.0, -97.0}
pos2 = {33.0, -73.0}
assert CommonVolume.area_km2(pos1, pos2, 400.0) == 0.0
end
test "partially overlapping circles -> positive area less than a full circle" do
pos1 = {33.0, -97.0}
# ~400 km apart -> roughly 1/2 circle overlap
pos2 = {33.0, -92.7}
area = CommonVolume.area_km2(pos1, pos2, 400.0)
full = :math.pi() * 400.0 * 400.0
assert area > 0.0
assert area < full
end
end
end

View file

@ -0,0 +1,62 @@
defmodule Microwaveprop.Propagation.RainScatterClassifierTest do
use ExUnit.Case, async: true
alias Microwaveprop.Propagation.RainScatterClassifier
defp inputs(overrides \\ %{}) do
base = %{
band_mhz: 10_000,
distance_km: 250.0,
radar: %{max_dbz: 42.0, heavy_rain_pixel_count: 12, coverage_pct: 80.0},
duct_either_endpoint: false
}
Map.merge(base, overrides)
end
describe "classify/1" do
test "10 GHz, short path, heavy rain in CV, no duct -> :likely_rainscatter" do
assert RainScatterClassifier.classify(inputs()) == :likely_rainscatter
end
test "duct at either endpoint overrides -> :tropo_duct even with rain" do
assert RainScatterClassifier.classify(inputs(%{duct_either_endpoint: true})) == :tropo_duct
end
test "no rain, no duct, 10 GHz -> :troposcatter" do
inp = inputs(%{radar: %{max_dbz: 5.0, heavy_rain_pixel_count: 0, coverage_pct: 90.0}})
assert RainScatterClassifier.classify(inp) == :troposcatter
end
test "path too long (> 800 km) -> :troposcatter regardless of rain" do
assert RainScatterClassifier.classify(inputs(%{distance_km: 1000.0})) == :troposcatter
end
test "band above rainscatter window (24 GHz) -> not rainscatter" do
result = RainScatterClassifier.classify(inputs(%{band_mhz: 24_000}))
assert result in [:tropo_duct, :troposcatter, :unknown]
refute result == :likely_rainscatter
end
test "band below 5 GHz -> not rainscatter (scattering efficiency drops fast)" do
result = RainScatterClassifier.classify(inputs(%{band_mhz: 1_296}))
refute result == :likely_rainscatter
end
test "no radar coverage at all -> :unknown (can't distinguish)" do
inp = inputs(%{radar: nil, duct_either_endpoint: false})
assert RainScatterClassifier.classify(inp) == :unknown
end
test "radar coverage is too spotty -> :unknown" do
inp = inputs(%{radar: %{max_dbz: 42.0, heavy_rain_pixel_count: 1, coverage_pct: 5.0}})
assert RainScatterClassifier.classify(inp) == :unknown
end
test "light rain (25-35 dBZ) counts as rainscatter-possible, not probable" do
inp = inputs(%{radar: %{max_dbz: 28.0, heavy_rain_pixel_count: 0, coverage_pct: 80.0}})
# 28 dBZ -> rainscatter_possible; absence of heavy cells is the discriminator
assert RainScatterClassifier.classify(inp) == :rainscatter_possible
end
end
end

View file

@ -0,0 +1,116 @@
defmodule Microwaveprop.Workers.CommonVolumeRadarWorkerTest do
use Microwaveprop.DataCase, async: false
use Oban.Testing, repo: Microwaveprop.Repo
alias Microwaveprop.Radio.Contact
alias Microwaveprop.Radio.ContactCommonVolumeRadar
alias Microwaveprop.Repo
alias Microwaveprop.Workers.CommonVolumeRadarWorker
setup do
Req.Test.set_req_test_from_context(%{async: false})
:ok
end
defp insert_contact(attrs \\ %{}) do
default = %{
station1: "W5XD",
station2: "K5TR",
qso_timestamp: ~U[2024-09-15 18:30:00Z],
mode: "CW",
band: Decimal.new("10000"),
grid1: "EM12",
grid2: "EM00",
pos1: %{"lat" => 32.9, "lon" => -97.0},
pos2: %{"lat" => 33.1, "lon" => -96.8},
distance_km: Decimal.new("30")
}
{:ok, contact} =
%Contact{}
|> Contact.changeset(Map.merge(default, attrs))
|> Repo.insert()
contact
end
describe "perform/1" do
test "marks radar_status :unavailable and does not insert a row when n0q frame is 404" do
contact = insert_contact()
Req.Test.stub(Microwaveprop.Weather.NexradClient, fn conn ->
Plug.Conn.send_resp(conn, 404, "not found")
end)
assert :ok =
perform_job(CommonVolumeRadarWorker, %{"contact_id" => contact.id})
refute Repo.get_by(ContactCommonVolumeRadar, contact_id: contact.id)
assert %Contact{radar_status: :unavailable} = Repo.get!(Contact, contact.id)
end
test "skips contacts without both positions and marks unavailable" do
contact = insert_contact(%{pos1: nil})
assert :ok =
perform_job(CommonVolumeRadarWorker, %{"contact_id" => contact.id})
refute Repo.get_by(ContactCommonVolumeRadar, contact_id: contact.id)
assert %Contact{radar_status: :unavailable} = Repo.get!(Contact, contact.id)
end
test "skips contacts beyond 2R so we don't waste a fetch" do
# 2R = 800 km. Put the stations 2000 km apart.
contact =
insert_contact(%{
pos1: %{"lat" => 33.0, "lon" => -97.0},
pos2: %{"lat" => 33.0, "lon" => -75.0}
})
assert :ok =
perform_job(CommonVolumeRadarWorker, %{"contact_id" => contact.id})
refute Repo.get_by(ContactCommonVolumeRadar, contact_id: contact.id)
assert %Contact{radar_status: :unavailable} = Repo.get!(Contact, contact.id)
end
end
describe "aggregate_stats/4 (pure)" do
test "returns zeroed stats when no pixels fall inside the common volume" do
# A 1-pixel-wide synthetic buffer far from the stations — shouldn't match.
# Simplest path: check stations 2000 km apart -> :empty bbox -> zero stats.
pos1 = {33.0, -97.0}
pos2 = {33.0, -75.0}
stats = CommonVolumeRadarWorker.aggregate_stats(<<0, 0, 0>>, 3, pos1, pos2, step: 1)
assert stats.pixel_count == 0
assert stats.rain_pixel_count == 0
assert stats.heavy_rain_pixel_count == 0
assert stats.max_dbz == nil
end
test "counts rain, heavy rain, and core pixels from a synthetic CV patch" do
# Drop two stations close together near the top-left of the image so
# a hand-built tiny buffer falls inside the common volume. We reuse
# NexradClient's lat/lon↔pixel mapping: (50 N, -126 W) is pixel (0, 0).
pos1 = {49.99, -125.99}
pos2 = {49.98, -125.98}
# 5x5 buffer, pixel values chosen so pixel_to_dbz/1 gives a known
# mix of 0 dBZ (no echo), ~30 dBZ (rain), ~40 dBZ (heavy), and
# ~50 dBZ (core). Pixel value 1 -> -30 dBZ, value 255 -> +95 dBZ.
# We just need a couple of non-zero pixels in the right spots.
# Value 123 ≈ 30 dBZ, 143 ≈ 40 dBZ, 163 ≈ 50 dBZ.
pixels =
<<0, 123, 0, 0, 0, 0, 143, 163, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0>>
stats = CommonVolumeRadarWorker.aggregate_stats(pixels, 5, pos1, pos2, step: 1)
assert stats.pixel_count >= 3
assert stats.rain_pixel_count >= 3
assert stats.heavy_rain_pixel_count >= 2
assert stats.core_pixel_count >= 1
assert stats.max_dbz >= 49.0
end
end
end

View file

@ -41,6 +41,7 @@ defmodule MicrowavepropWeb.ConnCase do
setup tags do
DataCase.setup_sandbox(tags)
DataCase.reset_score_files()
DataCase.stub_nexrad_default()
GridCache.clear()
ScoreCache.clear()
{:ok, conn: Phoenix.ConnTest.build_conn()}

View file

@ -32,9 +32,22 @@ defmodule Microwaveprop.DataCase do
setup tags do
Microwaveprop.DataCase.setup_sandbox(tags)
Microwaveprop.DataCase.reset_score_files()
Microwaveprop.DataCase.stub_nexrad_default()
:ok
end
@doc """
Installs a default 404 stub for the IEM NEXRAD client so tests that
trigger `CommonVolumeRadarWorker` via `enqueue_for_contact/1` don't
crash for lack of a Req.Test plug. Individual tests can override with
their own `Req.Test.stub/2`.
"""
def stub_nexrad_default do
Req.Test.stub(Microwaveprop.Weather.NexradClient, fn conn ->
Plug.Conn.send_resp(conn, 404, "not found")
end)
end
@doc """
Sets up the sandbox based on the test tags.
"""