Tutorial 2 — WorldView-3 UCSD (with bundle adjustment)

Tutorial 2 — WorldView-3 UCSD (with bundle adjustment)#

Variant of Tutorial 2 that adds bundle adjustment: the vendor RPC cameras are refined before orthorectification. Run both and compare the asp-plot reports.

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 (ASP calls orthorectification “mapprojection”).

  6. Run stereo and DEM generation.

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

from pathlib import Path

DATA = "../data/ucsd_stereo_21deg_12d"

!mkdir -p {DATA}

# Configuration shared across cells
T_PROJWIN = "476000 3637000 478000 3639000"                               # 2x2 km ROI: cliffs -> UCSD west
TR = 1.0                                                                  # orthorectification GSD; coarser than the native 0.315 m to speed processing
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.

scene_bbox runs ASP’s camera_footprint on each image and unions the results; utm_epsg picks the UTM zone at the bbox center. fetch_cop_dem.py then clips, mosaics, and reprojects the matching COP30 tiles.

from utils import scene_bbox, utm_epsg

BBOX = scene_bbox(
    (f"{DATA}/1040010007A93700_P001.NTF", f"{DATA}/1040010007A93700_P001.xml"),
    (f"{DATA}/1040010007CA4D00_P001.NTF", f"{DATA}/1040010007CA4D00_P001.xml"),
)
T_SRS = utm_epsg(BBOX)
print(f"Scene bbox: {BBOX}  UTM: {T_SRS}")

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 {BBOX} \
        --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.

--ip-per-tile 10 detects far fewer interest points than the default (about 50 per 1024² pixel tile). A two-camera solve needs far fewer, and matching time scales with the count.

%%time
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 4 \
        --ip-per-tile 10 \
        {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 1 m, coarser than the native 0.315 m GSD, to speed processing. --bundle-adjust-prefix ba/run uses our refined cameras; --t_projwin crops to the 2 × 2 km ROI.

%%time
for catid in ["1040010007A93700", "1040010007CA4D00"]:
    out = f"{DATA}/{catid}_P001_ortho_ba.tif"
    if Path(out).exists():
        print(f"{out} exists — skipping mapproject for {catid}. Delete it to reprocess.")
    else:
        !mapproject --threads 4 \
            --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 via --dem so the matcher knows the grid.

point2dem grids the point cloud at 4 m, roughly 4× the image GSD, a common rule of thumb for stereo DEM resolution.

Thread budget tuned for the 4-core Codespace floor: asp_mgm is memory-intensive and multithreads well, so ASP recommends few processes with many threads — here 1 process with 4 threads. On a larger machine, scale --processes to your core count divided by 4.

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

    !point2dem --tr 4.0 -r earth --errorimage {DATA}/stereo_ba/run-PC.tif

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

7. Full report with ICESat-2 alignment#

The asp-plot CLI produces a PDF report: input scenes, disparity, hillshades, dh vs. the reference DEM, and an ICESat-2 ATL06-SR comparison. By default it also runs pc_align to register the DEM to ICESat-2 and appends the alignment report.

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

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

The report is at: data/ucsd_stereo_21deg_12d/stereo_ba/ucsd_wv3_ba_report.pdf.

Compare against the no-BA run’s report at: data/ucsd_stereo_21deg_12d/stereo_ba/ucsd_wv3_report.pdf.