prop/lib/microwaveprop/weather/narr_client.ex

defmodule Microwaveprop.Weather.NarrClient do
  @moduledoc """
  Client for NCEP North American Regional Reanalysis (NARR) data served
  anonymously from NCEI at `https://www.ncei.noaa.gov/data/north-american-regional-reanalysis/`.

  Replaces the broken ERA5/CDS historical backfill with anonymous NCEI
  fetches for pre-HRRR (1979-01-01 → 2014-10-02) contact enrichment. See
  `docs/plans/2026-04-15-merra2-historical-backfill.md` for the full
  architecture — the filename is a historical artifact from an earlier
  MERRA-2 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]

  # NCEI's NARR archive publishes 1979-01-01 through 2014-10-01 21 UTC
  # (the last 3-hourly analysis). Dispatching fetches outside that window
  # just 404s, so callers filter on `in_coverage?/1` first.
  @coverage_start ~U[1979-01-01 00:00:00Z]
  @coverage_end ~U[2014-10-02 00:00:00Z]

  # 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"}
  ]

  # Map NCEP GRIB1 parameter codes to our semantic short names.
  # Using `-outputtab,code,...` (instead of `name`) makes cdo emit the raw
  # numeric param code, so the same parser works across cdo versions that
  # would otherwise produce either `varNN` (Mac / cdo 2.6) or `codeNN` +
  # named shortnames like `tpag10` / `2tag2` / `saip` (Debian / cdo 2.5.1).
  @cdo_code_to_name %{
    1 => "PRES",
    7 => "HGT",
    11 => "TMP",
    17 => "DPT",
    33 => "UGRD",
    34 => "VGRD",
    51 => "SPFH",
    54 => "PWAT",
    221 => "HPBL"
  }

  # cdo 2.5.1 on Debian ignores `-outputtab,code` for records whose GRIB1
  # parameter table entry has a shortname, emitting the shortname string
  # instead of the numeric code. Mapping verified against prod discard
  # output: 2tag2=PRES (985 mb surface pressure), saip=DPT (292 K 2-m
  # dewpoint), tpag10=HGT (983 m at 900 hPa geopotential height). UGRD
  # and VGRD aren't yet observed in that output — add them when they
  # surface in a real discard.
  @cdo_shortname_to_code %{
    "2tag2" => 1,
    "saip" => 17,
    "tpag10" => 7
  }

  @doc """
  Builds the NCEI HTTPS URL for the NARR analysis GRIB1 file at `valid_time`.

  Raises `ArgumentError` if `valid_time` is not on a 3-hourly analysis
  boundary — NARR analyses exist only at 00, 03, 06, 09, 12, 15, 18, and
  21 UTC.
  """
  @spec url_for(DateTime.t()) :: String.t()
  def url_for(valid_time) do
    validate_analysis_time!(valid_time)

    yyyymm = Calendar.strftime(valid_time, "%Y%m")
    yyyymmdd = Calendar.strftime(valid_time, "%Y%m%d")
    hhmm = Calendar.strftime(valid_time, "%H%M")

    "#{@base_url}/#{yyyymm}/#{yyyymmdd}/narr-a_221_#{yyyymmdd}_#{hhmm}_000.grb"
  end

  @doc """
  Builds the NCEI HTTPS URL for the NARR analysis inventory (`.inv`) file
  at `valid_time`. Same validation rules as `url_for/1`.
  """
  @spec url_for_inventory(DateTime.t()) :: String.t()
  def url_for_inventory(valid_time) do
    valid_time
    |> url_for()
    |> String.replace_suffix(".grb", ".inv")
  end

  @doc """
  Rounds a `DateTime` *down* to the nearest 3-hourly NARR analysis slot.

  NARR analyses are produced at 00, 03, 06, 09, 12, 15, 18, and 21 UTC.
  Minutes, seconds, and microseconds are zeroed out; the hour becomes
  `(hour div 3) * 3`.
  """
  @spec snap_to_analysis_hour(DateTime.t()) :: DateTime.t()
  def snap_to_analysis_hour(%DateTime{} = datetime) do
    snapped_hour = div(datetime.hour, 3) * 3
    %{datetime | hour: snapped_hour, minute: 0, second: 0, microsecond: {0, 0}}
  end

  @doc """
  True when `datetime` falls inside NCEI's NARR 3-hourly archive window
  (1979-01-01 through 2014-10-01 21 UTC). Callers use this to avoid
  dispatching fetches that would just 404.
  """
  @spec in_coverage?(DateTime.t()) :: boolean()
  def in_coverage?(%DateTime{} = datetime) do
    DateTime.compare(datetime, @coverage_start) != :lt and
      DateTime.before?(datetime, @coverage_end)
  end

  @doc "Exclusive end of the NARR coverage window (first timestamp NOT covered)."
  @spec coverage_end() :: DateTime.t()
  def coverage_end, do: @coverage_end

  @doc """
  Parses the raw text of a NARR `.inv` file into a map keyed by
  `{var, level}` with values of `{byte_offset, length}`.

  Each record's length is computed as `next_offset - this_offset`. The
  final record's length uses `file_size` as its "next offset".

  Blank lines and trailing whitespace are ignored. `var` and `level` are
  returned as-is (no trimming) so levels like `"2 m above gnd"` round-trip
  unchanged.
  """
  @spec parse_inventory(String.t(), non_neg_integer()) ::
          {:ok, %{{String.t(), String.t()} => {non_neg_integer(), non_neg_integer()}}}
  def parse_inventory(inv_text, file_size) when is_binary(inv_text) and is_integer(file_size) do
    records =
      inv_text
      |> String.split("\n", trim: true)
      |> Enum.map(&parse_inventory_line/1)
      |> Enum.reject(&is_nil/1)

    index =
      records
      |> Enum.zip(Enum.drop(records, 1) ++ [{nil, file_size, nil, nil}])
      |> Map.new(fn {{_num, offset, var, level}, {_next_num, next_offset, _next_var, _next_level}} ->
        {{var, level}, {offset, next_offset - offset}}
      end)

    {:ok, index}
  end

  @doc """
  Fetches the `.inv` text and the `.grb` `Content-Length` for `valid_time`,
  then delegates to `parse_inventory/2`.

  Does one `GET` against the inventory URL and one `HEAD` against the GRIB
  URL. Returns `{:ok, index}` on success; `{:error, reason}` if either HTTP
  call is non-200, transport-errored, or missing a `content-length` header.
  """
  @spec fetch_inventory(DateTime.t()) ::
          {:ok, %{{String.t(), String.t()} => {non_neg_integer(), non_neg_integer()}}}
          | {:error, term()}
  def fetch_inventory(%DateTime{} = valid_time) do
    inv_url = url_for_inventory(valid_time)
    grb_url = url_for(valid_time)

    with {:ok, inv_body} <- do_get_inventory(inv_url),
         {:ok, file_size} <- do_head_grib_size(grb_url) do
      parse_inventory(inv_body, file_size)
    end
  end

  @doc """
  Extracts a NarrProfile-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}` — a map with the NarrProfile fields
  (`:surface_temp_c`, `:surface_dewpoint_c`, `:surface_pressure_mb`,
  `:hpbl_m`, `:pwat_mm`, `:profile` list plus the derived fields from
  `SoundingParams.derive/1`) — 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 for direct `Repo.insert_all/3`
  against `NarrProfile`, 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
    Microwaveprop.Instrument.span([:narr, :fetch_profile], %{lat: lat, lon: lon}, fn ->
      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)
  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, env: %{}, 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,code,lev,value",
      "-remapnn,lon=#{lon}_lat=#{lat}",
      grb_path
    ]

    case System.cmd("cdo", args, env: %{}, stderr_to_stdout: true) do
      {stdout, 0} ->
        parse_cdo_outputtab(stdout)

      {output, code} ->
        {:error, "cdo extract exited #{code}: #{output}"}
    end
  end

  @doc false
  @spec parse_cdo_outputtab(String.t()) :: {:ok, map()} | {:error, String.t()}
  def parse_cdo_outputtab(stdout) do
    raw =
      stdout
      |> String.split("\n", trim: true)
      |> Enum.reduce(%{}, &accumulate_cdo_row/2)

    if raw == %{} do
      # Include the first 500 chars of cdo output so the error surfaces
      # what cdo actually printed — crucial for diagnosing year-specific
      # parameter-table quirks in NARR GRIB1.
      snippet = stdout |> String.slice(0, 500) |> String.replace("\n", " | ")
      {:error, "NARR cdo extract: no parseable values in cdo output (got: #{snippet})"}
    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

  # Rows are three whitespace-separated columns: <code> <lev> <value>. The
  # <code> token is one of:
  #   * a bare integer (cdo 2.6 on Mac / dev)
  #   * `codeNN` with "code" prefix (cdo 2.5.1 on Debian / prod)
  #   * a GRIB1 shortname like `2tag2` / `saip` / `tpag10` (cdo 2.5.1 again,
  #     for params whose GRIB1 parameter-table entry has a shortname)
  # Header, cdi warnings, and cdo progress lines lack the three-column shape
  # and fall through to :error.
  defp parse_cdo_row(line) do
    with [code_str, lev_str, value_str] <- String.split(line, ~r/\s+/, trim: true),
         {:ok, code_int} <- normalize_cdo_code_token(code_str),
         {:ok, name} <- Map.fetch(@cdo_code_to_name, code_int),
         {lev_int, ""} <- Integer.parse(lev_str),
         {value_float, ""} <- Float.parse(value_str) do
      {{name, lev_int}, value_float}
    else
      _ -> :error
    end
  end

  # Resolves the first column of a cdo `-outputtab,code,...` row to its
  # numeric GRIB1 parameter code. Handles all three cdo output dialects
  # we've seen in the wild (see `parse_cdo_row/1` for context).
  defp normalize_cdo_code_token("code" <> rest) do
    case Integer.parse(rest) do
      {code_int, ""} -> {:ok, code_int}
      _ -> :error
    end
  end

  defp normalize_cdo_code_token(token) do
    case Integer.parse(token) do
      {code_int, ""} ->
        {:ok, code_int}

      _ ->
        case Map.fetch(@cdo_shortname_to_code, token) do
          {:ok, code_int} -> {:ok, code_int}
          :error -> :error
        end
    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) ->
        {:ok, body}

      {:ok, %{status: status, body: body}} ->
        {:error, "NARR inventory HTTP #{status}: #{inspect(body)}"}

      {:error, reason} ->
        {:error, "NARR inventory request failed: #{inspect(reason)}"}
    end
  end

  defp do_head_grib_size(url) do
    case Req.head(url, [receive_timeout: 30_000] ++ req_options()) do
      {:ok, %{status: 200} = response} ->
        case content_length(response) do
          {:ok, size} -> {:ok, size}
          :error -> {:error, "NARR grib HEAD missing content-length"}
        end

      {:ok, %{status: status}} ->
        {:error, "NARR grib HEAD HTTP #{status}"}

      {:error, reason} ->
        {:error, "NARR grib HEAD failed: #{inspect(reason)}"}
    end
  end

  defp content_length(%Req.Response{} = response) do
    case response |> Req.Response.get_header("content-length") |> List.first() do
      nil ->
        :error

      value when is_binary(value) ->
        case Integer.parse(value) do
          {size, ""} when size >= 0 -> {:ok, size}
          _ -> :error
        end
    end
  end

  defp parse_inventory_line(line) do
    case String.split(line, ":") do
      [num, offset, _date, var, level | _rest] ->
        with {num_int, ""} <- Integer.parse(num),
             {offset_int, ""} <- Integer.parse(offset) do
          {num_int, offset_int, var, level}
        else
          _ -> nil
        end

      _ ->
        nil
    end
  end

  defp req_options do
    Application.get_env(:microwaveprop, :narr_req_options, [])
  end

  defp validate_analysis_time!(%DateTime{hour: hour, minute: 0, second: 0} = _vt) when hour in @valid_hours, do: :ok

  defp validate_analysis_time!(vt) do
    raise ArgumentError,
          "NARR analyses are only produced at 00/03/06/09/12/15/18/21 UTC on the hour, " <>
            "got #{inspect(vt)}"
  end
end