towerops/lib/towerops/lidar.ex

defmodule Towerops.Lidar do
  @moduledoc """
  Catalog and on-demand retrieval of LIDAR-derived DEMs covering Texas.

  We do not mirror raster data. The catalog stores tile metadata and the
  public HTTPS URL of each Cloud-Optimized GeoTIFF (USGS 3DEP, with TNRIS
  as fallback). Elevation reads stream byte ranges directly from those
  public buckets via GDAL's `/vsicurl/` driver.
  """

  import Ecto.Query

  alias Towerops.Lidar.Reader
  alias Towerops.Lidar.Sources.UsgsNed
  alias Towerops.Lidar.Tile
  alias Towerops.Repo

  @typedoc "{west, south, east, north} in WGS84 degrees"
  @type bbox :: {number(), number(), number(), number()}

  @doc """
  Lists catalog tiles whose footprint contains the given point, ordered
  best-first (highest resolution, then newest year).

  Filters out tiles flagged as unavailable.
  """
  @spec list_tiles_for_point(number(), number()) :: [Tile.t()]
  def list_tiles_for_point(lat, lon) do
    point = %Geo.Point{coordinates: {lon, lat}, srid: 4326}

    Tile
    |> where([t], t.availability == "available")
    |> where([t], fragment("ST_Intersects(?, ?)", t.footprint, ^point))
    |> order_by([t], asc: t.resolution_m, desc: t.year)
    |> Repo.all()
  end

  @doc """
  Lists catalog tiles whose footprint intersects the given bbox, ordered
  best-first (highest resolution, then newest year).
  """
  @spec list_tiles_for_area(bbox()) :: [Tile.t()]
  def list_tiles_for_area({west, south, east, north}) do
    envelope =
      %Geo.Polygon{
        coordinates: [
          [
            {west, south},
            {east, south},
            {east, north},
            {west, north},
            {west, south}
          ]
        ],
        srid: 4326
      }

    Tile
    |> where([t], t.availability == "available")
    |> where([t], fragment("ST_Intersects(?, ?)", t.footprint, ^envelope))
    |> order_by([t], asc: t.resolution_m, desc: t.year)
    |> Repo.all()
  end

  @doc """
  Returns the elevation in meters at the given (lat, lon).

  Picks the best (highest-resolution, newest) tile that covers the point.
  Falls through to the next-best tile if the chosen one returns nodata.
  Returns `{:error, :no_tile}` if no catalog tile covers the point, or
  `{:error, :nodata}` if every covering tile lacks data at that pixel.
  """
  @spec get_elevation(number(), number()) ::
          {:ok, float()} | {:error, :no_tile | :nodata | {:gdal, integer(), String.t()}}
  def get_elevation(lat, lon) do
    case list_tiles_for_point(lat, lon) do
      [] -> get_elevation_from_ned(lat, lon)
      tiles -> read_first_value(tiles, lat, lon)
    end
  end

  defp get_elevation_from_ned(lat, lon) do
    cond do
      not ned_fallback_enabled?() ->
        {:error, :no_tile}

      tile = UsgsNed.tile_for_point(lat, lon) ->
        Reader.elevation_at(tile, lat, lon)

      true ->
        {:error, :no_tile}
    end
  end

  defp read_first_value([], _lat, _lon), do: {:error, :nodata}

  defp read_first_value([tile | rest], lat, lon) do
    case Reader.elevation_at(tile, lat, lon) do
      {:ok, value} -> {:ok, value}
      {:error, :nodata} -> read_first_value(rest, lat, lon)
      {:error, _} = err when rest == [] -> err
      {:error, _} -> read_first_value(rest, lat, lon)
    end
  end

  @max_grid_cells 4_000_000

  @doc """
  Returns a 2D elevation grid covering the bbox at the given cell size
  (in degrees, EPSG:4326).

  Mosaics across all intersecting catalog tiles, filling nodata cells from
  next-best tiles. Caps at #{@max_grid_cells} cells to bound memory.
  """
  @spec get_elevation_grid(bbox(), number()) ::
          {:ok, map()}
          | {:error, :no_tile | :nodata | :grid_too_large | {:gdal, integer(), String.t()}}
  def get_elevation_grid({west, south, east, north} = bbox, cell_size_deg) do
    ncols = ceil((east - west) / cell_size_deg)
    nrows = ceil((north - south) / cell_size_deg)

    if ncols * nrows > @max_grid_cells do
      {:error, :grid_too_large}
    else
      case list_tiles_for_area(bbox) do
        [] -> mosaic_from_ned(bbox, cell_size_deg)
        tiles -> mosaic(tiles, bbox, cell_size_deg)
      end
    end
  end

  # When the curated 3DEP/TNRIS catalog has no tile for the bbox,
  # fall back to USGS NED (1/3 arc-second / ~10 m) which covers all
  # of CONUS, Alaska, Hawai'i, and the Caribbean US territories.
  # Tests disable this so they can assert on the bare :no_tile path
  # without spinning up GDAL stubs for synthetic NED URLs.
  defp mosaic_from_ned(bbox, cell_size_deg) do
    if ned_fallback_enabled?() do
      case UsgsNed.tiles_for_bbox(bbox) do
        [] -> {:error, :no_tile}
        tiles -> mosaic(tiles, bbox, cell_size_deg)
      end
    else
      {:error, :no_tile}
    end
  end

  defp ned_fallback_enabled? do
    Application.get_env(:towerops, :lidar_ned_fallback, true)
  end

  defp mosaic(tiles, bbox, cell_size_deg) do
    {grids, last_error} =
      Enum.reduce(tiles, {[], nil}, fn tile, {acc, err} ->
        case Reader.elevation_grid(tile, bbox, cell_size_deg) do
          {:ok, grid} -> {[grid | acc], err}
          {:error, reason} -> {acc, reason}
        end
      end)

    case Enum.reverse(grids) do
      [] -> {:error, last_error || :no_data}
      [first | rest] -> finalize_mosaic(first, rest)
    end
  end

  defp finalize_mosaic(base, additional) do
    nodata = base.nodata_value

    merged_cells =
      Enum.reduce(additional, base.cells, fn next, current ->
        merge_grid_cells(current, next.cells, nodata, next.nodata_value)
      end)

    if all_nodata?(merged_cells, nodata) do
      {:error, :nodata}
    else
      {:ok, %{base | cells: merged_cells}}
    end
  end

  defp merge_grid_cells(current_rows, next_rows, base_nodata, next_nodata) do
    current_rows
    |> Enum.zip(next_rows)
    |> Enum.map(fn {current_row, next_row} ->
      current_row
      |> Enum.zip(next_row)
      |> Enum.map(fn
        {cur, next} when cur == base_nodata and next != next_nodata -> next
        {cur, _next} -> cur
      end)
    end)
  end

  defp all_nodata?(rows, nodata) do
    Enum.all?(rows, fn row -> Enum.all?(row, &(&1 == nodata)) end)
  end
end