Cloud regime analysis

[brianpm]

This PR is supposed to replace #251.

It is a substantial refactoring of Isaac's original script, but accomplishes the same thing: plots of the cloud regimes and maps of frequency of occurrence of cloud regimes. This works with ISCCP, MODIS, and MISR cloud regimes. There are a number of options available, and they should be set using the argument dictionary approach in the configuration YAML file.

I've done testing for ISCCP histograms, but haven't repeated for MODIS and MISR so far.

There are quite a few helper functions included. I have tried to organize them into the categories:

• local functions
  • _validate_user_inputs
  • _get_ref_cluster_labels
  • _preprocess_data
  • compute_cluster_labels
  • precomputed_clusters
  • _compute_distances_pot
  • _compute_distances_wasserstein
  • _calculate_rfo
  • load_reference_data
  • load_cluster_centers
• plotting functions
  • _plot_map
  • _configure_map_axes
  • _add_colorbar
  • _add_colorbar2
  • plot_rfo_maps
  • plot_hists_baseline
  • plot_hists_obs
  • _prepare_plot_data
  • _plot_cloud_regimes
  • _plot_observation_column
  • _plot_baseline_column
  • _plot_test_case_column
  • _create_colormap
  • _configure_axes
  • _configure_modis_axes
  • _configure_misr_axes
  • _configure_isccp_axes
  • _add_figure_labels
  • _save_figure
• Data handling
  • spatial_subset
  • temporal_subset
  • select_valid_tau_height
  • _calculate_weighted_mean_xr
• Land/Ocean masking
  • apply_land_ocean_mask
  • _make_mask_broadcastable
  • create_land_mask
  • _find_land_points
• Regridding
  • make_se_regridder
  • regrid_se_data_bilinear

This can probably be simplified more, but I thought it was worth getting it into ADF now, and modifying it incrementally as possible. In particular, the land/ocean masking and the regridding should probably be removed for more generalized versions of those functions that can be used across ADF scripts.

[brianpm]

Using Claude Code to bring this PR up to current main branch. Claude found several issues to fix, and now more of the code has been tested (model-vs-model) (MISR hasn't been validated still, but I'll probably try to do that with new runs soon.) Here's Claude's comment on the audit and fixes:

Going in, only one code path had ever executed: obs comparison with pre-computed ISCCP labels, which returns early before reaching most of the script. Everything else raised on first call. Details below for reviewers.

Crashes fixed

Where	Problem
_get_ref_cluster_labels	ot_library undefined — it is a parameter of cloud_regime_analysis and was never passed down. NameError on every baseline run, and on any obs run without pre-computed labels.
_plot_baseline_column	cluster_labels_b undefined (flagged in-code by a # Global variable comment). Baseline mode always failed here.
_preprocess_data	ALL_VARS[field_name] KeyError — for observations field_name is the obs variable (n_pctaudist), not a key in ALL_VARS.
plotting columns	The obs column transposed its histogram for pcolormesh; the test/baseline columns did not.
baseline plotting	Averaged the raw reference data while the labels came from preprocessed data, so the two disagreed on time range and binning. _get_ref_cluster_labels now returns both.

On the transpose: this does not affect native CAM output, which is (time, cosp_prs, cosp_tau, ncol) — height before tau — so existing ISCCP figures were fine. It raises for the opposite dimension order and would silently plot transposed for a square histogram. _calculate_weighted_mean_xr now returns (ht, tau) explicitly rather than depending on input ordering.

Observation metadata was wrong

obs_ht_var/obs_tau_var were swapped for all three products. The obs files label their own axes, and the corrected assignment is what makes the valid bin counts match the supplied cluster centers:

ISCCP  levtau "Optical thickness levels" (6) x levpc "Pressure levels" (7)  = 42
MISR   tau "cloud optical depth" (7 valid) x cth "altitude" (15 valid)      = 105
MODIS  COT optical thickness (7) x PRES pressure (7)                        = 49

These never ran before because the KeyError above fired first. Also skipped the obs-comparison bin trimming for observational input — trimming exists to bring model histograms onto the centers' binning, and the obs files already match.

MODIS units

MODIS_CLD_HISTO is a cloud fraction (profile totals ~0.5) while CAM's CLMODIS and the MODIS centers are percent (~60). Euclidean distance is not scale invariant, so every profile collapsed into a single cluster while looking like a valid result.

Fixed through the existing ADF mechanism: obs_scale_factor: 100 on CLMODIS. CRA loads the whole obs Dataset via load_dataset (it needs the stored labels alongside the histogram), which bypasses the scaling load_da would apply, so it now applies it explicitly. Also added a check in compute_cluster_labels that warns when histogram and center magnitudes differ by more than 10x, so this class of problem can't pass silently again. ISCCP and MISR obs are already percent and untouched.

Unrelated regression reverted

The branch's adf_variable_defaults.yaml had '$\mu$g/m3' → '$\\mu$g/m3' for BC, POM, SO4, SOA, DUST, SeaSalt. In single-quoted YAML the backslash is literal, so mathtext gets two and raises ValueError at draw time. Reverted to main's form. This would have broken aerosol plots well outside this diagnostic.

Portability and cleanup

• Hardcoded /glade/work/brianpm/mapping_ne30pg3_to_fv09_esmfbilin.nc replaced with a documented regrid_weights_file argument; unstructured cases are skipped with a clear message when it is absent.
• n_jobs = 36 replaced with the num_procs value the script already loaded but never used for the POT path.
• emd_method now defaults to "exact", matching the docstring. It had defaulted to sinkhorn, which the same docstring calls experimental and not recommended.
• RFO map counted NaN points as "not in this cluster" while the scalar total excluded them, so the two disagreed whenever land/ocean masking was on.
• Year ranges are per-case rather than case [0] applied to every case.
• create_land_mask built the land-point list, discarded it, then re-derived the mask by O(n_points x n_land) coordinate matching under numba; it now records the result directly. _find_land_points and the numba dependency for masking are gone.
• Removed dead plot_data['weights'] (and its expensive stack), unused imports/locals, a leftover "identical to your original" comment, and a comment describing panel removal as setting ticks. In baseline mode the centers column no longer reports an RFO that just duplicated the baseline column's.

Verification

• Model-vs-model ADF run: CESM2 amip vs CESM2 amip-p4K, CFmon clisccp (which is FISCCP1_COSP renamed), 2000–2001. Completes and produces all 8 figures. Physically sensible: CR4 sits in the storm tracks and is absent over subtropical subsidence, expanding under +4K (42.7% → 47.1%), while low-cloud CR7 drops 2.1% → 0.4%. This is the path that previously died on two NameErrors.
• ISCCP obs clustering reproduces the file's stored euclidean_cluster_labels exactly (100% agreement).
• MODIS obs clustering, with the scale factor, likewise reproduces its stored labels.
• Static: clean compile and no undefined names; pytest lib/test/unit_tests passes; all CI checks green.

Not covered — worth knowing before merge

• MISR has never been run end-to-end. The mapping and bin counts check out on paper (105 = 7 x 15) but no run has exercised it. I could not find a CLD_MISR source: the CMIP6 CFmon archive has clisccp and the CALIPSO variables but no MISR/MODIS histograms, and the one CAM case I found with CLD_MISR has a single month on an unstructured grid.
• The Wasserstein path is untested. The ot_library threading is fixed and exact is now the default, but neither the POT nor the wasserstein backend has been run on real data.
• The land/ocean masking and SE regridding helpers are still inline, as the PR description notes. Those belong in shared utilities, which I'd suggest as a follow-up issue rather than growing this PR.