Mars Examination foR ClimAte Patterns
This analysis uses Mars Climate Sounder observations and known dust storm locations from the Mars Dust Storm Sequence Dataset (MDSSD) and Mars Dust Activity Database (MDAD).
MERCAP harmonizes atmospheric observations from MCS with MDSSD/MDAD storm catalogs. This framework enables targeted extraction and analysis of atmospheric data (dust opacity, temperature, and surface temperature) associated with specific dust storm instances. After extracting coincident atmospheric measurements, we can use that data to study how local and regional dust storms impact the Martian atmosphere and surface.
Further, we can use storm-free "control" measurements from the same location, season, and local time (pooled across Mars years) to estimate how much any given storm event perturbed the atmosphere.
MERCAP v2 replaces the per-storm database queries with static, hive-partitioned Parquet catalogs generated once by the numbered pipeline stages (mercap2/pipeline/). "Hive-partitioned" means partition-column values are encoded in the directory structure (e.g., ddr2/lat_bin=-30/mars_year=29/part-*.parquet); tools like pandas, pyarrow, dask, and duckdb recognize the key=value pattern and skip non-matching directories entirely when filtering, so reads touch only the spatial/temporal chunks needed. See the pyarrow partitioned-dataset docs for more.
Note on longitude conventions: all v2 catalogs and helpers use east-positive longitude on [-180, 180). Inputs on east-positive [0, 360) are convertible (mercap2/utils/gridding.py); west-positive longitudes (used by some products, e.g. raw MARCI) are not auto-detectable and, if desired when extending this codebase, must be converted to east-positive at ingest.
Create a conda environment, install the two git-only dependencies, then install mercap2 itself:
conda create -n mercap2 python=3.14
conda activate mercap2
git clone https://github.com/cloudspotting-on-mars/mcstools
pip install -e mcstools
git clone https://github.com/marek-slipski/mars_time
pip install -e mars_time
pip install -e .
Optional extras: pip install -e '.[dev]' (pytest, jupyter) and pip install -e '.[viz]' (cartopy etc. for the plotting scripts). All other dependencies are declared in pyproject.toml and installed automatically.
The pipeline is comprised of six numbered click CLIs in mercap2/pipeline/, run in order. Each stage reads the previous stage's outputs and writes its own under config.DATA_ROOT. Check out the docs/source/_static/diagrams/pipeline_diagram.md file for the data-flow diagram.
Conventions shared by every stage:
- Defaults live in
mercap2/config.py, grouped by stage. CLI options exist for paths and the most common knobs to turn, but grid sizes, exclusion rules, statistics, dtypes, etc. are config constants --smoke_testruns the stage end-to-end on a small slice of data (one flat file, one partition, one storm instance, ...). Use it as a sanity check before a full-scale run.--log_leveland--log_fpathcontrol console verbosity and optional file logging;--n_jobssets the joblib worker count (defaultconfig.N_JOBS).- Stages that write catalogs take
--existing_output {fail,append,clobber}(defaultfailto prevent overwriting) to control what happens when output already exists.
Parses the raw MDAD (NetCDF masks + master CSV) and MDSSD (IDL files) datasets into one GeoParquet catalog each: one row per storm instance (one storm outline on one MDGM) with the outline polygon, centroid, confidence level, and Mars timing columns. If config.IMPROVE_MARCI_TIMING_PRECISION is enabled (default), the MARCI options substitute more precise swath-level Ls/datetime values from the MARCI cumulative index (available from the PDS).
CURRENT_TIME=$(date "+%Y%m%d_%H%M%S") && \
python mercap2/pipeline/s00_parse_storm_catalogs.py \
--mdssd_head_dir /path/to/mdssd \
--mdad_head_dir /path/to/mdad \
--mdad_csv_fpath /path/to/mdad/mdad_master.csv \
--marci_lists_dir /path/to/marci_lists \
--marci_cumindex_fpath /path/to/marci/CUMINDEX.TAB \
--log_fpath /path/to/mercap2/storm_catalogs/s00_parse_storm_catalogs_${CURRENT_TIME}.log \
--write_csv
# Debugging
CURRENT_TIME=$(date "+%Y%m%d_%H%M%S") && \
python mercap2/pipeline/s00_parse_storm_catalogs.py \
--catalogs 'mdssd' \
--mdssd_head_dir /path/to/mdssd \
--marci_lists_dir /path/to/marci_lists \
--marci_cumindex_fpath /path/to/marci/CUMINDEX.TAB \
--n_jobs 1
Pull L2 MCS data. Traverses the MCS L2 flat-file tree (config.MCS_DATA_DIR, read via mcstools) and writes the DDR1 (per-profile, ~10M rows) and DDR2 (per-profile-level, ~1B rows) Parquet catalogs, applying the v1 validity conventions (GOOD_OBS_QUALS, DROP_INVALID_DDR1, MIN_DUST_PERMITTED) with logged drop counts. Validate with --smoke_test and/or a --start_date/--end_date slice before committing to the full multi-day run.
# Test it out
python mercap2/pipeline/s01_extract_catalogs.py --mcs_data_dir /path/to/mcs_flat_files --smoke_test
# Run full version
CURRENT_TIME=$(date "+%Y%m%d_%H%M%S") && \
python mercap2/pipeline/s01_extract_catalogs.py \
--mcs_data_dir /path/to/mcs_flat_files \
--log_fpath /path/to/mercap2/mcs_catalogs/s01_extract_catalogs_${CURRENT_TIME}.log \
--n_jobs 8 \
--existing_output clobber
Identify which profiles are disqualified from serving as storm-free controls: any profile whose coordinates fall inside a storm polygon (at/above min_conflev, in either catalog, optionally spatially buffered) within sol_buffer sols of that storm snapshot. Output is a single compact Parquet table, one row per profile. Knobs live in config.STORM_EXCLUSION.
CURRENT_TIME=$(date "+%Y%m%d_%H%M%S") && \
python mercap2/pipeline/s02_flag_control_exclusions.py \
--log_fpath /path/to/mercap2/storm_overlap_flags/s02_flag_storm_overlap_${CURRENT_TIME}.log
Aggregates the non-excluded ("control") measurements into background distributions per grid cell (lat_bin, lon_bin, ls_bin, ltst_window[, level]): counts, mean/std, median/MAD/IQR (raw and sigma-scaled), a quantile ladder, and year-balanced center statistics so a heavily sampled Mars year can't dominate a cell. Two streams share the grid: ddr2 (T, Dust per level) and ddr1 (surface quantities, no level axis). Cells failing min_count/min_mars_years are kept but marked invalid. Grid spacing, variables, and thresholds live in config.CLIMATOLOGY.
CURRENT_TIME=$(date "+%Y%m%d_%H%M%S") && \
python mercap2/pipeline/s03_build_climatology_cube.py \
--log_fpath /path/to/mercap2/climatology_cube/s03_build_climatology_cube_${CURRENT_TIME}.log
For every measurement, interpolates a background from the cube (within its spatial region, LTST window, and level) and writes a long-format perturbation catalog: raw value, background center/scale, delta, robust_z (delta / sigma_mad), empirical percentile_rank against the cell's quantile ladder, and background count/validity. Interpolation methods and background column choices live in config.PERTURBATION.
CURRENT_TIME=$(date "+%Y%m%d_%H%M%S") && \
python mercap2/pipeline/s04_compute_perturbations.py \
--existing_output clobber \
--log_fpath /path/to/mercap2/perturbation_catalogs/s04_compute_perturbations_${CURRENT_TIME}.log \
--n_jobs 8
For each storm instance, selects the perturbation rows inside the storm polygon, places each on a storm-relative local sol (see rel_sol_method in config.AGGREGATION; nighttime passes get a +0.5 sol reference offset since storms are detected in daytime MARCI images), and writes per-sol aggregates over sol_window_bounds (default -10 to +11) sols: the raw mean/median value, mean delta / robust z / percentile rank, tail fractions, and counts per (storm instance, LTST window, variable[, level], relative sol). The raw mean_value/median_value pool every selected row regardless of background validity (e.g., raw dust column opacity, a raw reading needs no climatology comparison to be meaningful); every other mean and the tail fractions cover only measurements with a valid background.
CURRENT_TIME=$(date "+%Y%m%d_%H%M%S") && \
python mercap2/pipeline/s05_aggregate_per_instance_data.py \
--n_jobs 16 \
--existing_output clobber \
--write_csv
All products are written under config.DATA_ROOT and read back with mercap2.utils.catalog_io.read_catalog (or pandas/pyarrow/etc. directly). Exact coordinates and times are always stored as columns; partition bins are a storage layout only, intentionally coarser than the analysis grid.
| Catalog | Stage | One row per | Layout |
|---|---|---|---|
storm_catalogs/ |
s00 | storm instance (storm x MDGM) | one GeoParquet file per catalog (mdad/mdssd) |
mcs_catalogs/ddr1/ |
s01 | MCS profile (~10M) | hive: lat_bin/mars_year |
mcs_catalogs/ddr2/ |
s01 | profile-level (~1B) | hive: lat_bin/mars_year |
storm_overlap_flags/ |
s02 | DDR1 profile | single Parquet file |
climatology_cube/{ddr1,ddr2}/ |
s03 | grid cell x variable | hive: variable/lat_bin |
perturbation_catalogs/{ddr1,ddr2}/ |
s04 | measurement x variable | hive: variable/lat_bin/mars_year |
storm_aggregates/{ddr1,ddr2}/ |
s05 | storm instance x LTST window x variable [x level] x relative sol | hive: storm_catalog/variable |
DDR2 rows carry duplicated per-profile metadata (coordinates, l_s, LTST, mars_year, sol, ltst_window). This is space inefficient but downstream stages never need a 1B-row join back to DDR1. Rows are sorted by l_s within files so Parquet row-group statistics prune seasonal windows at read time. Example read:
from mercap2 import config
from mercap2.utils.catalog_io import read_catalog
df = read_catalog(config.DDR2_CATALOG_DIR, columns=['T', 'Pres', 'l_s', 'Surf_lat', 'Surf_lon'],
filters=[('lat_bin', '=', -30), ('mars_year', '=', 29)])Partition bin sizes were determined experimentally with scripts/benchmark_partitions.py; rerun it if the read workloads change substantially and then tune the config.
mercap2/config.py is the primary location for tunable parameters, grouped by stage: catalog paths, PARTITIONING (hive bin sizes), LTST_WINDOWS (named day/night local-time windows), STORM_EXCLUSION (control disqualification rules), CLIMATOLOGY (grid spacing, variables, statistics, validity thresholds), PERTURBATION (interpolation and background column choices), and AGGREGATION (sol window, relative-sol method, tail percentiles). Prefer changing values there over adding per-script constants, and update the paths at the top of the config to match your local data locations before running the pipeline.
mercap2/viz/ holds the plotting CLIs, ported from mercap v1 onto the v2 Parquet catalogs. Two shared modules replace the filtersheet:
mercap2.utils.profile_search: single-storm-instance selection against DDR1/DDR2/s04 catalogs.get_storm_instanceresolves one storm row (by its positionalstorm_instance_id, e.g.'mdssd-006518', or by the natural key(storm_catalog, mdgm, storm_id, conflev, seq_id));select_rows_for_storm_instanceselects the measurements inside that instance's polygon and places them on its local-sol timeline.mercap2.viz.storm_plot_data_prep: tidy per-instance frames built from the s00 storm metadata and the s05 aggregates (load_storm_instance_metadata,load_aggregate_series,load_instance_coverage,join_instance_table,pivot_sol_columns), consumed bymercap2.utils.filtering.filter_storms(the metadata/coverage/area/lifecycle filters, also ported from v1) and bymercap2.viz.timeseries_utils(the reshape/pad/coverage logic shared by the two timeseries scripts below).
Every script below shares those loaders' --stream/--storm_catalog, --variable, --metric, and --ltst_window (--stream selects ddr1/ddr2; --storm_catalog selects mdad/mdssd) and takes the same metadata/coverage filter options (--conflev_min, --lat_min/--lat_max, --ls_min/--ls_max, --storm_len_max, --n_profiles_min, --area_min/--area_max, --storm_lifecycle_filter). Scripts fixed to one stream (the two timeseries scripts below) take no --stream. --area_min/--area_max bound storm area (km^2), inclusive; since area is always non-negative, pass -1 to either one to disable that bound. Plotting constants (figure sizes, DPI, fonts, colors, the DDR2 level label template, ...) live in mercap2/config_plotting.py.
Per-latitude-band time series of one s05 metric, MDAD and MDSSD side by side, with a shaded quantile band and a per-sol storm-count histogram. Always reads the DDR1 stream (--variable is limited to its surface variables): latitude bands claim the figure's rows, so a per-level figure belongs to the by-height script below.
python mercap2/viz/plot_timeseries_per_lat_bin_v3.py \
--variable Dust_column --metric mean_delta \
--sol_range -10,11 --output_path /path/to/mercap2/jgr_plots/dust_column_timeseries.png
The same per-latitude-band layout, but one row per DDR2 retrieval level instead of one row per latitude band (always reads the DDR2 stream, since levels are its retrieval axis). Each row is labeled with that level's pressure and altitude derived from the s05 aggregates' own mean_pres/mean_alt over the population plotted, so the altitudes follow the latitude band and storm selection (as opposed to a fixed table).
python mercap2/viz/plot_timeseries_by_height_v3.py \
--variable T --metric mean_delta --levels 30,24,18 \
--sol_range -20,21 --output_dir /path/to/mercap2/jgr_plots/timeseries_by_height
A before/during/after temperature (left column) and dust (right column) profile comparison for one storm instance: every raw storm profile in each relative-sol window, overlaid on a control distribution (--control_mode shaded_fill/density/lines). "Control" is s04's per-measurement background_center, attached 1:1 to each storm profile -- see the module docstring's porting notes. --plot_quantity value (default) plots the raw reading against its control background; --plot_quantity delta plots s04's already-computed perturbation instead and drops the control panel in favor of a zero-perturbation reference line, since a background's delta from itself is trivially zero. --add_geospatial_inset overlays a small cartopy globe (storm outline, center, and the "during" row's profile locations) on the middle row's dust panel.
python mercap2/viz/plot_enhanced_profile_comparison_v2.py \
--storm_instance_id mdssd-006518 \
--sol_range_before -12,-1 --sol_range_during -1,6 --sol_range_after 6,17 \
--control_mode shaded_fill --plot_quantity value --yaxis Alt --sharex --sharey \
--save_path_plot /path/to/mercap2/jgr_plots/profile_comparison_mdssd-006518.png \
--save_path_csv /path/to/mercap2/jgr_plots/profile_comparison_mdssd-006518.csv
To plot the storm's perturbation (delta from its s04 background) instead of the raw reading:
python mercap2/viz/plot_enhanced_profile_comparison_v2.py \
--storm_instance_id mdssd-006518 \
--sol_range_before -12,-1 --sol_range_during -1,6 --sol_range_after 6,17 \
--plot_quantity delta --yaxis Alt --sharex --sharey \
--save_path_plot /path/to/mercap2/jgr_plots/profile_comparison_delta_mdssd-006518.png
A 2x3 "collage" for one storm instance: a cartopy geospatial overview plus DDR1 dust/temperature perturbation traces on top, and DDR2 control/storm/difference dust heatmaps (pressure or altitude vs. relative sol) on the bottom, all restricted to the same background_valid rows so the difference panel is exactly the other two panels' difference. --add_inset draws a near-surface inset over the bottom-row heatmaps; --marci_fpath optionally overlays a MARCI RGB color-mosaic crop (the image the storm was originally traced from) in the geospatial subplot's corner.
python mercap2/viz/plot_collage_light_jgr2025_v3.py \
--storm_instance_id mdssd-006518 --output_fpath /path/to/mercap2/jgr_plots/collage_mdssd-006518.png \
--level_axis altitude --add_inset
The full 5x5, 22-panel storm collage (a rewrite of v1's gen_plot_collage_1 against the static v2 catalogs): the same geospatial overview and metadata block as the light collage, T/Dust storm/control/diff/robust-z/percentile-rank heatmaps for both variables, derived structure (log10(dust), sol-to-sol temperature difference, vertical temperature gradient) and DDR1 min/mean/max perturbation traces, and retrieval-quality/coverage panels (valid-vs-failed profile counts, per-profile failure-level box plots, and a combined profile-count/invalid-retrieval panel) built from the raw s01 catalogs so retrieval failures dropped by s04 are still visible. Unlike the light collage, it processes many storm instances at once (filtered the same way as the timeseries scripts) and parallelizes over them with joblib, writing one PNG per instance per requested LTST window.
python mercap2/viz/plot_storm_collages_v2.py \
--output_dir /path/to/mercap2/jgr_plots/storm_collages --ltst_window both \
--sol_range -20,21 --conflev_min 2 --n_profiles_min 10
mercap2/viz/unadapted_from_v1/ keeps the remaining v1 scripts (per-distance perturbation plots, older filtersheet-era timeseries/violin plots, MDGM outline burning) as reference; they still read the v1 database/filtersheet and have not been ported.
Unit tests live in test/ (one module per pipeline stage plus the utils) and run against small synthetic fixtures with no database, network, or large-file dependencies:
pytest test/
For end-to-end checks against real data, use each stage's --smoke_test flag instead.