Pléiades Neo Tri-Stereo Processing: Marseille Example

This notebook demonstrates a tri-stereo processing workflow for Airbus Pléiades Neo 3 imagery using the NASA Ames Stereo Pipeline (ASP), on the free Airbus sample dataset “Pléiades Neo Tristereo – Marseille” (acquired 2021-11-07, ~10:29 UTC).

It also demonstrates asp_plot’s support for Airbus DIMAP metadata: the stereo geometry, satellite position/orientation, and report tooling all work directly from the DIM_*.XML product files delivered with the imagery.

Tri-stereo scenes (PAN, 0.30 m GSD, each ~1 Gpix):

Scene (DATASET_NAME)	Time (UTC)	Incidence	Along-track view	Role
PNEO3_202111071029126_PAN_SEN	10:29:12	22.2°	+12.7°	fore
PNEO3_202111071029296_PAN_SEN	10:29:29	17.9°	+2.4°	nadir
PNEO3_202111071029456_PAN_SEN	10:29:45	19.5°	−7.7°	aft

Pairwise convergence angles (measured identically by the DIMAP metadata and by ASP bundle_adjust interest points):

• fore–nadir: 11.1°

• nadir–aft: 10.6°

• fore–aft: 21.7° (B:H 0.38) — the strongest pair geometry

The scene footprint covers central Marseille — the Vieux-Port, Notre-Dame de la Garde, and the coastline — a good mix of dense urban fabric, terrain relief, and water.

---

Data Retrieval and Preparation

The sample is delivered as a ~3.8 GB zip (Sample_Stereo_Marseille.zip) after filling out the Airbus sample-data form. Unzipped, the delivery keeps each scene in its own folder:

000251471_1_1_STE_A/
  IMG_01_PNEO3_PAN/   # nadir  (DIM_*.XML, RPC_*.XML, IMG_*_R1C1.JP2, ...)
  IMG_02_PNEO3_PAN/   # fore
  IMG_03_PNEO3_PAN/   # aft
  IMG_0*_PNEO3_MS-FS/ # multispectral (not used here)

JP2 → GeoTIFF conversion

ASP works best with tiled GeoTIFFs. Convert each PAN JP2 and copy its DIMAP camera file alongside with a matching basename (GDAL needs the JP2OpenJPEG driver — on conda-forge installs of GDAL ≥ 3.9 that is the separate libgdal-jp2openjpeg package):

$ cd pleiades_marseille
$ for sub in 02:202111071029126 01:202111071029296 03:202111071029456; do
    img=${sub%%:*}; stamp=${sub##*:}
    gdal_translate -of GTiff \
      -co TILED=YES -co BLOCKXSIZE=256 -co BLOCKYSIZE=256 \
      -co COMPRESS=LZW -co PREDICTOR=2 -co BIGTIFF=IF_SAFER \
      000251471_1_1_STE_A/IMG_${img}_PNEO3_PAN/IMG_PNEO3_*_P_R1C1.JP2 \
      ${stamp}_PAN.tif
    cp 000251471_1_1_STE_A/IMG_${img}_PNEO3_PAN/DIM_PNEO3_*.XML ${stamp}_PAN.XML
  done

The DIM_*.XML files are the exact linescan camera models consumed by ASP’s -t pleiades session. The ASP documentation specifically recommends the exact model over the RPC sidecars for Pléiades Neo, since the RPC models of a tri-stereo collect can be mutually inconsistent and inflate triangulation error. (A cam_test comparison of the two models for these scenes shows ~2.5 px of disagreement.)

Stereo Geometry Analysis

Before processing, analyze the acquisition geometry directly from the DIMAP metadata. With three scenes, asp_plot produces an N-scene overview plus one figure per pair:

# Set the base directory for your processing
directory = "~/Desktop/asp-plot-examples/pleiades_marseille/"

from asp_plot.stereo_geometry import StereoGeometryPlotter

sgp = StereoGeometryPlotter(directory)
scene_bounds = sgp.parser.get_catid_dicts()[0]["geom"].bounds
map_crs = "EPSG:32631"  # UTM 31N for Marseille

print(f"Scene footprint (lon/lat): {scene_bounds}")

_ = sgp.stereo_geom_plot()

Scene footprint (lon/lat): (5.307138681631315, 43.238916188302895, 5.431171743354403, 43.34072049054464)

The fore–aft pair (21.7° convergence, B:H 0.38) is the strongest geometry for DEM generation; the fore–nadir and nadir–aft pairs (~11°) trade height precision for reduced occlusion in the dense urban core.

Satellite Position and Orientation Data

The DIMAP metadata also carries the refined ephemeris (ECEF position/velocity) and attitude quaternions used by the exact camera model. asp_plot plots each scene’s orbital track and its roll/pitch/yaw relative to the local orbital (LVLH) frame.

Two DIMAP-specific notes:

• Airbus quaternions are scalar-first (Q0 is the scalar part), unlike the scalar-last layout of WorldView metadata.

• DIMAP reports no ephemeris or attitude covariance (unlike WorldView XML), so those panels are annotated as not provided.

The roll/pitch/yaw values directly reflect the tri-stereo pointing: roll ≈ −15.5° (across-track) for all three scenes, with pitch stepping from +12.5° (fore) through +2.4° (nadir) to −7.5° (aft).

sgp.satellite_position_orientation_plot()

---

Bundle Adjustment

Bundle adjust all three scenes together to produce one self-consistent set of cameras for the multi-view run below. ASP’s modern defaults (--camera-weight 0, --tri-weight 0.1) match the documented recommendation for Pléiades:

$ bundle_adjust -t pleiades \
    202111071029126_PAN.tif 202111071029296_PAN.tif 202111071029456_PAN.tif \
    202111071029126_PAN.XML 202111071029296_PAN.XML 202111071029456_PAN.XML \
    --datum WGS84 -o ba/ba

The refined DIMAP geometry is already excellent: mean reprojection errors land at 0.19–0.22 px and the optimized camera centers move only millimeters. Bundle adjustment is still worthwhile insurance for the joint multi-view triangulation below.

import contextily as ctx
from asp_plot.bundle_adjust import ReadBundleAdjustFiles, PlotBundleAdjustFiles

bundle_adjust_directory = "ba/"
stereo_directory = "stereo_mvs/"

ctx_kwargs = {
    "crs": map_crs,
    "source": ctx.providers.Esri.WorldImagery,
    "attribution_size": 0,
    "alpha": 0.5,
}

ba_files = ReadBundleAdjustFiles(directory, bundle_adjust_directory)
resid_initial_gdf, resid_final_gdf = ba_files.get_initial_final_residuals_gdfs(residuals_in_meters=True)

plotter = PlotBundleAdjustFiles(
    [resid_initial_gdf, resid_final_gdf],
    lognorm=True,
    title="Bundle Adjust Initial and Final Residuals",
)
plotter.plot_n_gdfs(
    column_name="mean_residual",
    cbar_label="Mean residual (px)",
    map_crs=map_crs,
    **ctx_kwargs,
)

WARNING:asp_plot.bundle_adjust:

Could not read stereopair metadata to get mean GSD. Residuals will remain in pixels.

WARNING:asp_plot.bundle_adjust:

Could not read stereopair metadata to get mean GSD. Residuals will remain in pixels.

Camera Position and Orientation Changes

The refined DIMAP geometry needed almost nothing. Read from the ba/ folder alone (*.adjust, camera_offsets.txt, triangulation_offsets.txt; no original cameras needed), the three camera centers moved by 1–9 mm and rotated by less than 4e-5°, and the triangulated points moved by a median of 7–11 mm. The figure keeps its axes in meters, so these bars are tiny by design — this is what a well-calibrated sensor looks like next to the 24–70 m shifts of the WorldView Atlanta example. The same figure is available from the bundle_adjust_cameras command-line tool and is a page of the asp_report PDF.

from asp_plot.bundle_adjust import ReadBundleAdjustCameras, PlotBundleAdjustCameras

cam_gdf = ReadBundleAdjustCameras(directory, bundle_adjust_directory).get_camera_optimization_gdf(
    map_crs=int(map_crs.split(":")[-1])
)
fig = PlotBundleAdjustCameras(cam_gdf, title="Camera changes from bundle_adjust (Marseille, Pléiades Neo tri-stereo)").summary_plot()

---

Multi-View Stereo: All Three Scenes

ASP supports true multi-view stereo: pass N images and N cameras, and the first image is the reference — disparities are computed from it to every other image, followed by a joint triangulation into a single point cloud. Here the nadir scene is the natural reference.

The full scenes are ~1 Gpix each, so we restrict correlation to a ~11,500 × 13,000 px window (~3.5 × 4 km) over central Marseille using --left-image-crop-win on the reference image — this keeps the exact camera models intact (no image cropping) while bounding the correlation work. (The crop-window pixel coordinates for a target lon/lat box can be derived from the Dataset_Extent vertices of the scene’s DIMAP file, which map the footprint corners to image row/column.)

$ parallel_stereo -t pleiades \
    --stereo-algorithm asp_mgm --subpixel-mode 9 \
    --left-image-crop-win 11552 10513 11523 12954 \
    --bundle-adjust-prefix ba/ba \
    202111071029296_PAN.tif 202111071029126_PAN.tif 202111071029456_PAN.tif \
    202111071029296_PAN.XML 202111071029126_PAN.XML 202111071029456_PAN.XML \
    stereo_mvs/run

# ~1.2 m DEM (4x the 0.3 m input GSD) + triangulation error image
$ point2dem --tr 1.2 --t_srs EPSG:32631 --errorimage stereo_mvs/run-PC.tif

Note that the ASP documentation discourages this mode in favor of pairwise runs merged with dem_mosaic; it is demonstrated here because a same-pass tri-stereo is the best case for it (identical illumination, rigid relative geometry).

---

Stereo Results Visualization

ASP keeps the per-pair intermediate products of a multi-view run (aligned scenes, match points, disparity maps) in its run-pair*/ subdirectories, with only the joint products (point cloud, DEM, triangulation error) at the top level. asp_plot detects this layout and renders the scenes, match points, and disparity one figure per pair — the reference (left) image is the same for every pair:

from asp_plot.scenes import ScenePlotter

scene_plotter = ScenePlotter(directory, stereo_directory, title="Input Scenes")
scene_plotter.plot_scenes()

[]

from asp_plot.stereo import StereoPlotter

stereo_plotter = StereoPlotter(directory, stereo_directory)
stereo_plotter.title = "Stereo Match Points"
stereo_plotter.plot_match_points()

WARNING:asp_plot.stereo:

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.

ASP DEM: /Users/ben/Desktop/asp-plot-examples/pleiades_marseille/stereo_mvs/run-DEM.tif

[]

stereo_plotter.title = "Disparity (pixels)"
stereo_plotter.plot_disparity(unit="pixels", quiver=True)

[]

The joint products of the multi-view run — the DEM with triangulation error, and detailed hillshades:

stereo_plotter.title = "Multi-View Stereo DEM Results (nadir reference, 3 scenes)"
stereo_plotter.plot_dem_results()

Plotting DEM results. This can take a minute for large inputs.

WARNING:asp_plot.utils:Could not find ('*diff.tif',) in /Users/ben/Desktop/asp-plot-examples/pleiades_marseille/stereo_mvs/. Some plots may be missing.

stereo_plotter.title = "Hillshade with Details"
stereo_plotter.plot_detailed_hillshade(subset_km=0.5)

---

Full Report Generation

Generate the standard asp_plot PDF report for the multi-view run (the ICESat-2 altimetry comparison runs automatically for Earth data; scene panels carry the “© Airbus DS” attribution). The Input Scenes, Match Points, and Disparity sections contain one figure per run-pair*/ pair, and the joint DEM, hillshade, and altimetry sections follow:

For a committed example of a multi-view report (same per-pair sections, WorldView data), see the WorldView_Atlanta_MVS report.

!asp_report \
    --directory ~/Desktop/asp-plot-examples/pleiades_marseille/ \
    --bundle-adjust-prefix ba/ \
    --stereo-directory stereo_mvs/ \
    --report-filename ~/Desktop/asp-plot-examples/pleiades_marseille/Pleiades_Neo_Marseille_MVS-asp-report.pdf