prop/lib/microwaveprop/weather/hrrr_client.ex
Graham McIntire a8bd29f78c
fix(logging): log silently dropped parse failures in grib extractor/wgrib2
- extractor: log warning when point falls outside Lambert grid bounds
- wgrib2: log warning when inventory lines are unparseable
- hrrr_client: optimize wgrib2_match_pattern with uniq_by + map_join
- gefs_client: same optimization
2026-05-07 12:15:26 -05:00

680 lines
22 KiB
Elixir
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

defmodule Microwaveprop.Weather.HrrrClient do
@moduledoc false
alias Microwaveprop.Propagation.Grid
alias Microwaveprop.Weather.Grib2.Extractor
alias Microwaveprop.Weather.Grib2.Wgrib2
require Logger
@hrrr_base_default "https://noaa-hrrr-bdp-pds.s3.amazonaws.com"
defp hrrr_base, do: Application.get_env(:microwaveprop, :hrrr_base_url, @hrrr_base_default)
# Per-point profile fetches (contact detail + skew-T): full troposphere +
# lower stratosphere. 25 levels from 1000→100 mb — fine-grained near the
# surface for duct detection, loosening above 700 mb for the skew-T plot.
@profile_pressure_levels [
1000,
975,
950,
925,
900,
875,
850,
825,
800,
775,
750,
725,
700,
650,
600,
550,
500,
450,
400,
350,
300,
250,
200,
150,
100
]
# Grid-wide fetches (PropagationGridWorker hot path): narrow to the
# near-surface band the scorer actually reads. `SoundingParams.derive`
# pulls `min_refractivity_gradient` from the lowest ~1 km, and native
# hybrid-sigma data (when present) takes priority for that metric
# anyway. Anything above 700 mb is dead weight that doubled the grid
# worker's peak memory when 92k points × levels × vars all decode
# through wgrib2 — enough to OOM the pod at 4 Gi. Keep the hot path
# on this slimmer list.
@grid_pressure_levels [
1000,
975,
950,
925,
900,
875,
850,
825,
800,
775,
750,
725,
700
]
@surface_messages [
%{var: "TMP", level: "2 m above ground"},
%{var: "DPT", level: "2 m above ground"},
%{var: "PRES", level: "surface"},
%{var: "HPBL", level: "surface"},
%{var: "PWAT", level: "entire atmosphere (considered as a single layer)"},
%{var: "UGRD", level: "10 m above ground"},
%{var: "VGRD", level: "10 m above ground"},
%{var: "TCDC", level: "entire atmosphere"},
%{var: "APCP", level: "surface"}
]
# --- Public API ---
@spec surface_messages() :: [%{var: String.t(), level: String.t()}]
def surface_messages, do: @surface_messages
@spec fetch_grid([{float(), float()}], DateTime.t(), keyword()) ::
{:ok, %{{float(), float()} => map()}} | {:error, term()}
def fetch_grid(points, run_time, opts \\ []) do
Microwaveprop.Instrument.span(
[:hrrr, :fetch_grid],
%{point_count: length(points), forecast_hour: Keyword.get(opts, :forecast_hour, 0)},
fn -> do_fetch_grid(points, run_time, opts) end
)
end
defp do_fetch_grid(points, run_time, opts) do
hour_dt = nearest_hrrr_hour(run_time)
date = DateTime.to_date(hour_dt)
hour = hour_dt.hour
forecast_hour = Keyword.get(opts, :forecast_hour, 0)
# Surface and pressure products live in separate GRIB files on
# NOAA's side (wrfsfcf vs wrfprsf), so the whole fetch → wgrib2
# extract pipeline is independent. Running them in parallel
# roughly halves per-forecast-hour wall time (both phases are
# ~30s each; sequential was ~60s, concurrent is max of the two).
sfc_task = Task.async(fn -> fetch_product_grid(date, hour, :surface, points, forecast_hour) end)
prs_task = Task.async(fn -> maybe_fetch_pressure_grid(date, hour, points, opts, forecast_hour) end)
# Match the 20-min per-forecast-hour ceiling in PropagationGridWorker
# so a stuck wgrib2 doesn't pin the chain past its timeout. Surface
# alone deciding the outcome — pressure is optional.
case Task.await(sfc_task, 20 * 60 * 1000) do
{:ok, sfc_grid} ->
prs_grid =
case Task.await(prs_task, 20 * 60 * 1000) do
{:ok, grid} ->
grid
{:error, reason} ->
Logger.warning("HRRR pressure fetch failed (continuing with surface only): #{inspect(reason)}")
%{}
end
# Merge and build profiles in one pass to avoid a full intermediate map copy.
# Use Map.merge to include any pressure-only points, then map directly to profiles.
profiles =
sfc_grid
|> Map.merge(prs_grid, fn _point, sfc_data, prs_data -> Map.merge(sfc_data, prs_data) end)
|> Map.new(fn {point, data} ->
profile = build_profile(data)
{point, Map.put(profile, :run_time, hour_dt)}
end)
{:ok, profiles}
{:error, _reason} = error ->
# Drain the pressure task so it doesn't leak past this scope.
_ = Task.await(prs_task, 20 * 60 * 1000)
error
end
end
@doc """
Per-point HRRR profile fetch returning the full 25 pressure levels
from 1000 → 100 mb. Treats `run_time` as the cycle hour; pass
`forecast_hour: N` to fetch the F+N forecast file from that cycle
(so `valid_time = run_time + N hours`). Default `forecast_hour: 0`
matches the analysis-only behaviour of the original 3-arg form.
"""
@spec fetch_profile(float(), float(), DateTime.t(), keyword()) :: {:ok, map()} | {:error, term()}
def fetch_profile(lat, lon, run_time, opts \\ []) do
forecast_hour = Keyword.get(opts, :forecast_hour, 0)
Microwaveprop.Instrument.span(
[:hrrr, :fetch_profile],
%{lat: lat, lon: lon, forecast_hour: forecast_hour},
fn -> do_fetch_profile(lat, lon, run_time, forecast_hour) end
)
end
defp do_fetch_profile(lat, lon, run_time, forecast_hour) do
hour_dt = nearest_hrrr_hour(run_time)
date = DateTime.to_date(hour_dt)
hour = hour_dt.hour
with {:ok, sfc_data} <- fetch_product(date, hour, :surface, lat, lon, forecast_hour) do
# Pressure is optional — old HRRR data may have corrupt messages
prs_data =
case fetch_product(date, hour, :pressure, lat, lon, forecast_hour) do
{:ok, data} ->
data
{:error, reason} ->
Logger.warning(
"HRRR pressure fetch failed for #{lat},#{lon} (continuing with surface only): #{inspect(reason)}"
)
%{}
end
merged = Map.merge(sfc_data, prs_data)
result = build_profile(merged)
{:ok, Map.merge(result, %{run_time: hour_dt, forecast_hour: forecast_hour})}
end
end
@doc """
Returns true if the HRRR cycle for `run_time` is fully published.
Probes the `wrfprsf18.grib2.idx` file (last forecast hour of the slower
pressure-level product — last to land for a given cycle). A 200 means
every earlier forecast hour is also on disk; anything else means the
cycle is still being uploaded or NOAA hasn't started it yet.
Used by `PropagationGridWorker` to prefer a fresh `now-1h` cycle and
fall back to the conservative `now-2h` slot only when the fresher
cycle hasn't fully published.
Override the probe in tests via the
`:microwaveprop, :hrrr_cycle_available_fn` Application env (a 1-arity
function from `DateTime` to `boolean`).
"""
@spec cycle_available?(DateTime.t()) :: boolean()
def cycle_available?(run_time) do
rounded = nearest_hrrr_hour(run_time)
case Application.get_env(:microwaveprop, :hrrr_cycle_available_fn) do
fun when is_function(fun, 1) -> fun.(rounded)
_ -> probe_cycle(rounded)
end
end
defp probe_cycle(%DateTime{} = run_time) do
url =
run_time
|> DateTime.to_date()
|> hrrr_url(run_time.hour, :pressure, 18)
|> Kernel.<>(".idx")
case Req.head(url, req_options()) do
{:ok, %{status: 200}} -> true
_ -> false
end
end
@doc """
Truncates a DateTime down to the nearest published HRRR cycle hour.
Always rounds DOWN — rounding to the nearest hour used to point at
the next cycle once we passed `:30`, which produced 404/403 errors
when that cycle had not yet been published by NOAA. Callers that need
a later valid_time should reach it via `forecast_hour` instead.
"""
@spec nearest_hrrr_hour(DateTime.t()) :: DateTime.t()
def nearest_hrrr_hour(dt) do
total_seconds = dt.minute * 60 + dt.second
DateTime.add(dt, -total_seconds, :second)
end
@spec hrrr_url(Date.t(), non_neg_integer(), :surface | :pressure, non_neg_integer()) :: String.t()
def hrrr_url(date, hour, product, forecast_hour \\ 0) do
date_str = Calendar.strftime(date, "%Y%m%d")
hour_str = hour |> Integer.to_string() |> String.pad_leading(2, "0")
fh_str = forecast_hour |> Integer.to_string() |> String.pad_leading(2, "0")
file =
case product do
:surface -> "hrrr.t#{hour_str}z.wrfsfcf#{fh_str}.grib2"
:pressure -> "hrrr.t#{hour_str}z.wrfprsf#{fh_str}.grib2"
end
"#{hrrr_base()}/hrrr.#{date_str}/conus/#{file}"
end
@spec parse_idx(String.t()) :: [%{msg: integer(), offset: integer(), var: String.t(), level: String.t()}]
def parse_idx(text) do
# Drop lines that don't start with "msg:offset:" — NOAA S3 sometimes
# returns an HTML error page as the idx body during brief outages,
# and the strict String.to_integer/1 path raised ArgumentError and
# killed the whole chain step. Skip unparseable lines and keep
# whatever well-formed rows are present.
text
|> String.split("\n")
|> Enum.reject(&(&1 == ""))
|> Enum.flat_map(&parse_idx_line/1)
end
defp parse_idx_line(line) do
parts = String.split(line, ":", parts: 8)
with [msg_s, off_s | _] <- parts,
{msg, ""} <- Integer.parse(msg_s),
{offset, ""} <- Integer.parse(off_s) do
[%{msg: msg, offset: offset, var: Enum.at(parts, 3), level: Enum.at(parts, 4)}]
else
_ -> []
end
end
@spec byte_ranges_for_messages([map()], [%{var: String.t(), level: String.t()}]) :: [
{non_neg_integer(), non_neg_integer()}
]
def byte_ranges_for_messages(idx_entries, wanted) do
Enum.flat_map(wanted, fn w ->
idx_entries
|> Enum.with_index()
|> Enum.flat_map(fn {entry, idx} ->
byte_range_for_entry(entry, idx, idx_entries, w)
end)
end)
end
defp byte_range_for_entry(entry, idx, idx_entries, wanted) do
if entry.var == wanted.var && entry.level == wanted.level do
end_offset =
case Enum.at(idx_entries, idx + 1) do
nil -> entry.offset + 10_000_000
next -> next.offset - 1
end
[{entry.offset, end_offset}]
else
[]
end
end
@spec merge_grid_data(map(), map()) :: map()
def merge_grid_data(sfc_grid, prs_grid) do
Map.merge(sfc_grid, prs_grid, fn _point, sfc_data, prs_data ->
Map.merge(sfc_data, prs_data)
end)
end
@spec build_profile(map()) :: map()
def build_profile(parsed) do
sfc_temp_k = parsed["TMP:2 m above ground"]
sfc_dpt_k = parsed["DPT:2 m above ground"]
sfc_pres_pa = parsed["PRES:surface"]
# Iterate the full profile list so per-contact fetches that include
# upper-air data still populate the whole skew-T trace. Grid fetches
# only populate the 1000→700 mb subset (see @grid_pressure_levels);
# the other levels just produce empty slots that this loop skips.
profile =
Enum.flat_map(@profile_pressure_levels, fn level ->
level_str = "#{level} mb"
tmp = parsed["TMP:#{level_str}"]
dpt = parsed["DPT:#{level_str}"]
hgt = parsed["HGT:#{level_str}"]
if tmp && dpt && hgt do
[
%{
"pres" => level * 1.0,
"tmpc" => tmp - 273.15,
"dwpc" => dpt - 273.15,
"hght" => hgt
}
]
else
[]
end
end)
# Fall back to lowest pressure level if surface fields are missing (near-coast points)
lowest_level = List.first(profile)
sfc_temp_c =
if sfc_temp_k, do: sfc_temp_k - 273.15, else: lowest_level && lowest_level["tmpc"]
sfc_dpt_c =
if sfc_dpt_k, do: sfc_dpt_k - 273.15, else: lowest_level && lowest_level["dwpc"]
%{
surface_temp_c: sfc_temp_c,
surface_dewpoint_c: sfc_dpt_c,
surface_pressure_mb: if(sfc_pres_pa, do: sfc_pres_pa / 100.0),
hpbl_m: parsed["HPBL:surface"],
pwat_mm: parsed["PWAT:entire atmosphere (considered as a single layer)"],
wind_u: parsed["UGRD:10 m above ground"],
wind_v: parsed["VGRD:10 m above ground"],
cloud_cover_pct: parsed["TCDC:entire atmosphere"],
precip_mm: parsed["APCP:surface"],
profile: profile
}
end
# --- Private ---
defp fetch_product(date, hour, product, lat, lon, forecast_hour) do
url = hrrr_url(date, hour, product, forecast_hour)
idx_url = url <> ".idx"
wanted =
case product do
:surface -> @surface_messages
:pressure -> pressure_messages(:profile)
end
Logger.info("HRRR fetching #{product} idx from #{idx_url}")
with {:ok, idx_text} <- fetch_idx(idx_url),
idx_entries = parse_idx(idx_text),
ranges = byte_ranges_for_messages(idx_entries, wanted),
_ = Logger.info("HRRR downloading #{length(ranges)} GRIB ranges for #{product}"),
{:ok, grib_binary} <- download_grib_ranges(url, ranges) do
Logger.info("HRRR #{product} downloaded #{byte_size(grib_binary)} bytes, extracting")
Extractor.extract_points(grib_binary, lat, lon)
end
end
defp fetch_product_grid(date, hour, product, _points, forecast_hour) do
url = hrrr_url(date, hour, product, forecast_hour)
idx_url = url <> ".idx"
wanted =
case product do
:surface -> @surface_messages
:pressure -> pressure_messages(:grid)
end
Logger.info("HRRR grid fetching #{product} idx from #{idx_url}")
with {:ok, idx_text} <- fetch_idx(idx_url),
idx_entries = parse_idx(idx_text),
ranges = byte_ranges_for_messages(idx_entries, wanted),
_ = Logger.info("HRRR grid downloading #{length(ranges)} GRIB ranges for #{product}"),
{:ok, grib_binary} <- download_grib_ranges(url, ranges) do
Logger.info("HRRR grid #{product} downloaded #{byte_size(grib_binary)} bytes, extracting")
extract_grid(grib_binary, wanted, product)
end
end
defp extract_grid(grib_binary, wanted, product) do
if Wgrib2.available?() do
match_pattern = wgrib2_match_pattern(wanted)
grid_spec = Grid.wgrib2_grid_spec()
Logger.info("HRRR grid using wgrib2 for #{product} extraction")
Wgrib2.extract_grid(grib_binary, match_pattern, grid_spec)
else
Logger.info("HRRR grid using Elixir decoder for #{product} extraction (slow)")
points = Grid.conus_points()
Extractor.extract_grid(grib_binary, points)
end
end
defp wgrib2_match_pattern(wanted) do
vars = wanted |> Enum.uniq_by(& &1.var) |> Enum.map_join("|", & &1.var)
":(#{vars}):"
end
defp maybe_fetch_pressure_grid(date, hour, points, opts, forecast_hour) do
if Keyword.get(opts, :pressure, true) do
fetch_product_grid(date, hour, :pressure, points, forecast_hour)
else
{:ok, %{}}
end
end
@doc """
Pressure-level GRIB message descriptors for a fetch variant.
* `:grid` — narrow list of 13 levels (1000→700 mb) used by the
propagation grid hot path. Covers the ducting band without
ballooning the per-forecast-hour memory footprint.
* `:profile` — full 25 levels (1000→100 mb) used for per-contact
skew-T fetches.
"""
@spec pressure_messages(:grid | :profile) :: [%{var: String.t(), level: String.t()}]
def pressure_messages(:grid) do
build_pressure_messages(@grid_pressure_levels)
end
def pressure_messages(:profile) do
build_pressure_messages(@profile_pressure_levels)
end
defp build_pressure_messages(levels) do
for level <- levels, var <- ["TMP", "DPT", "HGT"] do
%{var: var, level: "#{level} mb"}
end
end
# A published HRRR idx file is immutable for the lifetime of its
# model run (60 min), but the hourly propagation chain fetches the
# same `.idx` URL dozens of times across forecast-hour steps and
# products. Memoising in `Microwaveprop.Cache` removes 100+ seconds
# per chain from the `hrrr_fetch_idx` span.
@idx_cache_ttl_ms to_timeout(hour: 1)
defp fetch_idx(url) do
Microwaveprop.Cache.fetch_or_store({:hrrr_idx, url}, @idx_cache_ttl_ms, fn ->
Microwaveprop.Instrument.span([:hrrr, :fetch_idx], %{url: url}, fn -> do_fetch_idx(url) end)
end)
end
defp do_fetch_idx(url) do
case Req.get(url, req_options()) do
{:ok, %{status: 200, body: body}} -> {:ok, body}
{:ok, %{status: status}} -> {:error, "HRRR idx HTTP #{status}"}
{:error, reason} -> {:error, reason}
end
end
@doc """
Downloads the given byte ranges from a GRIB2 URL and concatenates
them in offset order. Merges adjacent/overlapping ranges and runs
requests in parallel. Backed by a disk cache keyed on
`(url, ranges)` in dev.
Public so the native-level client can reuse it rather than
duplicating the HTTP/parallelism/cache logic.
"""
@spec download_grib_ranges(String.t(), [{non_neg_integer(), non_neg_integer()}]) ::
{:ok, binary()} | {:error, term()}
def download_grib_ranges(_url, []), do: {:ok, <<>>}
def download_grib_ranges(url, ranges) do
Microwaveprop.Instrument.span(
[:hrrr, :download_grib_ranges],
%{range_count: length(ranges)},
fn -> do_download_grib_ranges(url, ranges) end
)
end
defp do_download_grib_ranges(url, ranges) do
cache_key = url_to_cache_key(url, ranges)
case read_cache(cache_key) do
{:ok, binary} ->
Logger.info("HRRR cache hit: #{cache_key}")
{:ok, binary}
:miss ->
download_and_cache_grib_ranges(url, ranges, cache_key)
end
end
defp download_and_cache_grib_ranges(url, ranges, cache_key) do
case download_grib_ranges_remote(url, ranges) do
{:ok, binary} = ok ->
write_cache(cache_key, binary)
ok
error ->
error
end
end
@doc """
Like `download_grib_ranges/2`, but writes the result to `dest_path`
instead of returning a binary. Each chunk is written to disk as it
arrives, so peak memory stays low even for ~530 MB native-level files.
Returns `:ok` or `{:error, reason}`.
"""
@spec download_grib_ranges_to_file(String.t(), [{non_neg_integer(), non_neg_integer()}], Path.t()) ::
:ok | {:error, term()}
def download_grib_ranges_to_file(_url, [], _dest_path), do: :ok
def download_grib_ranges_to_file(url, ranges, dest_path) do
case fetch_grib_ranges(url, ranges) do
{:error, _} = error ->
error
{:ok, results} ->
file = File.open!(dest_path, [:write, :binary])
try do
results
|> Enum.sort_by(fn {:ok, offset, _} -> offset end)
|> Enum.each(fn {:ok, _offset, body} -> IO.binwrite(file, body) end)
after
File.close(file)
end
:ok
end
end
defp download_grib_ranges_remote(url, ranges) do
case fetch_grib_ranges(url, ranges) do
{:error, _} = error ->
error
{:ok, results} ->
binary =
results
|> Enum.sort_by(fn {:ok, offset, _} -> offset end)
|> Enum.map(fn {:ok, _, body} -> body end)
|> IO.iodata_to_binary()
{:ok, binary}
end
end
defp fetch_grib_ranges(url, ranges) do
merged = merge_ranges(ranges)
supervisor = {:via, PartitionSupervisor, {Microwaveprop.TaskSupervisor, self()}}
results =
supervisor
|> Task.Supervisor.async_stream_nolink(
merged,
fn {start, stop} ->
case Req.get(url, [{:headers, [{"Range", "bytes=#{start}-#{stop}"}]} | req_options()]) do
{:ok, %{status: 206, body: body}} -> {:ok, start, body}
{:ok, %{status: status}} -> {:error, "HRRR grib HTTP #{status}"}
{:error, reason} -> {:error, reason}
end
end,
max_concurrency: 8,
timeout: 120_000
)
|> Enum.map(fn
{:ok, result} ->
result
{:exit, reason} ->
Logger.error("HRRR grib range download crashed: url=#{url} reason=#{inspect(reason)}")
{:error, reason}
end)
if error = Enum.find(results, &match?({:error, _}, &1)) do
error
else
{:ok, results}
end
end
defp hrrr_cache_dir, do: Application.get_env(:microwaveprop, :hrrr_cache_dir)
defp url_to_cache_key(url, ranges) do
# e.g. hrrr.20240922_conus_hrrr.t18z.wrfprsf00_r42.grib2
uri = URI.parse(url)
path_part = uri.path |> String.trim_leading("/") |> String.replace("/", "_")
range_hash = ranges |> :erlang.phash2() |> Integer.to_string(16)
"#{path_part}_r#{range_hash}.grib2"
end
defp read_cache(key) do
case hrrr_cache_dir() do
nil ->
:miss
dir ->
path = Path.join(dir, key)
case File.read(path) do
{:ok, binary} -> {:ok, binary}
{:error, _} -> :miss
end
end
end
defp write_cache(key, binary) do
case hrrr_cache_dir() do
nil ->
:ok
dir ->
File.mkdir_p!(dir)
path = Path.join(dir, key)
File.write!(path, binary)
Logger.info("HRRR cached: #{path} (#{byte_size(binary)} bytes)")
end
end
defp merge_ranges(ranges) do
ranges
|> Enum.sort_by(&elem(&1, 0))
|> Enum.reduce([], fn {s, e}, acc ->
case acc do
[{ps, pe} | rest] when s <= pe + 1 -> [{ps, max(pe, e)} | rest]
_ -> [{s, e} | acc]
end
end)
|> Enum.reverse()
end
defp req_options do
defaults = [receive_timeout: 120_000, retry: &retry?/2, max_retries: 5, retry_delay: &retry_delay/1]
overrides = Application.get_env(:microwaveprop, :hrrr_req_options, [])
Keyword.merge(defaults, overrides)
end
defp retry?(_request, response) do
case response do
%Req.Response{status: status} when status in [429, 500, 502, 503, 504] -> true
%{__exception__: true} -> true
_ -> false
end
end
defp retry_delay(n) do
base = Integer.pow(2, n) * 1_000
jitter = :rand.uniform(1_000)
base + jitter
end
end