prop/lib/microwaveprop/propagation/grid.ex

defmodule Microwaveprop.Propagation.Grid do
  @moduledoc """
  Propagation scoring grid definitions at 0.125° (~14 km).

  Two regions, disjoint by construction:

  * `conus_points/0` — HRRR-sourced CONUS cells (lat 25-50, lon -125 to -66).
  * `hrdps_only_points/0` — Canadian cells covered by HRDPS but **not** by
    HRRR. Bounded north at 60°N because SRTM elevation data stops at 60°N
    (Stage 8 of the HRDPS plan defers Arctic coverage). The mask is a
    rectangular bbox minus the HRRR overlap; over-water cells (Hudson Bay,
    Great Lakes spillover) get scored and ignored downstream.

  `point_source/1` returns `:hrrr | :hrdps | :neither` so callers can
  route fetch/score work to the correct NWP source.
  """

  @lat_min 25.0
  @lat_max 50.0
  @lon_min -125.0
  @lon_max -66.0
  @step 0.125

  @hrdps_lat_min 49.0
  @hrdps_lat_max 60.0
  @hrdps_lon_min -141.0
  @hrdps_lon_max -52.0

  @conus_points for(
                  lat <-
                    Enum.map(
                      0..round((@lat_max - @lat_min) / @step),
                      fn i -> Float.round(@lat_min + i * @step, 3) end
                    ),
                  lon <-
                    Enum.map(
                      0..round((@lon_max - @lon_min) / @step),
                      fn i -> Float.round(@lon_min + i * @step, 3) end
                    ),
                  do: {lat, lon}
                )

  @doc "Returns all grid points as `{lat, lon}` tuples covering CONUS at 0.125 degree spacing."
  @spec conus_points() :: [{float(), float()}]
  def conus_points do
    @conus_points
  end

  @doc "Returns the grid step size in degrees."
  @spec step() :: float()
  def step, do: @step

  @doc "Returns the CONUS bounding box as a map."
  @spec bounds() :: %{lat_min: float(), lat_max: float(), lon_min: float(), lon_max: float()}
  def bounds, do: %{lat_min: @lat_min, lat_max: @lat_max, lon_min: @lon_min, lon_max: @lon_max}

  @doc """
  True when the given position lies inside the CONUS bounding box
  (inclusive on all four edges). Callers use this to gate HRRR
  enrichment — contacts whose path origin is outside the grid would
  be enqueued forever because `hrrr_point_rs` silently writes no
  profiles for out-of-grid points.

  Accepts the canonical `%{"lat" => _, "lon" => _}` pos map. Returns
  `false` for nil, missing keys, or nil numeric values.
  """
  @spec contains?(map() | nil) :: boolean()
  def contains?(nil), do: false

  def contains?(%{"lat" => lat, "lon" => lon}) when is_number(lat) and is_number(lon) do
    lat >= @lat_min and lat <= @lat_max and lon >= @lon_min and lon <= @lon_max
  end

  def contains?(%{lat: lat, lon: lon}) when is_number(lat) and is_number(lon) do
    lat >= @lat_min and lat <= @lat_max and lon >= @lon_min and lon <= @lon_max
  end

  def contains?(_), do: false

  @doc "Returns the grid specification for wgrib2 -lola extraction."
  @spec wgrib2_grid_spec() :: %{
          lon_start: float(),
          lon_count: non_neg_integer(),
          lon_step: float(),
          lat_start: float(),
          lat_count: non_neg_integer(),
          lat_step: float()
        }
  def wgrib2_grid_spec do
    lon_count = round((@lon_max - @lon_min) / @step) + 1
    lat_count = round((@lat_max - @lat_min) / @step) + 1

    %{
      lon_start: @lon_min,
      lon_count: lon_count,
      lon_step: @step,
      lat_start: @lat_min,
      lat_count: lat_count,
      lat_step: @step
    }
  end

  @doc """
  Canadian-only grid points. Cells inside the HRDPS bbox (49-60°N,
  -141 to -52°W) but outside HRRR's CONUS bbox. Disjoint from
  `conus_points/0` by construction so the two grids never double-write
  the same `(lat, lon)`.
  """
  @spec hrdps_only_points() :: [{float(), float()}]
  def hrdps_only_points do
    for lat <- float_range(@hrdps_lat_min, @hrdps_lat_max, @step),
        lon <- float_range(@hrdps_lon_min, @hrdps_lon_max, @step),
        not in_conus_bbox?(lat, lon) do
      {Float.round(lat, 3), Float.round(lon, 3)}
    end
  end

  @doc """
  Returns the source NWP model for a `{lat, lon}` point:

    * `:hrrr` — inside the CONUS bbox; route to HrrrClient.
    * `:hrdps` — inside the Canadian extent (excluding HRRR overlap);
      route to HrdpsClient.
    * `:neither` — outside both. Callers should fall back to the
      IEMRE/RAOB enrichment paths.
  """
  @spec point_source({float(), float()} | %{lat: number(), lon: number()}) ::
          :hrrr | :hrdps | :neither
  def point_source({lat, lon}), do: classify(lat, lon)
  def point_source(%{lat: lat, lon: lon}), do: classify(lat, lon)

  defp classify(lat, lon) do
    cond do
      in_conus_bbox?(lat, lon) -> :hrrr
      in_hrdps_bbox?(lat, lon) -> :hrdps
      true -> :neither
    end
  end

  defp in_conus_bbox?(lat, lon) do
    lat >= @lat_min and lat <= @lat_max and lon >= @lon_min and lon <= @lon_max
  end

  defp in_hrdps_bbox?(lat, lon) do
    lat >= @hrdps_lat_min and lat <= @hrdps_lat_max and
      lon >= @hrdps_lon_min and lon <= @hrdps_lon_max
  end

  defp float_range(start, stop, step) do
    count = round((stop - start) / step) + 1
    # Avoid floating-point drift by rounding each value to 3 decimals.
    # The accumulation of start + i * step drifts beyond 0.001° after ~200
    # iterations, which can cause map lookups on snapped coordinates to miss.
    Enum.map(0..(count - 1), fn i -> Float.round(start + i * step, 3) end)
  end
end