Comprehensive WorldView Example with asp_plot

This notebook demonstrates the full capabilities of asp_plot for visualizing NASA Ames Stereo Pipeline (ASP) processing results from WorldView satellite imagery.

The example uses WorldView-3 data from Utqiagvik, Alaska (April 17, 2022) processed with ASP.

---

Processing Overview

This notebook covers:

1. Full report generation - Automated comprehensive PDF report from CLI call
2. Processing parameters - Extract and display ASP command history
3. Scene visualization - Display input satellite imagery
4. Stereo geometry - Analyze acquisition geometry and viewing angles
5. Bundle adjustment - Visualize camera optimization residuals
6. Stereo results - Display DEMs, hillshades, disparity maps, and match points
7. ICESat-2 validation - Compare DEMs with altimetry data
8. CSM camera analysis - Compare original and optimized camera models (Using example data from Salar de Uyuni)

Each section can be run independently for example modular analysis using the package.

View complete notebook with all outputs and figures

View the report generated by the asp_plot CLI call

1. Full Report Generation

The asp_plot CLI tool generates a comprehensive PDF report combining all visualizations below. This is the quickest way to get an overview of your stereo processing results.

directory = "/Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/"
ba_directory = "ba/"
stereo_directory = "stereo/"

!asp_plot \
  --directory $directory \
  --bundle_adjust_directory $ba_directory \
  --stereo_directory $stereo_directory \
  --map_crs EPSG:32604 \
  --dem_gsd 1 \
  --subset_km 1 \
  --add_basemap True \
  --plot_icesat True \
  --plot_geometry True

Processing ASP files in /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/


Reference DEM: /Users/ben/Dropbox/UW_Shean/COP/COP30_utqiagvik_lzw-adj_proj.tif


ASP DEM: /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/20220417_2252_1040010074793300_1040010075633C00-DEM_1m.tif


Reference DEM: /Users/ben/Dropbox/UW_Shean/COP/COP30_utqiagvik_lzw-adj_proj.tif


ASP DEM: /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/20220417_2252_1040010074793300_1040010075633C00-DEM_1m.tif

Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/00.png
Plotting DEM results. This can take a minute for large inputs.
WARNING:asp_plot.stereo:

Found a DEM of difference: /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/20220417_2252_1040010074793300_1040010075633C00-DEM_1m_COP30_utqiagvik_lzw-adj_proj_diff.tif.
Using that for difference map plotting.


Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/01.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/02.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/03.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/04.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/05.png

ICESat-2 ATL06 request processing for: all
{'poly': [{'lon': -156.83342896997698, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.26658548886508}], 'res': 20, 'len': 40, 'ats': 20, 'maxi': 6, 'samples': {'esa_worldcover': {'asset': 'esa-worldcover-10meter'}}, 'cnf': 'atl03_high', 'srt': -1, 'cnt': 10}
Existing file found, reading in: atl06sr_all.parquet
Filtering ATL06-SR all

ICESat-2 ATL06 request processing for: ground
{'poly': [{'lon': -156.83342896997698, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.26658548886508}], 'res': 20, 'len': 40, 'ats': 20, 'maxi': 6, 'samples': {'esa_worldcover': {'asset': 'esa-worldcover-10meter'}}, 'cnf': 'atl03_low', 'srt': -1, 'cnt': 5, 'atl08_class': 'atl08_ground'}
Existing file found, reading in: atl06sr_ground.parquet
Filtering ATL06-SR ground

Filtering ATL06 with 15 day pad, 45 day, 91 day pad, and seasonal pad around 2022-04-17 22:52:18.495475 for: all

Filtering ATL06 with 15 day pad, 45 day, 91 day pad, and seasonal pad around 2022-04-17 22:52:18.495475 for: ground

icesat_minus_dem not found in ATL06 dataframe: all. Running differencing first.

Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/06.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/07.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/08.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/09.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/10.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/11.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/12.png
Figure saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/tmp_asp_report_plots/13.png


Report saved to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/asp_plot_report_WV03_20220417_1040010074793300_1040010075633C00_20251024_190509.pdf

---

Individual Plots

The sections below demonstrate modular usage of asp_plot for detailed analysis and customization.

2. Processing Parameters

Extract and display the command-line parameters used for bundle adjustment, stereo processing, and DEM generation. This provides full traceability of processing settings.

%load_ext autoreload
%autoreload 2

from asp_plot.processing_parameters import ProcessingParameters

processing_parameters = ProcessingParameters(
    processing_directory=directory,
    bundle_adjust_directory=ba_directory,
    stereo_directory=stereo_directory
)
processing_parameters_dict = processing_parameters.from_log_files()

print(f"Processed on: {processing_parameters_dict['processing_timestamp']}\n")

print(f"Reference DEM: {processing_parameters_dict['reference_dem']}\n")

print(f"Bundle adjustment ({processing_parameters_dict['bundle_adjust_run_time']}):\n")
print(processing_parameters_dict["bundle_adjust"])

print(f"\nStereo ({processing_parameters_dict['stereo_run_time']}):\n")
print(processing_parameters_dict["stereo"])

print(f"\nPoint2dem ({processing_parameters_dict['point2dem_run_time']}):\n")
print(processing_parameters_dict["point2dem"])

Processed on: 2024-04-14 17:55:43

Reference DEM: /Users/ben/Dropbox/UW_Shean/COP/COP30_utqiagvik_lzw-adj_proj.tif

Bundle adjustment (0 hours and 7 minutes):

bundle_adjust -t dg --weight-image /nobackup/bpurint1/data/utqiagvik/WV/utqiagvik_wv_EE/2022/utqiagvik_10m_UTM4N_seaice_mask_0and1.tif --datum WGS84 --individually-normalize --normalize-ip-tiles --ip-per-tile 50 --matches-per-tile 10 --min-triangulation-angle 10 --mapproj-dem /Users/ben/Dropbox/UW_Shean/COP/COP30_utqiagvik_lzw-adj_proj.tif --propagate-errors --tri-weight 0.1 --tri-robust-threshold 0.1 --camera-weight 0 1040010074793300.r100.tif 1040010075633C00.r100.tif 1040010074793300.r100.xml 1040010075633C00.r100.xml -o ba/ba_50ips_10matches_dg_weight_image --threads 28

Stereo (3 hours and 41 minutes):

stereo --stereo-algorithm asp_mgm --corr-kernel 7 7 --subpixel-kernel 15 15 --cost-mode 4 --subpixel-mode 9 --corr-max-levels 5 --filter-mode 1 --erode-max-size 0 --individually-normalize --corr-memory-limit-mb 5000 --sgm-collar-size 256 --corr-tile-size 1024 --alignment-method none --corr-seed-mode 1 --compute-point-cloud-center-only --threads 24 1040010074793300_ortho_0.35m.tif 1040010075633C00_ortho_0.35m.tif ba/ba_50ips_10matches_dg_weight_image-1040010074793300.r100.adjusted_state.json ba/ba_50ips_10matches_dg_weight_image-1040010075633C00.r100.adjusted_state.json stereo_ba_50ips_10matches_dg_weight_image__ortho_0.35m_mode_asp_mgm_spm_9_corr_7_rfne_15_cost_4_refdem_COP30/20220417_2252_1040010074793300_1040010075633C00 /Users/ben/Dropbox/UW_Shean/COP/COP30_utqiagvik_lzw-adj_proj.tif

Point2dem (0 hours and 26 minutes):

point2dem --nodata-value -9999 --t_srs EPSG:32604 --threads 24 --propagate-errors --remove-outliers --remove-outliers-params 75.0 3.0 --errorimage --tr 1 -o stereo_ba_50ips_10matches_dg_weight_image__ortho_0.35m_mode_asp_mgm_spm_9_corr_7_rfne_15_cost_4_refdem_COP30/20220417_2252_1040010074793300_1040010075633C00_1m stereo_ba_50ips_10matches_dg_weight_image__ortho_0.35m_mode_asp_mgm_spm_9_corr_7_rfne_15_cost_4_refdem_COP30/20220417_2252_1040010074793300_1040010075633C00-PC.tif

3. Scene Plots

Visualize the input satellite imagery (orthorectified scenes) used for stereo processing. This helps verify image quality and overlap.

import os
from asp_plot.scenes import ScenePlotter

# Optional: create directory to save plots

# plots_directory = os.path.join(directory, "asp_plots")
# os.makedirs(plots_directory, exist_ok=True)

plotter = ScenePlotter(
  directory,
  stereo_directory,
  title="Input Scenes"
)

plotter.plot_scenes(
  # Optional parameters to set for saving (requires plots_directory to be created above):

  # save_dir=plots_directory,
  # fig_fn="stereo_scenes.png"
)

4. Scene Geometry Plots

Analyze stereo acquisition geometry from XML camera metadata. Displays:

• Skyplot: Satellite viewing angles (azimuth and elevation)
• Map view: Scene footprints and overlap area
• Stereo metrics: Convergence angle, base-to-height ratio, asymmetry

from asp_plot.stereo_geometry import StereoGeometryPlotter

geometry_plotter = StereoGeometryPlotter(
  directory
)

geometry_plotter.dg_geom_plot(
    # Optional parameters to set for saving (requires plots_directory to be created above):

    # save_dir=plots_directory,
    # fig_fn="scene_geometry.png"
)

5. Bundle Adjustment Plots

Visualize bundle adjustment optimization results:

• Residual maps: Initial vs. final camera optimization residuals
• Geodiff: Triangulated point differences vs. reference DEM
• Map-projected residuals: Interest point separation on reference DEM surface

Bundle adjustment refines camera models to improve stereo matching consistency.

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

map_crs = "EPSG:32604"  # UTM Zone 4N

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

ba_files = ReadBundleAdjustFiles(directory, ba_directory)
resid_initial_gdf, resid_final_gdf = ba_files.get_initial_final_residuals_gdfs(residuals_in_meters=True)
resid_mapprojected_gdf = ba_files.get_mapproj_residuals_gdf()
resid_triangulation_uncert_df = ba_files.get_propagated_triangulation_uncert_df()

resid_triangulation_uncert_df

	left_image	right_image	horiz_error_median	vert_error_median	horiz_error_mean	vert_error_mean	horiz_error_stddev	vert_error_stddev	num_meas
0	1040010074793300.r100.tif	1040010075633C00.r100.tif	3.181445	8.03649	3.182249	8.036751	0.012553	0.001756	7436

Bundle Adjustment Residuals (Log Scale)

Log scale visualization emphasizes patterns across wide value ranges.

plotter = PlotBundleAdjustFiles(
  [resid_initial_gdf, resid_final_gdf],
  lognorm=True,
  title="Bundle Adjust Initial and Final Residuals (Log Scale)"
)

plotter.plot_n_gdfs(
    column_name="mean_residual",
    cbar_label="Mean residual (px)",
    map_crs=map_crs,
    **ctx_kwargs
)

plotter.plot_n_gdfs(
    column_name="mean_residual_meters",
    cbar_label="Mean residual (m)",
    map_crs=map_crs,
    **ctx_kwargs
)

Bundle Adjustment Residuals (Linear Scale)

Linear scale shows absolute residual magnitudes more clearly.

plotter.lognorm = False
plotter.title = "Bundle Adjust Initial and Final Residuals (Linear Scale)"

plotter.plot_n_gdfs(
    column_name="mean_residual",
    cbar_label="Mean residual (px)",
    common_clim=False,
    map_crs=map_crs,
    **ctx_kwargs
)

plotter.plot_n_gdfs(
    column_name="mean_residual_meters",
    cbar_label="Mean residual (m)",
    common_clim=False,
    map_crs=map_crs,
    **ctx_kwargs
)

plotter.title = "Bundle Adjust Initial and Final Residuals (Common Linear Scale)"

plotter.plot_n_gdfs(
    column_name="mean_residual",
    cbar_label="Mean residual (px)",
    map_crs=map_crs,
    **ctx_kwargs
)

Bundle Adjustment Observation Statistics

Number of observations per point and height above datum.

plotter = PlotBundleAdjustFiles(
  [resid_final_gdf],
)

plotter.title = "Bundle Adjust Residual Number of Observations"

plotter.plot_n_gdfs(
    column_name="num_observations",
    cbar_label="Number of Observations",
    map_crs=map_crs,
    **ctx_kwargs
)

plotter.title = "Bundle Adjust Residuals Height Above Datum"

plotter.plot_n_gdfs(
    column_name="height_above_datum",
    cbar_label="Height above datum (m)",
    map_crs=map_crs,
    **ctx_kwargs
)

Map-Projected Residuals

Distance between interest points when projected onto the reference DEM surface.

plotter = PlotBundleAdjustFiles(
  [resid_mapprojected_gdf],
  title="Bundle Adjust Midpoint distance between\nfinal interest points projected onto reference DEM",
)

plotter.plot_n_gdfs(
    column_name="mapproj_ip_dist_meters",
    cbar_label="Interest point distance (m)",
    map_crs=map_crs,
    **ctx_kwargs
)

Geodiff vs. Reference DEM

Compares triangulated points (before and after bundle adjustment) to the reference DEM surface.

geodiff_initial_gdf, geodiff_final_gdf = ba_files.get_initial_final_geodiff_gdfs()

plotter = PlotBundleAdjustFiles(
  [geodiff_initial_gdf, geodiff_final_gdf],
  lognorm=False
)

plotter.title = "Bundle Adjust Initial and Final Geodiff vs. Reference DEM"

plotter.plot_n_gdfs(
    column_name="height_diff_meters",
    cbar_label="Height difference (m)",
    map_crs=map_crs,
    cmap="RdBu",
    clim=(-1, 1),
    **ctx_kwargs
)

plotter.title = "Bundle Adjust Initial and Final Geodiff vs. Reference DEM (Auto Scale)"

plotter.plot_n_gdfs(
    column_name="height_diff_meters",
    cbar_label="Height difference (m)",
    map_crs=map_crs,
    cmap="RdBu",
    symm_clim=True,
    **ctx_kwargs
)

6. Stereo Processing Results

Visualize the primary outputs of stereo processing:

• Hillshade: Shaded relief with intersection error overlays
• Match points: Interest point locations used for stereo correlation
• Disparity: Pixel offsets (in meters and pixels) between left/right images
• DEM results: Elevation, intersection error, and differences vs. reference DEM

from asp_plot.stereo import StereoPlotter

# Reference DEM can be passed to the StereoPlotter, but if it is not, it will
# be searched for in the stereo logs. If it's not found there, you will be warned,
# and no difference plots will be generated.

# reference_dem = "/Users/ben/Dropbox/UW_Shean/COP/COP30_utqiagvik_lzw-adj_proj.tif"

plotter = StereoPlotter(
  directory, 
  stereo_directory,
  # Optional args:
  # dem_fn="my-custom-dem-name.tif",
  dem_gsd=1,
  # reference_dem=reference_dem,
)

Reference DEM: /Users/ben/Dropbox/UW_Shean/COP/COP30_utqiagvik_lzw-adj_proj.tif


ASP DEM: /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/20220417_2252_1040010074793300_1040010075633C00-DEM_1m.tif

Hillshade with Details

Detailed hillshade with intersection error percentiles shown in inset.

plotter.title = "Hillshade with details"

plotter.plot_detailed_hillshade(
  intersection_error_percentiles=[16, 50, 84],
  subset_km=1,
)

Stereo Match Points

Locations of interest points used for stereo correlation.

plotter.title="Stereo Match Points"

plotter.plot_match_points()

Disparity (Meters)

Pixel offsets between stereo images, converted to ground distance in meters.

plotter.title = "Disparity (meters)"

plotter.plot_disparity(
  unit="meters",
  quiver=True,
)

Disparity (Pixels)

Pixel offsets between stereo images in original image pixel coordinates.

plotter.title = "Disparity (pixels)"

plotter.plot_disparity(
  unit="pixels",
  quiver=True,
)

DEM Results Summary

Multi-panel plot showing elevation, intersection error, and difference vs. reference DEM.

plotter.title = "Stereo DEM Results"

plotter.plot_dem_results()

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

WARNING:asp_plot.stereo:

Found a DEM of difference: /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/20220417_2252_1040010074793300_1040010075633C00-DEM_1m_COP30_utqiagvik_lzw-adj_proj_diff.tif.
Using that for difference map plotting.

# You can also use custom color limits for the DEM results plots:

# plotter.title = "Stereo DEM Results (forced color limits)"

# plotter.plot_dem_results(
#   el_clim=(0, 10),
#   ie_clim=(0, 0.2),
#   diff_clim=(-5, 5),
# )

7. ICESat-2 Altimetry Validation

Compare the ASP DEM with ICESat-2 ATL06-SR altimetry data from NASA's SlideRule API.

This section:

1. Requests ATL06-SR data for the DEM extent
2. Filters by processing level (all, ground, canopy, top-of-canopy)
3. Applies temporal filters around the acquisition date
4. Performs DEM-to-altimetry alignment using ASP's pc_align
5. Generates comparison plots and statistics

Note: Requires internet connection for SlideRule API access.

import glob
from asp_plot.altimetry import Altimetry

dem_fn = glob.glob(os.path.join(directory, "stereo*/*DEM_1m.tif"))[0]
aligned_dem_fn = glob.glob(os.path.join(directory, "stereo*/*DEM_1m*pc_align*.tif"))[0]

map_crs = "32604"  # UTM Zone 4N

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

icesat = Altimetry(
  directory=directory,
  dem_fn=dem_fn,
  # If the aligned_dem_fn is not provided, it will be calculated automatically
  # using ASP's pc_align tool, if it is available in your PATH.
  aligned_dem_fn=aligned_dem_fn
)

Request ATL06-SR Data

Fetch ICESat-2 data from SlideRule for multiple processing levels. Data is cached locally as parquet files.

icesat.request_atl06sr_multi_processing(
    save_to_parquet=True,
)

# There are many more custom processing options available.
# See SlideRule's Parameter documentation for details:
# https://slideruleearth.io/web/rtd/user_guide/icesat2.html#parameters

# icesat.request_atl06sr_multi_processing(
#     res=10,
#     len=20,
#     ats=20,
#     cnt=5,
#     maxi=5,
#     save_to_parquet=True,
#     processing_levels=["ground"],
# )

ICESat-2 ATL06 request processing for: all
{'poly': [{'lon': -156.83342896997698, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.26658548886508}], 'res': 20, 'len': 40, 'ats': 20, 'maxi': 6, 'samples': {'esa_worldcover': {'asset': 'esa-worldcover-10meter'}}, 'cnf': 'atl03_high', 'srt': -1, 'cnt': 10}
Existing file found, reading in: atl06sr_all.parquet
Filtering ATL06-SR all

ICESat-2 ATL06 request processing for: ground
{'poly': [{'lon': -156.83342896997698, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.26658548886508}], 'res': 20, 'len': 40, 'ats': 20, 'maxi': 6, 'samples': {'esa_worldcover': {'asset': 'esa-worldcover-10meter'}}, 'cnf': 'atl03_low', 'srt': -1, 'cnt': 5, 'atl08_class': 'atl08_ground'}

WARNING:sliderule.sliderule:Using simplified polygon (for CMR request only!), 5 points using tolerance of 0.0001

Filtering ATL06-SR ground

ICESat-2 ATL06 request processing for: canopy
{'poly': [{'lon': -156.83342896997698, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.26658548886508}], 'res': 20, 'len': 40, 'ats': 20, 'maxi': 6, 'samples': {'esa_worldcover': {'asset': 'esa-worldcover-10meter'}}, 'cnf': 'atl03_medium', 'srt': -1, 'cnt': 5, 'atl08_class': 'atl08_canopy'}

WARNING:sliderule.sliderule:Using simplified polygon (for CMR request only!), 5 points using tolerance of 0.0001

Filtering ATL06-SR canopy

ICESat-2 ATL06 request processing for: top_of_canopy
{'poly': [{'lon': -156.83342896997698, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.39688565394232}, {'lon': -156.40964712543294, 'lat': 71.26658548886508}, {'lon': -156.83342896997698, 'lat': 71.26658548886508}], 'res': 20, 'len': 40, 'ats': 20, 'maxi': 6, 'samples': {'esa_worldcover': {'asset': 'esa-worldcover-10meter'}}, 'cnf': 'atl03_medium', 'srt': -1, 'cnt': 5, 'atl08_class': 'atl08_top_of_canopy'}

WARNING:sliderule.sliderule:Using simplified polygon (for CMR request only!), 5 points using tolerance of 0.0001

Filtering ATL06-SR top_of_canopy

print("There are", icesat.atl06sr_processing_levels["all"].shape[0], "ATL06SR points available.")
print("There are", icesat.atl06sr_processing_levels_filtered["all"].shape[0], "filtered ATL06SR points available.")

There are 163350 ATL06SR points available.
There are 145856 filtered ATL06SR points available.

Filter ATL06-SR Data

Apply land cover and temporal filters to improve comparison quality.

# Filter by ESA WorldCover to remove water points
icesat.filter_esa_worldcover(filter_out="water")

print("After filtering out water, there are", icesat.atl06sr_processing_levels_filtered["ground"].shape[0], "filtered ground ATL06SR points available.")

After filtering out water, there are 45309 filtered ground ATL06SR points available.

# Apply temporal filters around acquisition date (15/45/91 day windows + seasonal)
icesat.predefined_temporal_filter_atl06sr()

# You can also create user defined temporal filters with `generic_temporal_filter_atl06sr`:

# icesat.generic_temporal_filter_atl06sr(
#     select_years=[2021, 2022, 2023],
#     select_months=[1, 2, 3, 10, 11, 12],
#     select_days=[1, 2, 3, 4, 5, 6, 7, 27, 28, 29, 30, 31]
# )

Filtering ATL06 with 15 day pad, 45 day, 91 day pad, and seasonal pad around 2022-04-17 22:52:18.495475 for: all

Filtering ATL06 with 15 day pad, 45 day, 91 day pad, and seasonal pad around 2022-04-17 22:52:18.495475 for: ground

Filtering ATL06 with 15 day pad, 45 day, 91 day pad, and seasonal pad around 2022-04-17 22:52:18.495475 for: canopy

Filtering ATL06 with 15 day pad, 45 day, 91 day pad, and seasonal pad around 2022-04-17 22:52:18.495475 for: top_of_canopy

print("There are", icesat.atl06sr_processing_levels_filtered["ground_seasonal"].shape[0], "filtered ground seasonal ATL06SR points available.")

There are 12505 filtered ground seasonal ATL06SR points available.

ICESat-2 Temporal Distribution

Visualize the temporal distribution of ATL06-SR points colored by acquisition date.

icesat.plot_atl06sr_time_stamps(
   key="ground",
   figsize=(15, 10),
   save_dir=None,
   fig_fn=None,
   **ctx_kwargs,
)

icesat.plot_atl06sr_time_stamps(
   key="canopy",
   figsize=(15, 10),
   save_dir=None,
   fig_fn=None,
   **ctx_kwargs,
)

ICESat-2 Spatial Distribution

Map view of ATL06-SR points colored by elevation and laser spot.

icesat.plot_atl06sr(
    key="ground_seasonal",
    map_crs=map_crs,
    cmap="inferno",
    plot_dem=False,
    **ctx_kwargs
)

icesat.plot_atl06sr(
    key="ground_45_day_pad",
    map_crs=map_crs,
    cmap="inferno",
    plot_dem=False,
    plot_beams=True,
    **ctx_kwargs
)

icesat.plot_atl06sr(
    key="canopy_seasonal",
    map_crs=map_crs,
    cmap="inferno",
    plot_dem=True,
    plot_beams=True,
    **ctx_kwargs
)

DEM vs. ICESat-2 Comparison

Map view showing DEM-altimetry differences before and after alignment.

icesat.mapview_plot_atl06sr_to_dem(
    key="ground_15_day_pad",
    **ctx_kwargs,
)

icesat_minus_dem not found in ATL06 dataframe: ground_15_day_pad. Running differencing first.

# Now for the pc_align aligned DEM version
icesat.mapview_plot_atl06sr_to_dem(
    key="ground_15_day_pad",
    plot_aligned=True,
    **ctx_kwargs,
)

# Show a quick histogram
icesat.histogram(
    key="ground_seasonal",
    plot_aligned=True,
)

PC Alignment Report

Perform pc_align with multiple temporal filters and compare results. This helps assess alignment consistency and chooses the best temporal filter.

icesat.alignment_report(
    processing_level="ground",
    minimum_points=500,
    agreement_threshold=0.25,
    write_out_aligned_dem=True,
    min_translation_threshold=0.1,
    key_for_aligned_dem="ground_15_day_pad",
)

Warning: Translation components vary by more than 25.0% across temporal filters. The translation applied to the aligned DEM may be inaccurate.

X shift range: 0.881 m (mean: 0.359 m)
Y shift range: 0.965 m (mean: 0.017 m)
Z shift range: 0.336 m (mean: -0.649 m)

Writing out: /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/20220417_2252_1040010074793300_1040010075633C00-DEM_1m_pc_align_translated.tif


Wrote out ground_15_day_pad aligned DEM to /Users/ben/Dropbox/UW_Shean/WV/2022/WV03_20220417_1040010074793300_1040010075633C00/stereo/20220417_2252_1040010074793300_1040010075633C00-DEM_1m_pc_align_translated.tif

icesat.alignment_report_df

	key	p16_beg	p50_beg	p84_beg	p16_end	p50_end	p84_end	x_shift	y_shift	z_shift	translation_magnitude
0	ground	0.393706	0.635813	0.840236	0.046457	0.153485	0.336891	0.223944	-0.068757	-0.593812	0.638350
1	ground_15_day_pad	0.641055	0.773392	0.921410	0.028533	0.094267	0.210883	0.155468	0.041450	-0.763709	0.780474
2	ground_45_day_pad	0.585919	0.759794	0.914996	0.033228	0.110376	0.247181	0.984863	0.071390	-0.477306	1.096756
3	ground_91_day_pad	0.560620	0.747609	0.931639	0.037673	0.127983	0.280320	0.103808	-0.461917	-0.813125	0.940912
4	ground_seasonal	0.566514	0.732346	0.886025	0.032193	0.108120	0.234966	0.324831	0.503479	-0.597760	0.846359

icesat.mapview_plot_atl06sr_to_dem(
    key="ground_15_day_pad",
    plot_aligned=True,
    **ctx_kwargs,
)

# Show a quick histogram
icesat.histogram(
    key="ground_seasonal",
    plot_aligned=True,
)

WIP: Profile plots

Planned functionality

# Collect only the coincident filtereded data again for profile plotting
# icesat.filter_atl06sr(
#     h_sigma_quantile=0.95,
#     mask_worldcover_water=True,
#     save_to_csv=False,
#     select_months=[4],
#     select_years=[2022],
# )

# icesat.plot_atl06sr(
#     title=f"Cleaned beam strengths (n={icesat.atl06sr_filtered.shape[0]})",
#     filtered=True,
#     plot_beams=True,
#     plot_dem=False,
#     map_crs=map_crs,
#     **ctx_kwargs
# )

# icesat.plot_atl06sr_dem_profiles(title="Profiles", only_strong_beams=True)

8. CSM Camera Model Analysis

Compare original and optimized CSM (Community Sensor Model) camera models from bundle_adjust or jitter_solve.

Visualizes:

• Position differences (X, Y, Z) along the satellite orbit
• Orientation angle differences (roll, pitch, yaw)
• Camera footprints on a map

Note: This example uses different test data (Salar de Uyuni, Bolivia) to demonstrate jitter correction.

from asp_plot.csm_camera import csm_camera_summary_plot

map_crs = "32619" # UTM 19S (for Uyuni)
title = "Jitter correction (Uyuni), Less Constrained"

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

# First set of cameras (scene 1)
original_camera = "../../tests/test_data/jitter/uyuni/csm-104001001427B900.r100.adjusted_state.json"
optimized_camera = "../../tests/test_data/jitter/uyuni/jitter_solved_run-csm-104001001427B900.r100.adjusted_state.json"
cam1_list = [original_camera, optimized_camera]

# Second set of cameras (scene 2)
original_camera = "../../tests/test_data/jitter/uyuni/csm-1040010014761800.r100.adjusted_state.json"
optimized_camera = "../../tests/test_data/jitter/uyuni/jitter_solved_run-csm-1040010014761800.r100.adjusted_state.json"
cam2_list = [original_camera, optimized_camera]

CSM Camera Comparison (Basic)

Basic comparison without trimming or scaling.

csm_camera_summary_plot(
    cam1_list,
    cam2_list,
    map_crs=map_crs,
    title=title,
    trim=False,
    figsize=(20, 15),
    add_basemap=True,
    **ctx_kwargs
)

CSM Camera Comparison (Advanced)

With trimming, log scaling, and percentile-based color limits to emphasize patterns.

csm_camera_summary_plot(
    cam1_list,
    cam2_list,
    map_crs=map_crs,
    title=title,
    trim=True,
    shared_scales=True,
    log_scale_positions=True,
    log_scale_angles=True,
    upper_magnitude_percentile=95,
    figsize=(20, 15),
    add_basemap=True,
    **ctx_kwargs
)

Single Camera Comparison

You can also analyze a single satellite by passing only cam1_list.

csm_camera_summary_plot(
    cam1_list,
    map_crs=map_crs,
    title=title,
    trim=True,
    shared_scales=True,
    log_scale_positions=True,
    log_scale_angles=True,
    upper_magnitude_percentile=95,
    figsize=(20, 15),
    add_basemap=True,
    **ctx_kwargs
)