diff --git a/lib/microwaveprop/buildings/ms_footprints.ex b/lib/microwaveprop/buildings/ms_footprints.ex new file mode 100644 index 00000000..1605a902 --- /dev/null +++ b/lib/microwaveprop/buildings/ms_footprints.ex @@ -0,0 +1,150 @@ +defmodule Microwaveprop.Buildings.MsFootprints do + @moduledoc """ + Bing-style quadkey math, dataset-index lookup, and tile downloader + for Microsoft's Global ML Building Footprints + (https://github.com/microsoft/GlobalMLBuildingFootprints). + + We use quadkey level 9 — each tile covers ~80 km × 65 km at DFW + latitude, which matches the rover Calculate radius. Files are + line-delimited GeoJSON inside csv.gz, hosted on Azure Blob Storage. + """ + + require Logger + + @dataset_links_url "https://minedbuildings.z5.web.core.windows.net/global-buildings/dataset-links.csv" + @region "UnitedStates" + + @doc """ + Returns the Bing/MS quadkey string for `lat, lon` at the given `zoom`. + + Encodes via Web Mercator pixel coordinates → tile XY → quadkey, per + https://learn.microsoft.com/en-us/bingmaps/articles/bing-maps-tile-system. + """ + @spec quadkey(float(), float(), pos_integer()) :: String.t() + def quadkey(lat, lon, zoom) when zoom in 1..23 do + {tx, ty} = tile_xy(lat, lon, zoom) + encode_quadkey(tx, ty, zoom) + end + + @doc """ + Returns the unique set of quadkeys at `zoom` whose tile boundary + intersects the bbox. Used to discover which MS tiles need to be + downloaded for a Calculate. + """ + @spec quadkeys_for_bbox(map(), pos_integer()) :: [String.t()] + def quadkeys_for_bbox(%{"south" => s, "north" => n, "west" => w, "east" => e}, zoom) do + {tx_min, ty_max} = tile_xy(s, w, zoom) + {tx_max, ty_min} = tile_xy(n, e, zoom) + + for tx <- tx_min..tx_max, ty <- ty_min..ty_max, uniq: true do + encode_quadkey(tx, ty, zoom) + end + end + + defp tile_xy(lat, lon, zoom) do + lat_clamped = lat |> max(-85.05112878) |> min(85.05112878) + n = :math.pow(2, zoom) + + x = (lon + 180.0) / 360.0 * n + sin_lat = :math.sin(lat_clamped * :math.pi() / 180.0) + y = (0.5 - :math.log((1 + sin_lat) / (1 - sin_lat)) / (4 * :math.pi())) * n + + {clamp_tile(trunc(x), zoom), clamp_tile(trunc(y), zoom)} + end + + defp clamp_tile(t, zoom) do + max_tile = trunc(:math.pow(2, zoom)) - 1 + t |> max(0) |> min(max_tile) + end + + defp encode_quadkey(tx, ty, zoom) do + Enum.reduce(zoom..1, "", fn i, acc -> + mask = Bitwise.bsl(1, i - 1) + digit = digit_for(tx, ty, mask) + acc <> Integer.to_string(digit) + end) + end + + defp digit_for(tx, ty, mask) do + x_bit = if Bitwise.band(tx, mask) == 0, do: 0, else: 1 + y_bit = if Bitwise.band(ty, mask) == 0, do: 0, else: 1 + Bitwise.bor(Bitwise.bsl(y_bit, 1), x_bit) + end + + @doc """ + Returns `%{quadkey => url}` for every UnitedStates row in the MS + dataset-links CSV. Result is memoized in `Microwaveprop.Cache` for + 24 hours; the index file is ~7 MB so we don't want to refetch it + on every Calculate. + """ + @spec dataset_index() :: %{String.t() => String.t()} + def dataset_index do + Microwaveprop.Cache.fetch_or_store(:ms_footprints_index, 24 * 60 * 60 * 1_000, fn -> + case Req.get(@dataset_links_url, retry: false, receive_timeout: 30_000) do + {:ok, %{status: 200, body: body}} -> parse_dataset_index(body) + other -> raise "ms footprints dataset index fetch failed: #{inspect(other)}" + end + end) + end + + defp parse_dataset_index(body) when is_binary(body) do + body + |> String.split("\n", trim: true) + |> Enum.drop(1) + |> Enum.flat_map(fn line -> + case String.split(line, ",", parts: 5) do + [@region, qk, url, _size, _date] -> [{qk, url}] + _ -> [] + end + end) + |> Map.new() + end + + @doc """ + Local cache directory for downloaded csv.gz files. Override via + `:microwaveprop, :buildings_cache_dir` to redirect into NFS in prod. + """ + @spec cache_dir() :: String.t() + def cache_dir do + Application.get_env(:microwaveprop, :buildings_cache_dir, "/data/buildings") + end + + @doc """ + Ensure `quadkey` is downloaded to disk; returns `{:ok, path}` or + `{:error, reason}`. Re-downloads are skipped if the cached file is + already present. + """ + @spec download_tile(String.t()) :: {:ok, String.t()} | {:error, term()} + def download_tile(quadkey) when is_binary(quadkey) do + path = Path.join(cache_dir(), "#{quadkey}.csv.gz") + + if File.exists?(path) do + {:ok, path} + else + do_download_tile(quadkey, path) + end + end + + defp do_download_tile(quadkey, path) do + case Map.fetch(dataset_index(), quadkey) do + :error -> + {:error, :unknown_quadkey} + + {:ok, url} -> + File.mkdir_p!(cache_dir()) + Logger.info("downloading MS footprint tile #{quadkey} from #{url}") + + case Req.get(url, retry: false, receive_timeout: 120_000) do + {:ok, %{status: 200, body: body}} -> + File.write!(path, body) + {:ok, path} + + {:ok, %{status: status}} -> + {:error, {:http, status}} + + {:error, reason} -> + {:error, reason} + end + end + end +end diff --git a/test/microwaveprop/buildings/ms_footprints_test.exs b/test/microwaveprop/buildings/ms_footprints_test.exs new file mode 100644 index 00000000..384ad19d --- /dev/null +++ b/test/microwaveprop/buildings/ms_footprints_test.exs @@ -0,0 +1,35 @@ +defmodule Microwaveprop.Buildings.MsFootprintsTest do + use ExUnit.Case, async: true + + alias Microwaveprop.Buildings.MsFootprints + + describe "quadkey/3" do + test "encodes a known DFW point at zoom 9" do + # (32.8, -97.0) — urban DFW — lies inside MS quadkey 023112330 at + # zoom 9, cross-checked against the published dataset index where + # 023112330 is the 109 MB tile covering the metro core. + assert MsFootprints.quadkey(32.8, -97.0, 9) == "023112330" + + # A polygon sampled directly from the 023112322 csv.gz file + # (lat 32.221708, lon -97.736116) confirms that tile's footprint. + assert MsFootprints.quadkey(32.221708, -97.736116, 9) == "023112322" + end + + test "Greenwich at zoom 1 hits root quadrants" do + assert MsFootprints.quadkey(0.5, 0.5, 1) == "1" + assert MsFootprints.quadkey(-0.5, -0.5, 1) == "2" + end + end + + describe "quadkeys_for_bbox/2" do + test "returns the cell containing a point and immediate neighbors at zoom 9" do + # Bbox covering ~80 km around DFW should span 4-9 zoom-9 quadkeys. + bbox = %{"south" => 32.0, "north" => 33.5, "west" => -97.8, "east" => -96.4} + keys = MsFootprints.quadkeys_for_bbox(bbox, 9) + + assert is_list(keys) + assert length(keys) > 1 + assert "023112330" in keys + end + end +end