prop/lib/microwaveprop/terrain/srtm.ex

defmodule Microwaveprop.Terrain.Srtm do
  @moduledoc false

  @samples 3601
  @void -32_768
  @base_url "https://elevation-tiles-prod.s3.amazonaws.com/skadi"

  @spec tile_filename(float(), float()) :: String.t()
  def tile_filename(lat, lon) do
    lat_floor = floor(lat)
    lon_floor = floor(lon)

    lat_prefix = if lat_floor >= 0, do: "N", else: "S"
    lon_prefix = if lon_floor >= 0, do: "E", else: "W"

    lat_str = lat_floor |> abs() |> Integer.to_string() |> String.pad_leading(2, "0")
    lon_str = lon_floor |> abs() |> Integer.to_string() |> String.pad_leading(3, "0")

    "#{lat_prefix}#{lat_str}#{lon_prefix}#{lon_str}.hgt"
  end

  @spec download_tile(float(), float(), String.t()) :: {:ok, String.t()} | {:error, term()}
  def download_tile(lat, lon, tiles_dir) do
    filename = tile_filename(lat, lon)
    lat_dir = String.slice(filename, 0, 3)
    url = "#{@base_url}/#{lat_dir}/#{filename}.gz"
    path = Path.join(tiles_dir, filename)

    case Req.get(url, req_options()) do
      {:ok, %{status: 200, body: body}} ->
        decompressed = :zlib.gunzip(body)
        File.write!(path, decompressed)
        {:ok, path}

      {:ok, %{status: 404}} ->
        {:error, :not_available}

      {:ok, %{status: status}} ->
        {:error, "SRTM download HTTP #{status}"}

      {:error, reason} ->
        {:error, "SRTM download error: #{inspect(reason)}"}
    end
  end

  @spec lookup(float(), float(), String.t()) ::
          {:ok, integer()} | {:error, :no_tile} | {:error, :void}
  def lookup(lat, lon, tiles_dir) do
    path = Path.join(tiles_dir, tile_filename(lat, lon))

    case :file.open(path, [:read, :binary, :raw]) do
      {:ok, fd} ->
        read_elevation(fd, lat, lon)

      {:error, :enoent} ->
        if File.dir?(tiles_dir) do
          case download_tile(lat, lon, tiles_dir) do
            {:ok, _path} ->
              case :file.open(path, [:read, :binary, :raw]) do
                {:ok, fd} -> read_elevation(fd, lat, lon)
                {:error, _} -> {:error, :no_tile}
              end

            {:error, _} ->
              {:error, :no_tile}
          end
        else
          {:error, :no_tile}
        end
    end
  end

  @spec fetch_elevation_profile(float(), float(), float(), float(), String.t(), pos_integer()) ::
          {:ok, list(map())} | {:error, term()}
  def fetch_elevation_profile(lat1, lon1, lat2, lon2, tiles_dir, n \\ 64) do
    pts = sample_path(lat1, lon1, lat2, lon2, n)
    dist_km = haversine_km(lat1, lon1, lat2, lon2)

    results =
      Enum.reduce_while(pts, {:ok, []}, fn pt, {:ok, acc} ->
        case lookup(pt.lat, pt.lon, tiles_dir) do
          {:ok, elev} ->
            entry = %{
              lat: pt.lat,
              lon: pt.lon,
              d: pt.d,
              elev: elev,
              dist_km: pt.d * dist_km
            }

            {:cont, {:ok, acc ++ [entry]}}

          {:error, reason} ->
            {:halt, {:error, reason}}
        end
      end)

    results
  end

  defp read_elevation(fd, lat, lon) do
    row = round((floor(lat) + 1 - lat) * (@samples - 1))
    col = round((lon - floor(lon)) * (@samples - 1))
    offset = (row * @samples + col) * 2

    result =
      case :file.pread(fd, offset, 2) do
        {:ok, <<elev::signed-big-integer-size(16)>>} when elev == @void ->
          {:error, :void}

        {:ok, <<elev::signed-big-integer-size(16)>>} ->
          {:ok, elev}

        _ ->
          {:error, :void}
      end

    :file.close(fd)
    result
  end

  defp req_options do
    defaults = [
      compressed: false,
      decode_body: false,
      retry: &retry?/2,
      max_retries: 3,
      retry_delay: &retry_delay/1
    ]

    overrides = Application.get_env(:microwaveprop, :srtm_req_options, [])
    Keyword.merge(defaults, overrides)
  end

  defp retry?(_request, response) do
    case response do
      %Req.Response{status: status} when status in [429, 500, 502, 503, 504] -> true
      %{__exception__: true} -> true
      _ -> false
    end
  end

  defp retry_delay(n) do
    base = Integer.pow(2, n) * 1_000
    jitter = :rand.uniform(1_000)
    base + jitter
  end

  defp sample_path(lat1, lon1, lat2, lon2, n) do
    for i <- 0..n do
      f = i / n

      %{
        lat: lat1 + f * (lat2 - lat1),
        lon: lon1 + f * (lon2 - lon1),
        d: f
      }
    end
  end

  defp haversine_km(lat1, lon1, lat2, lon2) do
    dlat = deg_to_rad(lat2 - lat1)
    dlon = deg_to_rad(lon2 - lon1)
    rlat1 = deg_to_rad(lat1)
    rlat2 = deg_to_rad(lat2)

    a =
      :math.sin(dlat / 2) ** 2 +
        :math.cos(rlat1) * :math.cos(rlat2) * :math.sin(dlon / 2) ** 2

    2 * 6371.0 * :math.asin(:math.sqrt(a))
  end

  defp deg_to_rad(deg), do: deg * :math.pi() / 180
end