Tutorial 2 — WorldView-3 UCSD

Tutorial 2 — WorldView-3 UCSD#

The full commercial-stereo recipe. Two WV3 panchromatic NTFs to a bundle-adjusted, orthorectified, ICESat-2-aligned ~1 m DEM.

We chose this dataset because:

  • Openly available. SpaceNet CORE3D imagery is on a public AWS bucket with --no-sign-request. No Vantor license, no API keys.

  • Well-conditioned geometry: 21° convergence, 0.37 base-to-height, 12-day temporal separation.

  • Mixed landscape: Black’s Beach sea cliffs, La Jolla Shores residential canyons, UCSD’s western campus.

  • ICESat-2 has good coverage at this latitude, so the alignment step shows real improvement.

Pair details#

Catid

Date

Off-nadir

Position

1040010007A93700

2015-02-12

8.4°

left

1040010007CA4D00

2015-02-24

12.9°

right

ROI: 2 × 2 km from Black’s Beach sea cliffs east through La Jolla Shores to UCSD’s western campus. UTM zone 11N (EPSG:32611).

What we’ll do#

  1. Download both NTFs and metadata from SpaceNet S3.

  2. Visualize the stereo geometry (skyplot, footprints).

  3. Clip a Copernicus DEM tile from AWS Open Data.

  4. Bundle adjust to refine the camera models.

  5. Run mapproject to orthorectify both images onto the COP-DEM grid, cropped to the ROI.

  6. Run stereo and DEM generation.

  7. Generate a full report including ICESat-2 alignment via pc_align.

Note

ASP calls this step “mapprojection” in its toolchain (the binary is mapproject, the bundle-adjust flag is --mapproj-dem). We use “orthorectification” — the more standard photogrammetry term — for the concept while keeping ASP’s tool names verbatim.

from pathlib import Path

DATA = "../data/ucsd_stereo_21deg_12d"

!mkdir -p {DATA}

# Configuration shared across cells
T_SRS = "EPSG:32611"                                                      # UTM 11N for San Diego
T_PROJWIN = "476000 3637000 478000 3639000"                               # 2x2 km ROI: cliffs -> UCSD west
TR = 0.315                                                                # output GSD (more nadir image's meanProductGSD)
REF_DEM = f"{DATA}/ref/cop30_ucsd.tif"

1. Download the imagery from SpaceNet#

Each WV3 scene comes as an NTF (the imagery) plus a TAR (metadata, including the IMD with RPC coefficients and an XML camera). The download script handles the AWS fetch + tar extract and renames files to short, human-readable names (<CATID>_P001.NTF, <CATID>_P001.xml).

!bash ../scripts/download_worldview_ucsd.sh {DATA}
!ls -lh {DATA}/*.NTF {DATA}/*.xml

2. Stereo geometry#

StereoGeometryPlotter reads the XML metadata and builds a skyplot (where each satellite was looking, in azimuth/elevation) plus a map view of the two image footprints.

from asp_plot.stereo_geometry import StereoGeometryPlotter

sgp = StereoGeometryPlotter(DATA)
print(f"UTM EPSG: {sgp.get_pair_utm_epsg()}")
print(f"Scene bounds (lon/lat): {sgp.get_scene_bounds()}")
sgp.dg_geom_plot()

3. Reference DEM from AWS Copernicus DEM#

Orthorectification needs a coarse reference DEM. Copernicus GLO-30 is on AWS Open Data; no API key, no auth.

fetch_cop_dem.py clips a bounding box, mosaics any tiles that overlap, and reprojects to the target SRS at the requested resolution.

Path(f"{DATA}/ref").mkdir(exist_ok=True)

if Path(REF_DEM).exists():
    print(f"{REF_DEM} exists — skipping COP30 fetch. Delete it to redownload.")
else:
    !python ../scripts/fetch_cop_dem.py \
        --bbox -117.27 32.84 -117.20 32.92 \
        --t-srs {T_SRS} \
        --tr 30 \
        --out {REF_DEM}

!ls -lh {REF_DEM}

4. Bundle adjustment#

Bundle-adjust on the full images (not the cropped ROI) so tie points are spread across the whole scene. The optimizer refines each camera’s position and pointing to minimize reprojection error.

Flags:

  • --ip-per-image 10000: dense interest points across the scene.

  • --tri-weight 0.1 --tri-robust-threshold 0.1: soft anchor on ground-point movement; keeps the solution from drifting.

  • --camera-weight 0: lets the cameras move freely (default would penalize that).

if Path(f"{DATA}/ba/run-final_residuals_pointmap.csv").exists():
    print(f"{DATA}/ba/run-final_residuals_pointmap.csv exists — skipping bundle_adjust. Delete {DATA}/ba/ to reprocess.")
else:
    !bundle_adjust \
        --threads 8 \
        --ip-per-image 10000 \
        --tri-weight 0.1 \
        --tri-robust-threshold 0.1 \
        --camera-weight 0 \
        {DATA}/1040010007A93700_P001.NTF \
        {DATA}/1040010007CA4D00_P001.NTF \
        {DATA}/1040010007A93700_P001.xml \
        {DATA}/1040010007CA4D00_P001.xml \
        -o {DATA}/ba/run

!ls {DATA}/ba/ | head

Plot the bundle-adjustment residuals before vs. after optimization.

from asp_plot.bundle_adjust import ReadBundleAdjustFiles, PlotBundleAdjustFiles

ba_reader = ReadBundleAdjustFiles(DATA, "ba")
initial_gdf, final_gdf = ba_reader.get_initial_final_residuals_gdfs()

plotter = PlotBundleAdjustFiles(
    [initial_gdf, final_gdf],
    lognorm=True,
    title="UCSD bundle adjust residuals",
)
plotter.plot_n_gdfs(column_name="mean_residual", cbar_label="Mean residual (px)")

5. Mapproject#

Resample both NTFs onto the COP-DEM grid at full WV3 GSD (0.315 m). --bundle-adjust-prefix ba/run uses our refined cameras; --t_projwin crops to the 2 × 2 km ROI.

for catid in ["1040010007A93700", "1040010007CA4D00"]:
    out = f"{DATA}/{catid}_P001_corr_map.tif"
    if Path(out).exists():
        print(f"{out} exists — skipping mapproject for {catid}. Delete it to reprocess.")
    else:
        !mapproject --threads 8 \
            --t_srs {T_SRS} --tr {TR} \
            --t_projwin {T_PROJWIN} \
            --bundle-adjust-prefix {DATA}/ba/run \
            {REF_DEM} \
            {DATA}/{catid}_P001.NTF \
            {DATA}/{catid}_P001.xml \
            {out}

6. Stereo + DEM#

Stereo on the orthorectified pair. The inputs are now nearly aligned in geographic space, so the matcher’s search range is small. Pass the reference DEM as the final argument so the matcher knows the grid.

Thread budget tuned for 8 cores: --processes 4 --threads-multiprocess 2 --threads-singleprocess 4 (4 × 2 = 8 in the parallel stage).

if Path(f"{DATA}/stereo/run-DEM.tif").exists():
    print(f"{DATA}/stereo/run-DEM.tif exists — skipping stereo + point2dem. Delete {DATA}/stereo/ to reprocess.")
else:
    !parallel_stereo \
        --threads-multiprocess 2 \
        --threads-singleprocess 4 \
        --processes 4 \
        --stereo-algorithm asp_mgm \
        --subpixel-mode 9 \
        --bundle-adjust-prefix {DATA}/ba/run \
        {DATA}/1040010007A93700_P001_corr_map.tif \
        {DATA}/1040010007CA4D00_P001_corr_map.tif \
        {DATA}/1040010007A93700_P001.xml \
        {DATA}/1040010007CA4D00_P001.xml \
        {DATA}/stereo/run \
        {REF_DEM}

    !point2dem --tr 1.0 -r earth --errorimage {DATA}/stereo/run-PC.tif

!ls -lh {DATA}/stereo/run-DEM.tif

7. Inspect with asp_plot#

Disparity map and high-resolution hillshade.

from asp_plot.stereo import StereoPlotter

plotter = StereoPlotter(
    DATA,
    "stereo/",
    reference_dem=REF_DEM,
    title="UCSD WV3 — stereo outputs",
)
plotter.plot_disparity(unit="meters", quiver=True)

plotter.plot_detailed_hillshade(subset_km=1)

ICESat-2 ATL06-SR vs the stereo DEM, before alignment. The asp_plot CLI step below runs pc_align and re-does this against the aligned DEM.

from utils import icesat2_check
icesat2_check(f"{DATA}/stereo/run-DEM.tif", directory=DATA)

8. Full report with ICESat-2 alignment#

The asp_plot CLI does what we just did manually, plus an ICESat-2 ATL06-SR comparison. With --pc_align True, it also runs an alignment step that registers the DEM to ICESat-2.

Expect this cell to take a few minutes; it fetches ATL06-SR data via SlideRule and runs pc_align.

if Path(f"{DATA}/stereo/ucsd_wv3_report.pdf").exists():
    print(f"{DATA}/stereo/ucsd_wv3_report.pdf exists — skipping asp_plot CLI. Delete it to regenerate.")
else:
    !asp_plot \
        --directory {DATA} \
        --bundle_adjust_directory ba \
        --stereo_directory stereo \
        --reference_dem {DATA}/ref/cop30_ucsd.tif \
        --subset_km 0.5 \
        --report_filename ucsd_wv3_report.pdf

Open data/ucsd_stereo_21deg_12d/stereo/ucsd_wv3_report.pdf from the file tree.

What’s next#

→ For more advanced WV processing (jitter correction with jitter_solve, no-orthorectification variants, multi-pair scene selection), see asp_plot’s WorldView notebooks.

→ For planetary stereo (Mars CTX/HiRISE/MOC, Lunar LRO NAC), see asp_plot’s planetary notebooks. The recipe is the same; the sensor-prep step changes.