diff --git a/notebooks/WorldView/run_benchmark.sh b/notebooks/WorldView/run_benchmark.sh new file mode 100755 index 0000000..b4bdaae --- /dev/null +++ b/notebooks/WorldView/run_benchmark.sh @@ -0,0 +1,132 @@ +#!/usr/bin/env bash +# Atlanta same-pass WorldView-2 scene-combination benchmark runs (uw-cryo/asp_plot#169). +# +# usage: run_benchmark.sh /path/to/atlanta_mvs +# +# Extends the runs behind worldview_spacenet_atlanta_mvs.ipynb (5-scene +# bundle_adjust, 3- and 5-scene MVS, pairs 13-10 / 13-16 / 10-16 and their +# dem_mosaic) to the full matrix scored by worldview_spacenet_atlanta_benchmark.ipynb: +# +# 1. the seven remaining pairs -> all ten pairs, convergence 5-32 deg +# 2. dem_mosaic of all ten pairs, and of the six pairs with convergence > 15 deg +# 3. a nested quad {13,10,16,21} -> the 2->3->4->5 scene curve +# (pair 13-10 < MVS3 {13,10,16} < quad < MVS5) and a wide-spread triple {13,8,21} +# +# Every step is skipped when its product already exists, so the script can be +# re-run after an interruption. Settings are identical to run_mvs.sh / +# run_pairwise.sh (asp_mgm, subpixel-mode 9, affineepipolar, the shared 5-scene +# ba/run prefix, point2dem 1.9 m EPSG:32616 with --errorimage). +# +# Cost on an 8-thread laptop: pair ~1 h, triple ~2 h, quad ~2.5 h -> about 11-12 h, +# ~30 GB. wv_correct is ~40 s per scene. +# +# Scenes (2009-12-22, SpaceNet AOI_6_Atlanta; along-track order 10, 8, 13, 16, 21): +# nadir8 10300100023BC100_P001 +# nadir10 1030010003CAF100_P002 +# nadir13 1030010002B7D800_P002 (MVS reference; the known crop window) +# nadir16 1030010002649200_P001 +# nadir21 1030010003127500_P001 +# +# Left-image crop windows: nadir13's is the reference crop used by every existing +# run; the others were derived from the bundle_adjust clean match points that fall +# inside that crop, padded 100 px (the benchmark notebook reproduces the derivation; +# nadir10's reproduces the window run_pairwise.sh used). +set -euo pipefail + +DATA=${1:?usage: run_benchmark.sh /path/to/atlanta_mvs} +cd "$DATA" +# Append, never prepend: the ASP release bundles its own python. +export PATH="$PATH:${ASP_BIN:-$HOME/asp/dev/bin}" +THREADS=${THREADS:-8} +# The raw nadir13 / nadir10 L1B tiles live with the earlier two-scene example. +RAW_FALLBACK_DIR=${RAW_FALLBACK_DIR:-../atlanta_stereo_22deg_0d} + +N8=10300100023BC100_P001 +N10=1030010003CAF100_P002 +N13=1030010002B7D800_P002 +N16=1030010002649200_P001 +N21=1030010003127500_P001 + +crop_win () { + case "$1" in + "$N13") echo "5879 13107 12981 11894" ;; + "$N10") echo "5540 13223 13206 12931" ;; + "$N8") echo "5461 9789 13354 12980" ;; + "$N16") echo "6083 8756 12879 11461" ;; + "$N21") echo "6467 8177 12495 10687" ;; + *) echo "unknown scene $1" >&2; exit 1 ;; + esac +} + +log () { echo "=== [$(date)] $*"; } + +# 0. CCD correction. The *_corr.tif intermediates were purged after the earlier +# runs (two of them were symlinks into the two-scene example, now dangling). +for cid in $N13 $N10 $N8 $N16 $N21; do + if [ ! -s "${cid}_corr.tif" ]; then + raw="${cid}.tif" + [ -f "$raw" ] || raw="${RAW_FALLBACK_DIR}/${cid}.tif" + log "wv_correct $cid (from $raw)" + rm -f "${cid}_corr.tif" + wv_correct --threads "$THREADS" "$raw" "${cid}.xml" "${cid}_corr.tif" + fi +done + +# run_stereo [ ...] -- first scene is the reference +run_stereo () { + local out=$1; shift + local left=$1 + local imgs="" cams="" + for cid in "$@"; do + imgs="$imgs ${cid}_corr.tif" + cams="$cams ${cid}.xml" + done + log "stereo $out ($*)" + if [ ! -f "$out/run-PC.tif" ]; then + # shellcheck disable=SC2086 + parallel_stereo --stereo-algorithm asp_mgm --subpixel-mode 9 \ + --alignment-method affineepipolar \ + --left-image-crop-win $(crop_win "$left") \ + --bundle-adjust-prefix ba/run \ + $imgs $cams "$out/run" + fi + if [ ! -f "$out/run-DEM.tif" ]; then + point2dem --tr 1.9 --t_srs EPSG:32616 --errorimage "$out/run-PC.tif" + fi +} + +# mosaic [ ...] +mosaic () { + local prefix=$1; shift + if [ ! -f "${prefix}-DEM.tif" ]; then + log "dem_mosaic $prefix" + dem_mosaic "$@" -o "$prefix" + mv "${prefix}-tile-0.tif" "${prefix}-DEM.tif" + fi +} + +# 1. The seven remaining pairs (13-10, 13-16 and 10-16 exist from run_pairwise.sh). +# Left image = the reference-side scene so the known crops are reused where possible. +run_stereo stereo_pair_13_8 "$N13" "$N8" +run_stereo stereo_pair_13_21 "$N13" "$N21" +run_stereo stereo_pair_10_8 "$N10" "$N8" +run_stereo stereo_pair_10_21 "$N10" "$N21" +run_stereo stereo_pair_8_16 "$N8" "$N16" +run_stereo stereo_pair_8_21 "$N8" "$N21" +run_stereo stereo_pair_16_21 "$N16" "$N21" + +# 2. Pairwise + mosaic at five scenes: every pair, and only the well-converged ones +# (> 15 deg: 13-8 16.3, 8-16 21.5, 13-10 21.8, 8-21 26.8, 10-16 26.9, 10-21 32.3). +ALL10="stereo_pair_13_10 stereo_pair_13_16 stereo_pair_10_16 stereo_pair_13_8 stereo_pair_13_21 \ + stereo_pair_10_8 stereo_pair_10_21 stereo_pair_8_16 stereo_pair_8_21 stereo_pair_16_21" +WIDE6="stereo_pair_13_8 stereo_pair_8_16 stereo_pair_13_10 stereo_pair_8_21 stereo_pair_10_16 stereo_pair_10_21" +# shellcheck disable=SC2046 +mosaic pairwise10_mosaic $(for d in $ALL10; do echo "$d/run-DEM.tif"; done) +# shellcheck disable=SC2046 +mosaic pairwise_wide6_mosaic $(for d in $WIDE6; do echo "$d/run-DEM.tif"; done) + +# 3. Scene count: the nested quad, and a triple that spans the pass (16.3 / 10.5 / 26.8 deg). +run_stereo stereo_mvs4 "$N13" "$N10" "$N16" "$N21" +run_stereo stereo_mvs3_wide "$N13" "$N8" "$N21" + +log "ALL DONE" diff --git a/notebooks/WorldView/worldview_spacenet_atlanta_benchmark.ipynb b/notebooks/WorldView/worldview_spacenet_atlanta_benchmark.ipynb new file mode 100644 index 0000000..b271b07 --- /dev/null +++ b/notebooks/WorldView/worldview_spacenet_atlanta_benchmark.ipynb @@ -0,0 +1,1278 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "7efbc1a0", + "metadata": {}, + "source": [ + "## WorldView Scene-Combination Benchmark: SpaceNet Atlanta\n", + "\n", + "**Given N same-pass scenes, which combinations and which processing flow give the best DEM?** This notebook is the systematic version of the two-DEM comparison in [`worldview_spacenet_atlanta_mvs.ipynb`](https://asp-plot.readthedocs.io/en/latest/examples/notebooks/worldview_spacenet_atlanta_mvs.html) ([issue #169](https://github.com/uw-cryo/asp_plot/issues/169)). Everything is scored with `asp_plot.dem_benchmark.DEMBenchmark` against **one** ICESat-2 ATL06-SR sample, so every number below is directly comparable.\n", + "\n", + "Three questions, three axes of the experiment:\n", + "\n", + "1. **Convergence angle** — all ten pairs among the five scenes span 5–32°. How does single-pair accuracy depend on it, and where does the pair-selection rule of thumb (20–30°) sit on this data?\n", + "2. **Number of scenes** — a nested chain, pair 13-10 ⊂ MVS {13,10,16} ⊂ MVS {13,10,16,21} ⊂ MVS {all five}: does each added scene earn its cost?\n", + "3. **Processing flow** — joint multi-view triangulation vs. independent pairs merged with `dem_mosaic`, at three and at five scenes; and whether merging only the well-converged pairs rescues the mosaic.\n", + "\n", + "> **Status.** The runs are produced by [`run_benchmark.sh`](run_benchmark.sh) next to this notebook (~11–12 h on a laptop; see the *Processing* section). The scoring cells score **whatever DEMs exist on disk** and list the rest as pending, so the notebook is valid at every stage — as committed it covers the six DEMs from the MVS example (three pairs, their mosaic, and the 3- and 5-scene MVS runs)." + ] + }, + { + "cell_type": "markdown", + "id": "ac8f087f", + "metadata": {}, + "source": [ + "---\n", + "\n", + "## Testbed\n", + "\n", + "Five same-pass WorldView-2 scenes over Hartsfield–Jackson Atlanta International Airport (2009-12-22, SpaceNet AOI 6, public on `s3://spacenet-dataset/AOIs/AOI_6_Atlanta/`), all `wv_correct`ed and bundle-adjusted **together** (`ba/run`, five scenes), so every stereo run below shares one camera network and any DEM difference comes from the stereo strategy, not the cameras. See the MVS notebook for the download, `wv_correct`, and `bundle_adjust` steps.\n", + "\n", + "| Name | CATID / tile | Along-track order | Role |\n", + "|---|---|---|---|\n", + "| nadir10 | `1030010003CAF100_P002` | 1 | |\n", + "| nadir8 | `10300100023BC100_P001` | 2 | |\n", + "| nadir13 | `1030010002B7D800_P002` | 3 | reference scene of every MVS run (crop window `5879 13107 12981 11894`) |\n", + "| nadir16 | `1030010002649200_P001` | 4 | |\n", + "| nadir21 | `1030010003127500_P001` | 5 | |\n", + "\n", + "`bundle_adjust` reports the convergence angle of every image pair that had matches, which is the first axis of the experiment:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "964ddb7d", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-28T15:22:32.467587Z", + "iopub.status.busy": "2026-08-28T15:22:32.467328Z", + "iopub.status.idle": "2026-08-28T15:22:33.274803Z", + "shell.execute_reply": "2026-08-28T15:22:33.274286Z" + } + }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
pairconvergence (deg)p25p75num_matchesrun_dirdem_exists
013-165.1235244.8317925.4798537451stereo_pair_13_16True
116-215.3616785.3272805.3947986201stereo_pair_16_21False
28-105.4563455.0991225.8768437305stereo_pair_8_10False
313-2110.50738010.24135110.8176214813stereo_pair_13_21False
413-816.32520215.90007816.6732665550stereo_pair_13_8False
58-1621.46548921.40153521.5182223707stereo_pair_8_16False
613-1021.79153021.78685921.7955574725stereo_pair_13_10True
78-2126.83444326.72578126.9202762188stereo_pair_8_21False
810-1626.93490326.62633827.2948472822stereo_pair_10_16True
910-2132.29331932.01771432.6229301860stereo_pair_10_21False
\n", + "
" + ], + "text/plain": [ + " pair convergence (deg) p25 p75 num_matches \\\n", + "0 13-16 5.123524 4.831792 5.479853 7451 \n", + "1 16-21 5.361678 5.327280 5.394798 6201 \n", + "2 8-10 5.456345 5.099122 5.876843 7305 \n", + "3 13-21 10.507380 10.241351 10.817621 4813 \n", + "4 13-8 16.325202 15.900078 16.673266 5550 \n", + "5 8-16 21.465489 21.401535 21.518222 3707 \n", + "6 13-10 21.791530 21.786859 21.795557 4725 \n", + "7 8-21 26.834443 26.725781 26.920276 2188 \n", + "8 10-16 26.934903 26.626338 27.294847 2822 \n", + "9 10-21 32.293319 32.017714 32.622930 1860 \n", + "\n", + " run_dir dem_exists \n", + "0 stereo_pair_13_16 True \n", + "1 stereo_pair_16_21 False \n", + "2 stereo_pair_8_10 False \n", + "3 stereo_pair_13_21 False \n", + "4 stereo_pair_13_8 False \n", + "5 stereo_pair_8_16 False \n", + "6 stereo_pair_13_10 True \n", + "7 stereo_pair_8_21 False \n", + "8 stereo_pair_10_16 True \n", + "9 stereo_pair_10_21 False " + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import os\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "\n", + "# Set the base directory for your processing\n", + "directory = \"~/Desktop/asp-plot-examples/atlanta_mvs/\"\n", + "directory = os.path.expanduser(directory)\n", + "\n", + "SCENES = {\n", + " \"1030010003CAF100_P002\": \"10\",\n", + " \"10300100023BC100_P001\": \"8\",\n", + " \"1030010002B7D800_P002\": \"13\",\n", + " \"1030010002649200_P001\": \"16\",\n", + " \"1030010003127500_P001\": \"21\",\n", + "}\n", + "\n", + "# ba/run-convergence_angles.txt: left right 25% 50% 75% num_matches (degrees)\n", + "conv = pd.read_csv(\n", + " f\"{directory}ba/run-convergence_angles.txt\",\n", + " sep=r\"\\s+\",\n", + " comment=\"#\",\n", + " names=[\"left\", \"right\", \"p25\", \"p50\", \"p75\", \"num_matches\"],\n", + ")\n", + "for col in (\"left\", \"right\"):\n", + " conv[col] = conv[col].str.replace(\"_corr.tif\", \"\", regex=False).map(SCENES)\n", + "\n", + "\n", + "def pair_dir(a, b):\n", + " \"\"\"stereo_pair__ as produced by run_pairwise.sh / run_benchmark.sh (either order).\"\"\"\n", + " for left, right in ((a, b), (b, a)):\n", + " d = f\"stereo_pair_{left}_{right}\"\n", + " if os.path.isdir(f\"{directory}{d}\"):\n", + " return d\n", + " return f\"stereo_pair_{a}_{b}\"\n", + "\n", + "\n", + "conv[\"pair\"] = conv.apply(lambda r: f\"{r.left}-{r.right}\", axis=1)\n", + "conv[\"run_dir\"] = conv.apply(lambda r: pair_dir(r.left, r.right), axis=1)\n", + "conv[\"dem_exists\"] = conv[\"run_dir\"].map(lambda d: os.path.exists(f\"{directory}{d}/run-DEM.tif\"))\n", + "conv = conv.sort_values(\"p50\").reset_index(drop=True)\n", + "conv[[\"pair\", \"p50\", \"p25\", \"p75\", \"num_matches\", \"run_dir\", \"dem_exists\"]].rename(\n", + " columns={\"p50\": \"convergence (deg)\"}\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "dafc7d0a", + "metadata": {}, + "source": [ + "Ten pairs, evenly spread from 5° to 32°, with three of them already run for the MVS example (5.1°, 21.8°, 26.9°). The match count falls with convergence, as expected for the same interest-point detector seeing an increasingly different view.\n", + "\n", + "### Crop windows for the new pairs\n", + "\n", + "Every existing run crops the reference scene (nadir13) to the airport window and lets ASP find the corresponding region in the other image. A pair whose *left* image is not nadir13 needs its own `--left-image-crop-win`, which is derived the same way `run_pairwise.sh` derived nadir10's: take the `bundle_adjust` clean match points between nadir13 and that scene, keep those inside the nadir13 crop, and pad the bounding box of their partners by 100 px. nadir10's window reproduces the one used for `stereo_pair_10_16`, which validates the method; the other three are what `run_benchmark.sh` uses." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "1811930d", + "metadata": { + "execution": { + "iopub.execute_input": "2026-08-28T15:22:33.276153Z", + "iopub.status.busy": "2026-08-28T15:22:33.276044Z", + "iopub.status.idle": "2026-08-28T15:22:34.362533Z", + "shell.execute_reply": "2026-08-28T15:22:34.362109Z" + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "WARNING:asp_plot.stereo:\n", + "\n", + "No reference DEM found in log files. Please supply the reference DEM you used during stereo processing (or another reference DEM) if you would like to see some difference maps.\n", + "\n", + "\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "ASP DEM: /Users/ben/Desktop/asp-plot-examples/atlanta_mvs/stereo_mvs3/run-DEM.tif\n", + "\n", + "nadir10 775 of 4725 clean matches in the reference crop -> --left-image-crop-win 5540 13223 13206 12931\n", + "nadir8 890 of 5550 clean matches in the reference crop -> --left-image-crop-win 5461 9789 13354 12980\n", + "nadir16 1176 of 7451 clean matches in the reference crop -> --left-image-crop-win 6083 8756 12879 11461\n", + "nadir21 785 of 4813 clean matches in the reference crop -> --left-image-crop-win 6467 8177 12495 10687\n" + ] + } + ], + "source": [ + "from asp_plot.stereo import StereoPlotter\n", + "\n", + "XOFF, YOFF, W, H = 5879, 13107, 12981, 11894 # nadir13 reference crop\n", + "REF = \"1030010002B7D800_P002\"\n", + "\n", + "stereo_plotter = StereoPlotter(directory, \"stereo_mvs3\")\n", + "crops = {\"13\": (XOFF, YOFF, W, H)}\n", + "for cid, name in SCENES.items():\n", + " if cid == REF:\n", + " continue\n", + " match_fn = f\"{directory}ba/run-{REF}_corr__{cid}_corr-clean.match\"\n", + " match_df = stereo_plotter.get_match_point_df(match_point_fn=match_fn)\n", + " inwin = match_df[match_df.x1.between(XOFF, XOFF + W) & match_df.y1.between(YOFF, YOFF + H)]\n", + " pad = 100\n", + " x0, y0 = int(inwin.x2.min()) - pad, int(inwin.y2.min()) - pad\n", + " x1, y1 = int(inwin.x2.max()) + pad, int(inwin.y2.max()) + pad\n", + " crops[name] = (x0, y0, x1 - x0, y1 - y0)\n", + " print(\n", + " f\"nadir{name:<3} {len(inwin):5d} of {len(match_df):5d} clean matches in the reference crop\"\n", + " f\" -> --left-image-crop-win {x0} {y0} {x1 - x0} {y1 - y0}\"\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "cab5b78c", + "metadata": {}, + "source": [ + "---\n", + "\n", + "## Processing\n", + "\n", + "[`run_benchmark.sh`](run_benchmark.sh) runs everything the matrix needs that the MVS example did not, with the settings of the existing runs (`asp_mgm`, `--subpixel-mode 9`, `affineepipolar`, the shared `ba/run` prefix, `point2dem --tr 1.9 --t_srs EPSG:32616 --errorimage`). Every step is skipped when its product exists, so it can be re-run after an interruption:\n", + "\n", + "```bash\n", + "$ bash notebooks/WorldView/run_benchmark.sh ~/Desktop/asp-plot-examples/atlanta_mvs\n", + "```\n", + "\n", + "| Step | Products | Answers |\n", + "|---|---|---|\n", + "| `wv_correct` × 5 | `_corr.tif` | (the intermediates were purged after the earlier runs) |\n", + "| the seven remaining pairs | `stereo_pair_13_8`, `_13_21`, `_10_8`, `_10_21`, `_8_16`, `_8_21`, `_16_21` | axis 1: NMAD vs. convergence over all ten pairs |\n", + "| `dem_mosaic` of all ten pairs / of the six pairs > 15° | `pairwise10_mosaic-DEM.tif`, `pairwise_wide6_mosaic-DEM.tif` | axis 3: the five-scene pairwise flows |\n", + "| nested quad `{13,10,16,21}` | `stereo_mvs4/` | axis 2: the 2→3→4→5 chain |\n", + "| wide-spread triple `{13,8,21}` | `stereo_mvs3_wide/` | axis 2/3: does the *choice* of three scenes matter as much as their number? |\n", + "\n", + "About 11–12 hours on an 8-thread laptop (pair ≈ 1 h, triple ≈ 2 h, quad ≈ 2.5 h) and ~30 GB." + ] + }, + { + "cell_type": "markdown", + "id": "cdb9965e", + "metadata": {}, + "source": [ + "---\n", + "\n", + "## Scoring every DEM against one ICESat-2 sample\n", + "\n", + "`DEMBenchmark` scores each candidate with the report's recipe — the cached ATL06-SR parquet replayed (no SlideRule request), ESA WorldCover water returns dropped, a 3σ outlier cut per DEM — and a `pc_align --compute-translation-only` of its own, whose products go under `dem_benchmark/