diff --git a/lib/microwaveprop/weather/era5_batch_client.ex b/lib/microwaveprop/weather/era5_batch_client.ex index e1f6853a..2bd81be4 100644 --- a/lib/microwaveprop/weather/era5_batch_client.ex +++ b/lib/microwaveprop/weather/era5_batch_client.ex @@ -337,9 +337,13 @@ defmodule Microwaveprop.Weather.Era5BatchClient do # ERA5 pressure-level moisture is specific humidity (kg/kg). Convert to # dewpoint in °C via the Magnus-Tetens inverse already used for HRRR native # profiles. `level_pa` is the pressure level itself in Pa. + # ThetaE.dewpoint_from_spfh/2 returns Kelvin. Profile `dwpc` entries + # are Celsius (the downstream SoundingParams.compute_refractivity_profile + # passes p["dwpc"] straight to the Buck saturation-vapor-pressure + # equation which expects °C). Convert here. defp dewpoint_from_spfh(nil, _level_pa), do: nil defp dewpoint_from_spfh(spfh, _level_pa) when spfh <= 0, do: nil - defp dewpoint_from_spfh(spfh, level_pa), do: ThetaE.dewpoint_from_spfh(spfh, level_pa) + defp dewpoint_from_spfh(spfh, level_pa), do: ThetaE.dewpoint_from_spfh(spfh, level_pa) - 273.15 # -- Bulk insert ----------------------------------------------------------- diff --git a/lib/microwaveprop/weather/narr_client.ex b/lib/microwaveprop/weather/narr_client.ex index a18c9be8..33822458 100644 --- a/lib/microwaveprop/weather/narr_client.ex +++ b/lib/microwaveprop/weather/narr_client.ex @@ -9,14 +9,44 @@ defmodule Microwaveprop.Weather.NarrClient do architecture — the filename is a historical artifact from an earlier MERRA-2 plan. - Currently provides URL/time helpers, the inventory parser, and the - inventory HTTP fetch. The byte-range `fetch_profile_at/2` pipeline - lands in a follow-up task of the same plan. + Provides URL/time helpers, the inventory parser, the inventory HTTP + fetch, the cdo-based composite-file profile extractor, and the + end-to-end byte-range `fetch_profile_at/2` pipeline. """ + alias Microwaveprop.Weather.SoundingParams + alias Microwaveprop.Weather.ThetaE + @base_url "https://www.ncei.noaa.gov/data/north-american-regional-reanalysis/access/3-hourly" @valid_hours [0, 3, 6, 9, 12, 15, 18, 21] + # Surface records we need from the inventory. Tuples of {var, level} keyed + # exactly as `parse_inventory/2` returns them. Pressure-level records are + # added dynamically below from whatever (TMP|HGT|SPFH):"NN mb" tuples the + # parsed inventory contains. + @surface_records [ + {"HPBL", "sfc"}, + {"PRES", "sfc"}, + {"TMP", "2 m above gnd"}, + {"DPT", "2 m above gnd"}, + {"PWAT", "atmos col"} + ] + + # cdo emits NCEP GRIB1 parameters as `varNNN`. This table maps the cdo + # parameter token to the semantic NCEP short name, verified empirically + # in the spike (see plan doc cdo lookup table). + @cdo_var_to_name %{ + "var1" => "PRES", + "var7" => "HGT", + "var11" => "TMP", + "var17" => "DPT", + "var33" => "UGRD", + "var34" => "VGRD", + "var51" => "SPFH", + "var54" => "PWAT", + "var221" => "HPBL" + } + @doc """ Builds the NCEI HTTPS URL for the NARR analysis GRIB1 file at `valid_time`. @@ -110,6 +140,252 @@ defmodule Microwaveprop.Weather.NarrClient do end end + @doc """ + Extracts an `era5_profiles`-shaped attribute map from a composite NARR + GRIB1 file at the given lat/lon point. + + Pure-ish: shells out to `cdo -outputtab,name,lev,value -remapnn,lon=…_lat=…` + on the local file. No network. + + Returns `{:ok, attrs}` where `attrs` has the same shape + `Era5BatchClient.build_profile_attrs/5` produces, or `{:error, reason}` + on missing file, cdo failure, or empty parser output. + """ + @spec extract_profile_from_file(Path.t(), float(), float()) :: + {:ok, map()} | {:error, term()} + def extract_profile_from_file(grb_path, lat, lon) when is_binary(grb_path) and is_number(lat) and is_number(lon) do + if File.exists?(grb_path) do + run_cdo_outputtab(grb_path, lat, lon) + else + {:error, "NARR cdo extract: file not found: #{grb_path}"} + end + end + + @doc """ + Fetches all the NARR records we need for `valid_time` via byte-range + `Range:` GETs against the NCEI grib file, merges them into a composite + GRIB1 with `cdo -merge`, then extracts the profile at `{lat, lon}` + using `extract_profile_from_file/3`. + + Returns `{:ok, profile_attrs}` shaped like `Era5BatchClient.build_profile_attrs/5` + or `{:error, reason}`. Temp files are always cleaned up via `try/after`. + """ + @spec fetch_profile_at(DateTime.t(), {float(), float()}) :: + {:ok, map()} | {:error, term()} + def fetch_profile_at(%DateTime{} = valid_time, {lat, lon}) when is_number(lat) and is_number(lon) do + grb_url = url_for(valid_time) + + with {:ok, index} <- fetch_inventory(valid_time), + {:ok, record_keys} <- pick_records(index) do + do_fetch_profile(grb_url, index, record_keys, lat, lon) + end + end + + defp do_fetch_profile(grb_url, index, record_keys, lat, lon) do + tag = "narr_#{System.unique_integer([:positive])}" + tmp_dir = System.tmp_dir!() + merged_path = Path.join(tmp_dir, "#{tag}_merged.grb") + + record_paths = + Enum.with_index(record_keys, fn key, i -> + {key, Path.join(tmp_dir, "#{tag}_#{i}.grb")} + end) + + try do + with :ok <- fetch_records(grb_url, index, record_paths), + :ok <- merge_records(record_paths, merged_path) do + extract_profile_from_file(merged_path, lat, lon) + end + after + Enum.each(record_paths, fn {_key, path} -> File.rm(path) end) + File.rm(merged_path) + end + end + + # Builds the list of {var, level} keys we want to pull from the inventory: + # the fixed surface set plus every (TMP|HGT|SPFH) at "NN mb" present in + # the index. Returns `{:error, _}` if any of the required surface records + # are missing from the inventory. + defp pick_records(index) do + missing_surface = Enum.reject(@surface_records, &Map.has_key?(index, &1)) + + if missing_surface == [] do + pressure_keys = + index + |> Map.keys() + |> Enum.filter(fn + {var, level} when var in ["TMP", "HGT", "SPFH"] -> + String.ends_with?(level, " mb") + + _ -> + false + end) + + {:ok, @surface_records ++ pressure_keys} + else + {:error, "NARR inventory missing required surface records: #{inspect(missing_surface)}"} + end + end + + defp fetch_records(grb_url, index, record_paths) do + Enum.reduce_while(record_paths, :ok, fn {key, path}, _acc -> + {offset, length} = Map.fetch!(index, key) + + case fetch_byte_range(grb_url, offset, length, path) do + :ok -> {:cont, :ok} + {:error, reason} -> {:halt, {:error, reason}} + end + end) + end + + defp fetch_byte_range(url, offset, length, dest_path) do + range_header = "bytes=#{offset}-#{offset + length - 1}" + + case Req.get(url, [headers: [{"range", range_header}], receive_timeout: 60_000] ++ req_options()) do + {:ok, %{status: status, body: body}} when status in [200, 206] and is_binary(body) -> + File.write(dest_path, body) + + {:ok, %{status: status}} -> + {:error, "NARR byte-range HTTP #{status} for #{range_header}"} + + {:error, reason} -> + {:error, "NARR byte-range request failed: #{inspect(reason)}"} + end + end + + defp merge_records(record_paths, merged_path) do + inputs = Enum.map(record_paths, fn {_key, path} -> path end) + args = ["-merge"] ++ inputs ++ [merged_path] + + case System.cmd("cdo", args, stderr_to_stdout: true) do + {_output, 0} -> :ok + {output, code} -> {:error, "cdo -merge exited #{code}: #{output}"} + end + end + + defp run_cdo_outputtab(grb_path, lat, lon) do + args = [ + "-outputtab,name,lev,value", + "-remapnn,lon=#{lon}_lat=#{lat}", + grb_path + ] + + case System.cmd("cdo", args, stderr_to_stdout: false) do + {stdout, 0} -> + parse_cdo_outputtab(stdout) + + {output, code} -> + {:error, "cdo extract exited #{code}: #{output}"} + end + end + + defp parse_cdo_outputtab(stdout) do + raw = + stdout + |> String.split("\n", trim: true) + |> Enum.map(&String.trim/1) + |> Enum.filter(&String.starts_with?(&1, "var")) + |> Enum.reduce(%{}, &accumulate_cdo_row/2) + + if raw == %{} do + {:error, "NARR cdo extract: no parseable values in cdo output"} + else + {:ok, build_profile_attrs(raw)} + end + end + + defp accumulate_cdo_row(line, acc) do + case parse_cdo_row(line) do + {key, value} -> Map.put(acc, key, value) + :error -> acc + end + end + + defp parse_cdo_row(line) do + with [var_token, lev_str, value_str] <- String.split(line, ~r/\s+/, trim: true), + {:ok, name} <- Map.fetch(@cdo_var_to_name, var_token), + {lev_int, ""} <- Integer.parse(lev_str), + {value_float, ""} <- Float.parse(value_str) do + {{name, lev_int}, value_float} + else + _ -> :error + end + end + + defp build_profile_attrs(raw) do + profile = build_pressure_profile(raw) + + base = %{ + profile: profile, + surface_temp_c: kelvin_to_celsius(raw[{"TMP", 2}]), + surface_dewpoint_c: kelvin_to_celsius(raw[{"DPT", 2}]), + surface_pressure_mb: pa_to_mb(raw[{"PRES", 0}]), + hpbl_m: raw[{"HPBL", 0}], + pwat_mm: raw[{"PWAT", 0}] + } + + params = SoundingParams.derive(profile) + + if params do + Map.merge(base, %{ + surface_refractivity: params.surface_refractivity, + min_refractivity_gradient: params.min_refractivity_gradient, + ducting_detected: params.ducting_detected || false, + duct_characteristics: params.duct_characteristics + }) + else + Map.merge(base, %{ + surface_refractivity: nil, + min_refractivity_gradient: nil, + ducting_detected: false, + duct_characteristics: nil + }) + end + end + + defp build_pressure_profile(raw) do + raw + |> Map.keys() + |> Enum.filter(fn + {"TMP", lev} when lev > 100 -> true + _ -> false + end) + |> Enum.map(fn {"TMP", lev} -> lev end) + |> Enum.sort(:desc) + |> Enum.map(&pressure_level_entry(&1, raw)) + |> Enum.reject(&is_nil/1) + end + + defp pressure_level_entry(lev_pa, raw) do + with t_k when is_float(t_k) <- raw[{"TMP", lev_pa}], + h_m when is_float(h_m) <- raw[{"HGT", lev_pa}], + spfh when is_float(spfh) <- raw[{"SPFH", lev_pa}] do + %{ + "pres" => lev_pa / 100.0, + "tmpc" => kelvin_to_celsius(t_k), + "hght" => h_m, + "dwpc" => dewpoint_from_spfh(spfh, lev_pa) + } + else + _ -> nil + end + end + + # ThetaE.dewpoint_from_spfh/2 returns Kelvin (the function adds 273.15 + # at the end of the Magnus-Tetens inversion). Profile entries store + # `dwpc` in *Celsius* — the `c` is load-bearing — and downstream + # SoundingParams.compute_refractivity_profile passes p["dwpc"] straight + # to the Buck saturation-vapor-pressure function which expects °C. So + # convert here. + defp dewpoint_from_spfh(spfh, _lev_pa) when not is_number(spfh) or spfh <= 0, do: nil + defp dewpoint_from_spfh(spfh, lev_pa), do: ThetaE.dewpoint_from_spfh(spfh, lev_pa * 1.0) - 273.15 + + defp kelvin_to_celsius(nil), do: nil + defp kelvin_to_celsius(k) when is_number(k), do: k - 273.15 + + defp pa_to_mb(nil), do: nil + defp pa_to_mb(pa) when is_number(pa), do: pa / 100.0 + defp do_get_inventory(url) do case Req.get(url, [receive_timeout: 30_000] ++ req_options()) do {:ok, %{status: 200, body: body}} when is_binary(body) -> diff --git a/test/fixtures/narr/narr_dfw_2010-06-15_12z.grb b/test/fixtures/narr/narr_dfw_2010-06-15_12z.grb new file mode 100644 index 00000000..d20eb656 Binary files /dev/null and b/test/fixtures/narr/narr_dfw_2010-06-15_12z.grb differ diff --git a/test/microwaveprop/weather/narr_client_test.exs b/test/microwaveprop/weather/narr_client_test.exs index 832eb59a..3eaf654e 100644 --- a/test/microwaveprop/weather/narr_client_test.exs +++ b/test/microwaveprop/weather/narr_client_test.exs @@ -182,4 +182,96 @@ defmodule Microwaveprop.Weather.NarrClientTest do assert {:error, _reason} = NarrClient.fetch_inventory(~U[2010-06-15 12:00:00Z]) end end + + describe "extract_profile_from_file/3" do + @grb_fixture "test/fixtures/narr/narr_dfw_2010-06-15_12z.grb" + + test "extracts a profile_attrs map from the spike-captured composite GRIB1 fixture" do + assert {:ok, attrs} = + NarrClient.extract_profile_from_file(@grb_fixture, 32.9, -97.0) + + # Surface vars (verified in the spike against real DFW 2010-06-15 12Z analysis) + # 296.94 K + assert_in_delta attrs.surface_temp_c, 23.79, 0.5 + # 294.49 K + assert_in_delta attrs.surface_dewpoint_c, 21.34, 0.5 + # 99383.6 Pa + assert_in_delta attrs.surface_pressure_mb, 993.84, 1.0 + assert_in_delta attrs.hpbl_m, 703.5, 5.0 + assert_in_delta attrs.pwat_mm, 45.1, 1.0 + + # Profile array — fixture only contains 850 mb, so length is 1 + assert length(attrs.profile) == 1 + [level_850] = attrs.profile + assert_in_delta level_850["pres"], 850.0, 0.1 + # 291.52 K + assert_in_delta level_850["tmpc"], 18.37, 0.5 + assert_in_delta level_850["hght"], 1547.8, 5.0 + # Derived from SPFH=0.01040 kg/kg at 850 mb via Magnus-Tetens. + # Hand-check: e ≈ q*P/(0.622+0.378*q) ≈ 1410 Pa = 14.10 hPa, + # ln(14.10/6.1078) ≈ 0.836, Td_C = 243.04 * 0.836 / (17.625 - 0.836) ≈ 12.1 °C. + # Stored as Celsius — Kelvin would be ~285 and would crash the + # downstream Buck equation in SoundingParams. + assert_in_delta level_850["dwpc"], 12.1, 1.0 + + # Derived fields — SoundingParams.derive returns these (may be nil if profile is too short) + assert Map.has_key?(attrs, :surface_refractivity) + assert Map.has_key?(attrs, :min_refractivity_gradient) + assert Map.has_key?(attrs, :ducting_detected) + assert Map.has_key?(attrs, :duct_characteristics) + end + + test "returns {:error, _} when the file doesn't exist" do + assert {:error, _} = + NarrClient.extract_profile_from_file("/tmp/does_not_exist.grb", 32.9, -97.0) + end + end + + describe "fetch_profile_at/2" do + @inv_fixture "test/fixtures/narr/narr-a_221_20100615_1200_000.inv" + @grb_fixture "test/fixtures/narr/narr_dfw_2010-06-15_12z.grb" + @file_size 56_221_430 + + test "byte-range fetches records, merges via cdo, and returns the profile_attrs" do + inv_body = File.read!(@inv_fixture) + grb_body = File.read!(@grb_fixture) + + Req.Test.stub(NarrClient, fn conn -> + cond do + String.ends_with?(conn.request_path, ".inv") and conn.method == "GET" -> + Plug.Conn.send_resp(conn, 200, inv_body) + + String.ends_with?(conn.request_path, ".grb") and conn.method == "HEAD" -> + conn + |> Plug.Conn.put_resp_header("content-length", Integer.to_string(@file_size)) + |> Plug.Conn.send_resp(200, "") + + String.ends_with?(conn.request_path, ".grb") and conn.method == "GET" -> + # Byte-range request — return the entire fixture body. The composite + # fixture happens to contain ALL the records the byte-range fetcher + # asks for (just at different offsets than the real prod file), so + # cdo -merge will produce the same composite either way for the + # purpose of this test. The test asserts the OUTPUT of extract, + # not the byte-range mechanics. + Plug.Conn.send_resp(conn, 200, grb_body) + + true -> + Plug.Conn.send_resp(conn, 500, "unexpected #{conn.method} #{conn.request_path}") + end + end) + + assert {:ok, attrs} = + NarrClient.fetch_profile_at(~U[2010-06-15 12:00:00Z], {32.9, -97.0}) + + assert_in_delta attrs.surface_temp_c, 23.79, 0.5 + end + + test "returns {:error, _} when fetch_inventory fails" do + Req.Test.stub(NarrClient, fn conn -> + Plug.Conn.send_resp(conn, 404, "not found") + end) + + assert {:error, _} = NarrClient.fetch_profile_at(~U[2010-06-15 12:00:00Z], {32.9, -97.0}) + end + end end