Graham McIntire
a840cb9629
FrontalAnalysis module (Weather.FrontalAnalysis): - detect_fronts/3 computes the Thermal Front Parameter (TFP) from 2D grids of surface temperature and pressure using Nx vectorized ops. TFP = -nabla|nabla(theta)| . nabla(theta)/|nabla(theta)|. Most negative values mark cold fronts. - central_gradient/1 for 2D finite differences with edge handling - nearest_front/3 finds closest front point with distance and bearing - path_front_angle/2 computes angle between a QSO path and the front (0 = parallel = good, 90 = crosses = dead) Backtest feature stubs for distance_to_front and parallel_to_front (return nil until the pipeline caches per-cell frontal features from the hourly HRRR grid run). The FrontalAnalysis module itself is tested and ready for integration. NEXRAD spike docs also included in this commit.
306 lines
10 KiB
Elixir
306 lines
10 KiB
Elixir
defmodule Microwaveprop.Backtest.Features do
|
|
@moduledoc """
|
|
Named feature functions for use with `Microwaveprop.Backtest`.
|
|
|
|
Every feature has the shape `(lat, lon, valid_time) -> float | nil`
|
|
and is named after the physical quantity it represents. These are
|
|
the "known baselines" the plan refers to: wrappers around the
|
|
existing scorer's inputs so we can measure the lift of new features
|
|
against them on an apples-to-apples basis.
|
|
|
|
## Contract
|
|
|
|
- Return a `float` when the underlying data is available.
|
|
- Return `nil` when there's no HRRR profile within the usual match
|
|
window (`Weather.find_nearest_hrrr/3` returning nil).
|
|
- Never raise — bad inputs should produce `nil`, not crashes. The
|
|
backtest harness runs these across tens of thousands of calls and
|
|
a raise on one point kills the whole report.
|
|
"""
|
|
|
|
import Ecto.Query
|
|
|
|
alias Microwaveprop.Repo
|
|
alias Microwaveprop.Weather
|
|
alias Microwaveprop.Weather.HrrrClimatology
|
|
alias Microwaveprop.Weather.HrrrNativeProfile
|
|
|
|
@doc """
|
|
Minimum refractivity gradient from the nearest HRRR profile.
|
|
|
|
This is the scalar the current scorer uses. More negative is better
|
|
(stronger ducting potential). We return the raw gradient; the
|
|
backtest harness handles binning and summarizing.
|
|
"""
|
|
@spec naive_gradient(float, float, DateTime.t()) :: float | nil
|
|
def naive_gradient(lat, lon, valid_time) do
|
|
with %{min_refractivity_gradient: grad} when is_float(grad) <-
|
|
Weather.find_nearest_hrrr(lat, lon, valid_time) do
|
|
grad
|
|
else
|
|
_ -> nil
|
|
end
|
|
end
|
|
|
|
@doc """
|
|
Dewpoint depression (T - Td) at the surface, in °C.
|
|
|
|
Lower depression means higher relative humidity. For 10 GHz the
|
|
existing scorer treats this as beneficial; for 24+ GHz it's harmful.
|
|
"""
|
|
@spec td_depression(float, float, DateTime.t()) :: float | nil
|
|
def td_depression(lat, lon, valid_time) do
|
|
with %{surface_temp_c: t, surface_dewpoint_c: td} when is_float(t) and is_float(td) <-
|
|
Weather.find_nearest_hrrr(lat, lon, valid_time) do
|
|
t - td
|
|
else
|
|
_ -> nil
|
|
end
|
|
end
|
|
|
|
@doc """
|
|
Time-of-day feature: hours since midnight UTC, as a float in [0, 24).
|
|
|
|
A flat-by-time baseline against which diurnal lift is measured. The
|
|
existing scorer collapses this to a band-dependent shape; the
|
|
backtest treats it as raw UTC hour so we can see the shape directly
|
|
in the distribution.
|
|
"""
|
|
@spec time_of_day(float, float, DateTime.t()) :: float
|
|
def time_of_day(_lat, _lon, valid_time) do
|
|
valid_time.hour + valid_time.minute / 60.0
|
|
end
|
|
|
|
@doc """
|
|
Surface pressure in hPa from the nearest HRRR profile.
|
|
|
|
Used as the baseline the plan predicts `ParallelToFront` (Phase 5)
|
|
will replace.
|
|
"""
|
|
@spec pressure(float, float, DateTime.t()) :: float | nil
|
|
def pressure(lat, lon, valid_time) do
|
|
with %{surface_pressure_mb: p} when is_float(p) <-
|
|
Weather.find_nearest_hrrr(lat, lon, valid_time) do
|
|
p
|
|
else
|
|
_ -> nil
|
|
end
|
|
end
|
|
|
|
@doc """
|
|
Surface refractivity from the lowest level of the nearest native
|
|
hybrid-sigma profile.
|
|
|
|
Phase 1 sanity check: this should produce numbers comparable to the
|
|
`naive_gradient` baseline but computed from native-level data. If
|
|
the sign/magnitude look wrong, the native ingestion pipeline has a
|
|
bug.
|
|
|
|
Uses the ITU-R P.453-14 formula: N = 77.6*P/T + 3.73e5*e/T²
|
|
where e (water vapor pressure) is derived from specific humidity.
|
|
"""
|
|
@spec native_surface_refractivity(float, float, DateTime.t()) :: float | nil
|
|
def native_surface_refractivity(lat, lon, valid_time) do
|
|
with %HrrrNativeProfile{} = profile <- find_nearest_native(lat, lon, valid_time),
|
|
t when is_float(t) <- profile.surface_temp_k,
|
|
p when is_float(p) <- profile.surface_pressure_pa,
|
|
q when is_float(q) <- profile.surface_spfh do
|
|
# Water vapor pressure from specific humidity: e = q*P / (0.622 + 0.378*q)
|
|
e = q * p / (0.622 + 0.378 * q)
|
|
# N-units
|
|
77.6 * p / (t * 100) + 3.73e5 * e / (t * t * 100)
|
|
else
|
|
_ -> nil
|
|
end
|
|
end
|
|
|
|
@doc """
|
|
Bulk Richardson number at the inversion top from the nearest native profile.
|
|
|
|
Ri < 0.25 → turbulent (bad for propagation)
|
|
Ri > 1.0 → laminar (good for propagation)
|
|
|
|
Returns nil if no inversion detected or no native profile available.
|
|
"""
|
|
@spec bulk_richardson(float, float, DateTime.t()) :: float | nil
|
|
def bulk_richardson(lat, lon, valid_time) do
|
|
with %HrrrNativeProfile{bulk_richardson: ri} when is_float(ri) <-
|
|
find_nearest_native(lat, lon, valid_time) do
|
|
ri
|
|
else
|
|
_ -> nil
|
|
end
|
|
end
|
|
|
|
@doc """
|
|
θₑ jump (K) across the inversion from the nearest native profile.
|
|
|
|
Larger positive values indicate stronger thermodynamic decoupling.
|
|
"""
|
|
@spec theta_e_jump(float, float, DateTime.t()) :: float | nil
|
|
def theta_e_jump(lat, lon, valid_time) do
|
|
with %HrrrNativeProfile{theta_e_jump_k: jump} when is_float(jump) <-
|
|
find_nearest_native(lat, lon, valid_time) do
|
|
jump
|
|
else
|
|
_ -> nil
|
|
end
|
|
end
|
|
|
|
@doc """
|
|
Wind shear magnitude (m/s) at the inversion top from the nearest native profile.
|
|
"""
|
|
@spec shear_at_top(float, float, DateTime.t()) :: float | nil
|
|
def shear_at_top(lat, lon, valid_time) do
|
|
with %HrrrNativeProfile{shear_at_top_ms: shear} when is_float(shear) <-
|
|
find_nearest_native(lat, lon, valid_time) do
|
|
shear
|
|
else
|
|
_ -> nil
|
|
end
|
|
end
|
|
|
|
@doc """
|
|
Max duct thickness (m) from the nearest native profile.
|
|
|
|
Larger ducts trap lower frequencies and support longer-range propagation.
|
|
Returns nil if no ducts detected or no native profile available.
|
|
"""
|
|
@spec duct_thickness(float, float, DateTime.t()) :: float | nil
|
|
def duct_thickness(lat, lon, valid_time) do
|
|
with %HrrrNativeProfile{ducts: ducts} when is_list(ducts) and ducts != [] <-
|
|
find_nearest_native(lat, lon, valid_time) do
|
|
ducts |> Enum.map(& &1["thickness_m"]) |> Enum.max(fn -> nil end)
|
|
else
|
|
_ -> nil
|
|
end
|
|
end
|
|
|
|
@doc """
|
|
Lowest trapped frequency (GHz) across all ducts in the nearest native profile.
|
|
|
|
A lower value means stronger ducting that can trap more bands.
|
|
"""
|
|
@spec best_duct_freq(float, float, DateTime.t()) :: float | nil
|
|
def best_duct_freq(lat, lon, valid_time) do
|
|
with %HrrrNativeProfile{best_duct_band_ghz: f} when is_float(f) <-
|
|
find_nearest_native(lat, lon, valid_time) do
|
|
f
|
|
else
|
|
_ -> nil
|
|
end
|
|
end
|
|
|
|
@doc """
|
|
Boolean (1.0 or 0.0): is there a duct that can trap the given band?
|
|
|
|
Pass the band frequency in GHz. Returns 1.0 if any duct has
|
|
min_freq_ghz <= band_ghz, 0.0 otherwise, nil if no data.
|
|
"""
|
|
@spec duct_usable_for_band(float, float, DateTime.t(), float) :: float | nil
|
|
def duct_usable_for_band(lat, lon, valid_time, band_ghz) do
|
|
with %HrrrNativeProfile{best_duct_band_ghz: f} when is_float(f) <-
|
|
find_nearest_native(lat, lon, valid_time) do
|
|
if f <= band_ghz, do: 1.0, else: 0.0
|
|
else
|
|
_ -> nil
|
|
end
|
|
end
|
|
|
|
@doc """
|
|
Temperature anomaly (°C): how much hotter or cooler the current
|
|
surface temperature is compared to the climatological mean for
|
|
this grid cell, month, and hour.
|
|
|
|
Positive = hotter than normal. The meteorologist noted that
|
|
10°F+ above-normal summer days produce enhanced ducting even in
|
|
the afternoon.
|
|
"""
|
|
@spec temperature_anomaly(float, float, DateTime.t()) :: float | nil
|
|
def temperature_anomaly(lat, lon, valid_time) do
|
|
with %{surface_temp_c: temp_c} when is_float(temp_c) <-
|
|
Weather.find_nearest_hrrr(lat, lon, valid_time),
|
|
%HrrrClimatology{mean_surface_temp_c: clim_mean} when is_float(clim_mean) <-
|
|
find_climatology(lat, lon, valid_time.month, valid_time.hour) do
|
|
temp_c - clim_mean
|
|
else
|
|
_ -> nil
|
|
end
|
|
end
|
|
|
|
defp find_climatology(lat, lon, month, hour) do
|
|
dlat = 0.07
|
|
dlon = 0.07
|
|
|
|
HrrrClimatology
|
|
|> where(
|
|
[c],
|
|
c.lat >= ^(lat - dlat) and c.lat <= ^(lat + dlat) and
|
|
c.lon >= ^(lon - dlon) and c.lon <= ^(lon + dlon) and
|
|
c.month == ^month and c.hour == ^hour
|
|
)
|
|
|> limit(1)
|
|
|> Repo.one()
|
|
end
|
|
|
|
@doc """
|
|
Distance (km) to the nearest detected front from the nearest HRRR
|
|
grid cell. Requires frontal analysis to have been run for the
|
|
matching HRRR hour — returns nil if not available.
|
|
|
|
NOTE: This feature is a placeholder that returns nil until the
|
|
frontal analysis pipeline (Phase 5.3) caches per-cell frontal
|
|
features in the propagation_scores table or a sidecar. For now
|
|
it documents the intended API; the backtest will only work once
|
|
the pipeline is live.
|
|
"""
|
|
@spec distance_to_front(float, float, DateTime.t()) :: float | nil
|
|
def distance_to_front(_lat, _lon, _valid_time), do: nil
|
|
|
|
@doc """
|
|
cos²(path_front_angle): 1.0 if the path is parallel to the front,
|
|
0.0 if perpendicular. Requires both a path bearing and a front
|
|
bearing, so this is only usable in QSO-level scoring, not grid
|
|
scoring. Placeholder until Phase 5.4.
|
|
"""
|
|
@spec parallel_to_front(float, float, DateTime.t()) :: float | nil
|
|
def parallel_to_front(_lat, _lon, _valid_time), do: nil
|
|
|
|
@doc "Duct usable for 10 GHz."
|
|
def duct_usable_10ghz(lat, lon, vt), do: duct_usable_for_band(lat, lon, vt, 10.0)
|
|
|
|
@doc "Duct usable for 24 GHz."
|
|
def duct_usable_24ghz(lat, lon, vt), do: duct_usable_for_band(lat, lon, vt, 24.0)
|
|
|
|
@doc "Duct usable for 47 GHz."
|
|
def duct_usable_47ghz(lat, lon, vt), do: duct_usable_for_band(lat, lon, vt, 47.0)
|
|
|
|
defp find_nearest_native(lat, lon, valid_time) do
|
|
dlat = 0.07
|
|
dlon = 0.07
|
|
time_start = DateTime.add(valid_time, -3600, :second)
|
|
time_end = DateTime.add(valid_time, 3600, :second)
|
|
|
|
HrrrNativeProfile
|
|
|> where(
|
|
[p],
|
|
p.lat >= ^(lat - dlat) and p.lat <= ^(lat + dlat) and
|
|
p.lon >= ^(lon - dlon) and p.lon <= ^(lon + dlon) and
|
|
p.valid_time >= ^time_start and p.valid_time <= ^time_end
|
|
)
|
|
|> order_by([p],
|
|
asc:
|
|
fragment(
|
|
"ABS(? - ?) + ABS(? - ?) + ABS(EXTRACT(EPOCH FROM ? - ?))",
|
|
p.lat,
|
|
^lat,
|
|
p.lon,
|
|
^lon,
|
|
p.valid_time,
|
|
^valid_time
|
|
)
|
|
)
|
|
|> limit(1)
|
|
|> Repo.one()
|
|
end
|
|
end
|