Graham McIntire
0f0e5e8d43
Phase 3 NEXRAD spike: IEM n0q composite available at 5-min cadence back to 2022+. Compression-ratio proxy shows afternoon images have 13-81% more texture than dawn (directionally correct), but the n0q product thresholds out the faint clear-air returns needed for BL stability detection. Parked until MRMS or Level III products can be investigated. See docs/research/nexrad_spike.md. Phase 6: hrrr_climatology table aggregating surface_temp_c by (lat, lon, month, hour) from the 42M+ hrrr_profiles grid-point rows. mix hrrr_climatology builds it via a single SQL GROUP BY + upsert. Backtest.Features.temperature_anomaly computes current_temp minus climatological mean — the meteorologist's "temperature deviation above normal" predictor for summer afternoon enhancement.
283 lines
9.1 KiB
Elixir
283 lines
9.1 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 "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
|