prop/lib/microwaveprop/canopy/bulk_fetch.ex

defmodule Microwaveprop.Canopy.BulkFetch do
  @moduledoc """
  One-shot prefetch of Lang et al. 2020 10 m global canopy-height COGs
  (https://share.phys.ethz.ch/~pf/nlangdata/ETH_GlobalCanopyHeight_10m_2020_version1/3deg_cogs/),
  sliced into per-degree uint8 tiles compatible with `Microwaveprop.Canopy`.

  Run on prod (mix tasks don't work there) via:

      bin/microwaveprop rpc 'Microwaveprop.Canopy.BulkFetch.run(32.8, -97.0, 200)'

  Source COGs are 3° × 3° mosaics, ~150 MB each. Each is sliced into
  the 1° × 1° uint8 .canopy tiles the runtime expects in
  `/data/canopy/`. Requires `gdal_translate` in the runtime image.
  """

  alias Microwaveprop.Canopy

  require Logger

  @cog_base_url "https://share.phys.ethz.ch/~pf/nlangdata/ETH_GlobalCanopyHeight_10m_2020_version1/3deg_cogs"
  @samples 3600

  @doc """
  Downloads + slices canopy tiles covering a `radius_mi` circle around
  `(lat, lon)`. Skips per-degree tiles already on disk. Returns
  `%{downloaded, sliced, cached, failed}` counters.
  """
  @spec run(float(), float(), pos_integer(), keyword()) :: map()
  def run(lat, lon, radius_mi, opts \\ []) do
    cache_dir = Keyword.get(opts, :cache_dir, default_cache_dir())
    staging_dir = Keyword.get(opts, :staging_dir, default_staging_dir())
    File.mkdir_p!(cache_dir)
    File.mkdir_p!(staging_dir)

    lat_deg = radius_mi / 69.0
    lon_deg = radius_mi / (69.0 * max(:math.cos(lat * :math.pi() / 180.0), 0.1))

    bbox = %{
      south: lat - lat_deg,
      north: lat + lat_deg,
      west: lon - lon_deg,
      east: lon + lon_deg
    }

    degrees = degree_tiles(bbox)
    cog_tiles = Enum.uniq_by(degrees, &cog_tile_key/1)

    Logger.info(
      "canopy bulk fetch: lat=#{lat} lon=#{lon} radius_mi=#{radius_mi} " <>
        "degree_tiles=#{length(degrees)} cog_tiles=#{length(cog_tiles)}"
    )

    state = %{downloaded: 0, sliced: 0, cached: 0, failed: 0}

    state =
      Enum.reduce(cog_tiles, state, fn deg, acc ->
        process_cog_tile(deg, staging_dir, cache_dir, degrees, acc)
      end)

    Logger.info(
      "canopy bulk fetch DONE: downloaded=#{state.downloaded} sliced=#{state.sliced} " <>
        "cached=#{state.cached} failed=#{state.failed}"
    )

    Map.put(state, :cache_dir, cache_dir)
  end

  defp process_cog_tile({cog_lat, cog_lon}, staging_dir, cache_dir, all_degrees, acc) do
    cog_path = Path.join(staging_dir, cog_filename(cog_lat, cog_lon))

    case ensure_cog(cog_path, cog_lat, cog_lon) do
      {:ok, :downloaded} ->
        slice_cog(cog_path, cog_lat, cog_lon, cache_dir, all_degrees, %{acc | downloaded: acc.downloaded + 1})

      {:ok, :cached} ->
        slice_cog(cog_path, cog_lat, cog_lon, cache_dir, all_degrees, acc)

      {:error, reason} ->
        Logger.warning("canopy: skipping COG #{cog_filename(cog_lat, cog_lon)}: #{inspect(reason)}")
        %{acc | failed: acc.failed + 1}
    end
  end

  defp ensure_cog(path, cog_lat, cog_lon) do
    if File.exists?(path) and File.stat!(path).size > 1_000_000 do
      {:ok, :cached}
    else
      url = "#{@cog_base_url}/#{cog_filename(cog_lat, cog_lon)}"
      Logger.info("canopy: downloading COG #{cog_filename(cog_lat, cog_lon)}")

      case Req.get(url, into: File.stream!(path), retry: false, receive_timeout: 600_000) do
        {:ok, %{status: 200}} -> {:ok, :downloaded}
        {:ok, %{status: status}} -> {:error, {:http, status}}
        {:error, reason} -> {:error, reason}
      end
    end
  end

  defp slice_cog(cog_path, cog_lat, cog_lon, cache_dir, all_degrees, acc) do
    degrees_in_cog = Enum.filter(all_degrees, fn d -> cog_tile_key(d) == {cog_lat, cog_lon} end)
    Enum.reduce(degrees_in_cog, acc, &slice_degree(&1, cog_path, cache_dir, &2))
  end

  defp slice_degree({lat, lon}, cog_path, cache_dir, acc) do
    out_path = Path.join(cache_dir, Canopy.tile_filename(lat, lon))

    cond do
      File.exists?(out_path) -> %{acc | cached: acc.cached + 1}
      match?(:ok, slice_one(cog_path, lat, lon, out_path)) -> %{acc | sliced: acc.sliced + 1}
      true -> %{acc | failed: acc.failed + 1}
    end
  end

  defp slice_one(cog_path, lat, lon, out_path) do
    # gdal_translate -of ENVI -ot Byte -outsize 3600 3600 \
    #   -projwin <ulx> <uly> <lrx> <lry> input.tif /tmp/N32W097.canopy
    # ENVI puts the raw bytes in `out_path` plus sidecar `.hdr` files
    # which we delete. Width before height = `-outsize 3600 3600`.
    ulx = lon * 1.0
    uly = (lat + 1) * 1.0
    lrx = (lon + 1) * 1.0
    lry = lat * 1.0

    args = [
      "-q",
      "-of",
      "ENVI",
      "-ot",
      "Byte",
      "-outsize",
      Integer.to_string(@samples),
      Integer.to_string(@samples),
      "-projwin",
      Float.to_string(ulx),
      Float.to_string(uly),
      Float.to_string(lrx),
      Float.to_string(lry),
      cog_path,
      out_path
    ]

    case System.cmd("gdal_translate", args, env: %{}, stderr_to_stdout: true) do
      {_, 0} ->
        # Discard ENVI sidecar metadata.
        Enum.each([".hdr", ".aux.xml"], fn ext -> _ = File.rm(out_path <> ext) end)
        :ok

      {output, status} ->
        Logger.warning("canopy: gdal_translate failed (status #{status}): #{output}")
        _ = File.rm(out_path)
        {:error, {:gdal, status}}
    end
  end

  defp degree_tiles(bbox) do
    south = floor(bbox.south)
    north = floor(bbox.north)
    west = floor(bbox.west)
    east = floor(bbox.east)

    for lat <- south..north, lon <- west..east, do: {lat, lon}
  end

  # Lang COGs are bucketed in 3° × 3° steps anchored at multiples of 3.
  defp cog_tile_key({lat, lon}) do
    lat_anchor = lat - Integer.mod(lat, 3)
    lon_anchor = lon - Integer.mod(lon, 3)
    {lat_anchor, lon_anchor}
  end

  defp cog_filename(lat, lon) do
    lat_prefix = if lat >= 0, do: "N", else: "S"
    lon_prefix = if lon >= 0, do: "E", else: "W"
    lat_str = lat |> abs() |> Integer.to_string() |> String.pad_leading(2, "0")
    lon_str = lon |> abs() |> Integer.to_string() |> String.pad_leading(3, "0")

    "ETH_GlobalCanopyHeight_10m_2020_#{lat_prefix}#{lat_str}#{lon_prefix}#{lon_str}_Map.tif"
  end

  defp default_cache_dir, do: Application.get_env(:microwaveprop, :canopy_tiles_dir, "/data/canopy")

  defp default_staging_dir, do: Application.get_env(:microwaveprop, :canopy_staging_dir, "/data/canopy/staging")
end