prop/lib/microwaveprop/weather/era5_client.ex
Graham McIntire fc3eb58910
Refresh propagation algo against expanded prod corpus (2026-04-13)
Direct queries against prod (75M HRRR profiles, 16,864 soundings, 392k
surface obs, new hrrr_native_profiles table) drive these changes:

* Reinstate shallow-BL bonus as an HPBL multiplier on the refractivity
  score (1.10× <200m → 0.78× ≥2000m). The Apr 11 "shallow BL bonus
  removed" conclusion was an artifact of the prior matching strategy;
  with the cleaner contact↔HRRR join (n=680) the binned data goes 230 km
  avg at HPBL <200m vs 100 km at ≥2000m, monotonic across all bins.
* Add 6 missing bands (142, 145, 288, 322, 403, 411 GHz) so contacts in
  those bands stop being silently dropped. Coefficients extrapolate
  ITU-R P.676/P.838 trends from the existing 134/241 GHz entries.
* Bump ERA5 poll timeout 10 min → 1 hour and Era5MonthBatchWorker
  max_attempts 3 → 5 so CDS slowness stops discarding tiles. Also dump
  the full failure body when CDS omits the message field — the prior
  "ERA5 job failed: nil" was masking real error reasons in oban_jobs.
* Add scripts/recalibrate_algo.py so this analysis can be re-run any
  time new data lands without manual SQL. Reads PROP_PROD_DB_URL from
  .envrc, drops a Markdown report into docs/algo-reports/.
* Append a dated Part 2c to algo.md documenting the corpus expansion,
  the honest accounting of historical HRRR coverage (only ~1,020 of 58k
  contacts are precision-matchable), and the empty era5/rtma/climatology
  tables. Update the gaseous-absorption and rain-attenuation tables to
  include the new bands.

Test suite: 1,335 tests, 0 failures.
2026-04-13 10:59:01 -05:00

313 lines
11 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.Era5Client do
@moduledoc """
Client for the Copernicus Climate Data Store (CDS) API.
Fetches ERA5 reanalysis data (pressure-level profiles + single-level surface fields)
for historical contact enrichment where HRRR is unavailable (pre-2014).
The CDS API is asynchronous: submit a request, poll for completion, download result.
Results are GRIB2 format, decoded with the existing GRIB2 infrastructure.
"""
alias Microwaveprop.Weather.Grib2.Extractor
alias Microwaveprop.Weather.SoundingParams
require Logger
@cds_url "https://cds.climate.copernicus.eu/api"
@poll_interval_ms 5_000
# CDS routinely takes 30+ min to return month-tile jobs under load. The prior
# 120-attempt (10 min) ceiling was burning through Era5MonthBatchWorker
# max_attempts before CDS even produced a result. 720 attempts × 5s = 1 hour
# bounds the worker but still gives CDS enough room.
@max_poll_attempts 720
@pressure_levels ~w(1000 975 950 925 900 875 850 825 800 775 750 725 700)
@single_level_vars ~w(
2m_temperature
2m_dewpoint_temperature
surface_pressure
10m_u_component_of_wind
10m_v_component_of_wind
total_column_water_vapour
boundary_layer_height
)
@pressure_level_vars ~w(temperature dewpoint_temperature geopotential)
@doc false
@spec pressure_levels() :: [String.t()]
def pressure_levels, do: @pressure_levels
@doc false
@spec single_level_vars() :: [String.t()]
def single_level_vars, do: @single_level_vars
@doc false
@spec pressure_level_vars() :: [String.t()]
def pressure_level_vars, do: @pressure_level_vars
@doc false
@spec cds_url() :: String.t()
def cds_url, do: @cds_url
@doc false
@spec poll_interval_ms() :: pos_integer()
def poll_interval_ms, do: @poll_interval_ms
@doc false
@spec max_poll_attempts() :: pos_integer()
def max_poll_attempts, do: @max_poll_attempts
@doc """
Fetch ERA5 profile for a single point and time.
Returns {:ok, profile_attrs} or {:error, reason}.
The profile_attrs map has the same shape as HRRR profiles for interoperability.
"""
@spec fetch_profile(float(), float(), DateTime.t()) :: {:ok, map()} | {:error, term()}
def fetch_profile(lat, lon, timestamp) do
# Round to ERA5 grid (0.25°)
rlat = Float.round(lat * 4) / 4
rlon = Float.round(lon * 4) / 4
# ERA5 hourly — round to nearest hour
valid_time = DateTime.truncate(timestamp, :second)
valid_time = %{valid_time | minute: 0, second: 0}
# Build a small bounding box (0.25° around the point)
area = [rlat + 0.25, rlon - 0.25, rlat - 0.25, rlon + 0.25]
with {:ok, single_data} <- fetch_single_levels(valid_time, area),
{:ok, pressure_data} <- fetch_pressure_levels(valid_time, area) do
build_profile(rlat, rlon, valid_time, single_data, pressure_data)
end
end
@doc """
Fetch ERA5 single-level (surface) data for a time and area.
"""
@spec fetch_single_levels(DateTime.t(), [number()]) :: {:ok, binary()} | {:error, term()}
def fetch_single_levels(valid_time, area) do
request = %{
"product_type" => ["reanalysis"],
"variable" => @single_level_vars,
"year" => [Calendar.strftime(valid_time, "%Y")],
"month" => [Calendar.strftime(valid_time, "%m")],
"day" => [Calendar.strftime(valid_time, "%d")],
"time" => [Calendar.strftime(valid_time, "%H:00")],
"area" => area,
"data_format" => "grib"
}
do_submit_and_download("reanalysis-era5-single-levels", request)
end
@doc """
Fetch ERA5 pressure-level data for a time and area.
"""
@spec fetch_pressure_levels(DateTime.t(), [number()]) :: {:ok, binary()} | {:error, term()}
def fetch_pressure_levels(valid_time, area) do
request = %{
"product_type" => ["reanalysis"],
"variable" => @pressure_level_vars,
"pressure_level" => @pressure_levels,
"year" => [Calendar.strftime(valid_time, "%Y")],
"month" => [Calendar.strftime(valid_time, "%m")],
"day" => [Calendar.strftime(valid_time, "%d")],
"time" => [Calendar.strftime(valid_time, "%H:00")],
"area" => area,
"data_format" => "grib"
}
do_submit_and_download("reanalysis-era5-pressure-levels", request)
end
@doc false
@spec submit_and_download(String.t(), map()) :: {:ok, binary()} | {:error, term()}
def submit_and_download(dataset, request) do
do_submit_and_download(dataset, request)
end
defp do_submit_and_download(dataset, request) do
api_key = api_key()
if is_nil(api_key) do
{:error, "ERA5_CDS_API_KEY not configured"}
else
case submit_request(dataset, request, api_key) do
{:ok, job_id} -> poll_and_download(job_id, api_key)
{:error, reason} -> {:error, reason}
end
end
end
defp submit_request(dataset, request, api_key) do
url = "#{@cds_url}/retrieve/v1/processes/#{dataset}/execution"
case Req.post(url,
json: %{"inputs" => request},
headers: [{"PRIVATE-TOKEN", api_key}],
receive_timeout: 30_000
) do
{:ok, %{status: status, body: body}} when status in [200, 201, 202] ->
job_id = body["jobID"] || body["request_id"] || body["jobId"]
if job_id do
Logger.info("ERA5: submitted job #{job_id} for #{dataset}")
{:ok, job_id}
else
{:error, "ERA5: no job_id in response: #{inspect(body)}"}
end
{:ok, %{status: status, body: body}} ->
{:error, "ERA5 CDS HTTP #{status}: #{inspect(body)}"}
{:error, reason} ->
{:error, "ERA5 CDS request failed: #{inspect(reason)}"}
end
end
defp poll_and_download(job_id, api_key, attempt \\ 0)
defp poll_and_download(_job_id, _api_key, attempt) when attempt >= @max_poll_attempts do
{:error, "ERA5: poll timeout after #{@max_poll_attempts} attempts"}
end
defp poll_and_download(job_id, api_key, attempt) do
url = "#{@cds_url}/retrieve/v1/jobs/#{job_id}"
case Req.get(url, headers: [{"PRIVATE-TOKEN", api_key}], receive_timeout: 10_000) do
{:ok, %{status: 200, body: %{"status" => "successful"}}} ->
results_url = "#{url}/results"
download_result(results_url, api_key)
{:ok, %{status: 200, body: %{"status" => status}}}
when status in ["accepted", "running"] ->
Process.sleep(@poll_interval_ms)
poll_and_download(job_id, api_key, attempt + 1)
{:ok, %{status: 200, body: %{"status" => "failed"} = body}} ->
# CDS often omits "message" — fall back to dumping the whole body so the
# failure isn't summarised as `nil` in oban_jobs.errors.
reason = body["message"] || body["error"] || inspect(body)
{:error, "ERA5 job #{job_id} failed: #{reason}"}
{:ok, %{status: status, body: body}} ->
{:error, "ERA5 poll HTTP #{status}: #{inspect(body)}"}
{:error, reason} ->
{:error, "ERA5 poll failed: #{inspect(reason)}"}
end
end
defp download_result(url, api_key) do
case Req.get(url, headers: [{"PRIVATE-TOKEN", api_key}], receive_timeout: 120_000) do
{:ok, %{status: 200, body: %{"asset" => %{"value" => %{"href" => href}}}}} ->
# New CDS returns a JSON with download link
download_file(href, api_key)
{:ok, %{status: 200, headers: headers, body: body}} ->
# Direct binary download
content_type = %Req.Response{headers: headers} |> Req.Response.get_header("content-type") |> List.first("")
if String.contains?(content_type, "grib") or is_binary(body) do
{:ok, body}
else
{:error, "ERA5: unexpected response format"}
end
{:ok, %{status: status, body: body}} ->
{:error, "ERA5 download HTTP #{status}: #{inspect(body)}"}
{:error, reason} ->
{:error, "ERA5 download failed: #{inspect(reason)}"}
end
end
defp download_file(href, api_key) do
case Req.get(href, headers: [{"PRIVATE-TOKEN", api_key}], receive_timeout: 120_000) do
{:ok, %{status: 200, body: body}} -> {:ok, body}
{:ok, %{status: status}} -> {:error, "ERA5 file download HTTP #{status}"}
{:error, reason} -> {:error, "ERA5 file download failed: #{inspect(reason)}"}
end
end
defp build_profile(lat, lon, valid_time, single_grib, pressure_grib) do
# Decode GRIB2 data using existing infrastructure
with {:ok, single_points} <- extract_point(single_grib, lat, lon),
{:ok, pressure_points} <- extract_point(pressure_grib, lat, lon) do
# ERA5 temperatures are in Kelvin, convert to Celsius
temp_c = kelvin_to_celsius(single_points["TMP:2 m above ground"])
dewpoint_c = kelvin_to_celsius(single_points["DPT:2 m above ground"])
# ERA5 surface pressure is in Pa, convert to mb
pressure_mb = pa_to_mb(single_points["PRES:surface"])
hpbl_m = single_points["HPBL:surface"]
# ERA5 TCWV is in kg/m² which equals mm
pwat_mm = single_points["TCWV:entire atmosphere"]
# Build pressure-level profile
profile =
@pressure_levels
|> Enum.map(fn level_str ->
level = String.to_integer(level_str)
pres_key = "#{level} mb"
%{
"pres" => level * 1.0,
"tmpc" => kelvin_to_celsius(pressure_points["TMP:#{pres_key}"]),
"dwpc" => kelvin_to_celsius(pressure_points["DPT:#{pres_key}"]),
"hght" => geopotential_to_height(pressure_points["HGT:#{pres_key}"])
}
end)
|> Enum.reject(fn p -> is_nil(p["tmpc"]) end)
# Derive refractivity and ducting from profile
params = SoundingParams.derive(profile)
attrs = %{
valid_time: valid_time,
lat: lat,
lon: lon,
profile: profile,
hpbl_m: hpbl_m,
pwat_mm: pwat_mm,
surface_temp_c: temp_c,
surface_dewpoint_c: dewpoint_c,
surface_pressure_mb: pressure_mb
}
attrs =
if params do
Map.merge(attrs, %{
surface_refractivity: params.surface_refractivity,
min_refractivity_gradient: params.min_refractivity_gradient,
ducting_detected: params.ducting_detected || false,
duct_characteristics: params.duct_characteristics
})
else
attrs
end
{:ok, attrs}
end
rescue
e -> {:error, "ERA5 GRIB2 decode failed: #{Exception.message(e)}"}
end
defp extract_point(grib_data, lat, lon) do
case Extractor.extract_points(grib_data, lat, lon) do
{:ok, fields} when map_size(fields) > 0 -> {:ok, fields}
{:ok, _} -> {:error, "no data at #{lat},#{lon}"}
{:error, reason} -> {:error, reason}
end
end
defp kelvin_to_celsius(nil), do: nil
defp kelvin_to_celsius(k), do: k - 273.15
defp pa_to_mb(nil), do: nil
defp pa_to_mb(pa), do: pa / 100.0
defp geopotential_to_height(nil), do: nil
defp geopotential_to_height(gp), do: gp / 9.80665
defp api_key do
System.get_env("ERA5_CDS_API_KEY")
end
end