towerops/lib/towerops/coverages/raster.ex

defmodule Towerops.Coverages.Raster do
  @moduledoc """
  Writes coverage outputs to disk: a Float32 GeoTIFF holding the raw
  RSSI grid (one band, dBm values, NaN for "no coverage"), and a
  cnHeat-style colored PNG for direct overlay on Leaflet.

  Both outputs are written to
  `priv/static/coverage/<organization_id>/<coverage_id>/` and served by
  the existing `Plug.Static` pipeline.

  Conversion uses GDAL command-line tools (already installed via
  `gdal-bin` in production and `brew install gdal` in dev).
  """

  alias Towerops.Coverages.Coverage

  require Logger

  @nan_sentinel 1.0e30

  # cnHeat-ish palette: green strong, yellow good, orange marginal, red poor.
  @palette """
  -50 0 200 0 230
  -65 200 230 0 230
  -75 255 165 0 220
  -85 220 60 60 200
  -95 100 0 0 160
  nv 0 0 0 0
  """

  @doc """
  Writes both rasters to a per-coverage directory under
  `priv/static/coverage/`. Returns `{:ok, %{tif: rel_path, png: rel_path,
  bbox: {min_lat, min_lon, max_lat, max_lon}}}`.

  `pixels` is a flat list of dBm values (or `:nan` for no-coverage)
  laid out row-major, north-to-south, west-to-east, with `ncols` per
  row and `nrows` rows total. `bbox` is `{min_lat, min_lon, max_lat,
  max_lon}` in WGS84.
  """
  @spec write(
          Coverage.t(),
          [number() | :nan],
          %{
            ncols: pos_integer(),
            nrows: pos_integer(),
            bbox: {number(), number(), number(), number()}
          },
          String.t()
        ) ::
          {:ok, %{tif: String.t(), png: String.t(), bbox: tuple()}} | {:error, term()}
  def write(coverage, pixels, dims, suffix \\ "")

  def write(%Coverage{} = coverage, pixels, %{ncols: ncols, nrows: nrows, bbox: bbox}, suffix)
      when is_list(pixels) and is_binary(suffix) do
    expected = ncols * nrows

    if length(pixels) == expected do
      do_write(coverage, pixels, ncols, nrows, bbox, suffix)
    else
      {:error, {:pixel_count_mismatch, %{expected: expected, got: length(pixels)}}}
    end
  end

  defp do_write(coverage, pixels, ncols, nrows, {min_lat, min_lon, max_lat, max_lon} = bbox, suffix) do
    out_dir = output_dir(coverage)
    File.mkdir_p!(out_dir)

    base = "rssi" <> suffix
    raw_path = Path.join(out_dir, base <> ".f32")
    vrt_path = Path.join(out_dir, base <> ".vrt")
    tif_path = Path.join(out_dir, base <> ".tif")
    png_path = Path.join(out_dir, base <> ".png")
    palette_path = Path.join(out_dir, "palette.txt")

    with :ok <- write_raw(raw_path, pixels),
         :ok <- write_vrt(vrt_path, raw_path, ncols, nrows, min_lon, max_lat, max_lon, min_lat),
         :ok <- run_gdal_translate(vrt_path, tif_path),
         :ok <- File.write(palette_path, @palette),
         :ok <- run_gdaldem(tif_path, palette_path, png_path) do
      _ = File.rm(raw_path)
      _ = File.rm(vrt_path)
      _ = File.rm(palette_path)

      {:ok,
       %{
         tif: relative_static_path(coverage, base <> ".tif"),
         png: relative_static_path(coverage, base <> ".png"),
         bbox: bbox
       }}
    end
  end

  @doc "Removes the coverage's output directory if it exists."
  @spec cleanup(Coverage.t()) :: :ok
  def cleanup(%Coverage{} = coverage) do
    out_dir = output_dir(coverage)
    _ = if File.dir?(out_dir), do: File.rm_rf!(out_dir)
    :ok
  end

  @doc """
  Reads the RSSI value (dBm) at a single WGS84 point from the
  coverage's GeoTIFF via `gdallocationinfo`.

  Returns:
    * `{:ok, dbm}` — pixel had a valid value
    * `{:ok, :no_coverage}` — pixel was the NoData sentinel
    * `{:error, :outside_bbox}` — point is outside the raster's extent
    * `{:error, :not_ready}` — coverage hasn't been computed yet
    * `{:error, term}` — gdal failure
  """
  @spec query_rssi(Coverage.t(), float(), float()) ::
          {:ok, float() | :no_coverage} | {:error, term()}
  def query_rssi(%Coverage{status: status, raster_path: nil}, _lat, _lon) when status != "ready", do: {:error, :not_ready}

  def query_rssi(%Coverage{raster_path: nil}, _lat, _lon), do: {:error, :not_ready}

  def query_rssi(%Coverage{} = coverage, lat, lon) when is_number(lat) and is_number(lon) do
    if inside_bbox?(coverage, lat, lon) do
      tif = absolute_raster_path(coverage)

      if File.exists?(tif) do
        run_gdallocationinfo(tif, lon, lat)
      else
        {:error, :raster_missing}
      end
    else
      {:error, :outside_bbox}
    end
  end

  defp inside_bbox?(
         %Coverage{bbox_min_lat: min_lat, bbox_max_lat: max_lat, bbox_min_lon: min_lon, bbox_max_lon: max_lon},
         lat,
         lon
       )
       when is_number(min_lat) and is_number(max_lat) and is_number(min_lon) and is_number(max_lon) do
    lat >= min_lat and lat <= max_lat and lon >= min_lon and lon <= max_lon
  end

  defp inside_bbox?(_coverage, _lat, _lon), do: false

  defp absolute_raster_path(%Coverage{raster_path: rel}) do
    # rel looks like "/coverage/<org>/<id>/rssi*.tif"; strip the URL prefix
    # and rejoin against the configured storage dir (which already corresponds
    # to /coverage in the static plug).
    suffix =
      rel
      |> Path.relative_to("/")
      |> String.replace_prefix("coverage/", "")

    Path.join(coverage_dir(), suffix)
  end

  # gdallocationinfo with -wgs84 -valonly prints just the pixel value or
  # a blank line if outside, plus an error code line.
  defp run_gdallocationinfo(tif, lon, lat) do
    args = ["-valonly", "-wgs84", tif, Float.to_string(lon * 1.0), Float.to_string(lat * 1.0)]

    case System.cmd("gdallocationinfo", args, stderr_to_stdout: true) do
      {output, 0} ->
        parse_gdal_value(output)

      {output, code} ->
        {:error, {:gdallocationinfo_failed, code, String.slice(output, 0, 200)}}
    end
  end

  defp parse_gdal_value(output) do
    trimmed = output |> String.trim() |> String.split("\n") |> List.first() |> to_string()

    cond do
      trimmed == "" ->
        {:error, :outside_bbox}

      String.contains?(trimmed, "outside") ->
        {:error, :outside_bbox}

      true ->
        case Float.parse(trimmed) do
          {v, _} when v >= @nan_sentinel * 0.99 -> {:ok, :no_coverage}
          {v, _} -> {:ok, v}
          :error -> {:error, {:bad_gdal_output, trimmed}}
        end
    end
  end

  defp output_dir(%Coverage{organization_id: org_id, id: id}) do
    Path.join([coverage_dir(), to_string(org_id), to_string(id)])
  end

  # Resolves :coverage_storage_dir config to an absolute path.
  #
  #   * `{:otp_app, "rel/path"}` → `Application.app_dir/2` (dev/test default
  #     so rasters live under priv/static and Mix releases pick them up).
  #   * `"/abs/path"` → used as-is (production points at the NFS mount).
  defp coverage_dir do
    case Application.get_env(:towerops, :coverage_storage_dir, {:towerops, "priv/static/coverage"}) do
      {app, rel} when is_atom(app) and is_binary(rel) -> Application.app_dir(app, rel)
      path when is_binary(path) -> path
    end
  end

  defp relative_static_path(coverage, filename) do
    Path.join([
      "/coverage",
      to_string(coverage.organization_id),
      to_string(coverage.id),
      filename
    ])
  end

  defp write_raw(path, pixels) do
    bin =
      pixels
      |> Enum.map(&pixel_to_float/1)
      |> Enum.map(fn f -> <<f::float-32-little>> end)
      |> IO.iodata_to_binary()

    File.write(path, bin)
  end

  defp pixel_to_float(:nan), do: @nan_sentinel
  defp pixel_to_float(n) when is_number(n), do: n * 1.0

  # Build a minimal VRT that wraps the raw Float32 file. We give it a
  # geotransform (north-up, square pixels in degrees) so gdal_translate
  # can write a properly georeferenced GeoTIFF.
  defp write_vrt(vrt_path, raw_path, ncols, nrows, min_lon, max_lat, max_lon, min_lat) do
    pixel_size_x = (max_lon - min_lon) / ncols
    pixel_size_y = (max_lat - min_lat) / nrows

    vrt = """
    <VRTDataset rasterXSize="#{ncols}" rasterYSize="#{nrows}">
      <SRS>EPSG:4326</SRS>
      <GeoTransform>#{min_lon}, #{pixel_size_x}, 0.0, #{max_lat}, 0.0, #{-pixel_size_y}</GeoTransform>
      <VRTRasterBand dataType="Float32" band="1" subClass="VRTRawRasterBand">
        <SourceFilename relativeToVRT="1">#{Path.basename(raw_path)}</SourceFilename>
        <ImageOffset>0</ImageOffset>
        <PixelOffset>4</PixelOffset>
        <LineOffset>#{ncols * 4}</LineOffset>
        <ByteOrder>LSB</ByteOrder>
        <NoDataValue>#{@nan_sentinel}</NoDataValue>
      </VRTRasterBand>
    </VRTDataset>
    """

    File.write(vrt_path, vrt)
  end

  defp run_gdal_translate(vrt_path, tif_path) do
    case System.cmd(
           "gdal_translate",
           [
             "-q",
             "-of",
             "GTiff",
             "-co",
             "COMPRESS=DEFLATE",
             "-co",
             "TILED=YES",
             "-a_nodata",
             "#{@nan_sentinel}",
             vrt_path,
             tif_path
           ],
           stderr_to_stdout: true
         ) do
      {_, 0} -> :ok
      {output, code} -> {:error, {:gdal_translate_failed, code, output}}
    end
  end

  defp run_gdaldem(tif_path, palette_path, png_path) do
    case System.cmd(
           "gdaldem",
           [
             "color-relief",
             "-of",
             "PNG",
             "-alpha",
             "-nearest_color_entry",
             tif_path,
             palette_path,
             png_path
           ],
           stderr_to_stdout: true
         ) do
      {_, 0} -> :ok
      {output, code} -> {:error, {:gdaldem_failed, code, output}}
    end
  end
end