Add plot scripts for map and CDF for single field

Author: mkstratos

detclim/plot_ensemble_cdf.py

@@ -0,0 +1,95 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import seaborn as sns
+import xarray as xr
+
+import detclim.results_plot as rplt
+
+
+def main(dvar):
+    scratch_dir = Path("/lustre/storm/nwp501/proj-shared/mkelleher/e3sm_scratch")
+    ctl = "ctl"
+    test = "clubb_c1-5p0pct"
+
+    cases = {
+        ctl: "2025-09-16.F2010.ne30pg2_r05_oECv3_control",
+        test: "2025-09-16.F2010.ne30pg2_r05_oECv3_clubb_c1_2p520000",
+    }
+
+    ens_dirs = {_ens: Path(scratch_dir, _name, "run") for _ens, _name in cases.items()}
+
+    file_pattern = "{case}.eam.h0.aavg.nc"
+    ens_file = {
+        _ens: Path(ens_dirs[_ens], file_pattern.format(case=_name))
+        for _ens, _name in cases.items()
+    }
+    data = {_ens: xr.open_dataset(_file) for _ens, _file in ens_file.items()}
+
+    _figwidth = 10 / 2.54
+    _figheight = _figwidth * 0.9
+
+    data_ref = data[ctl][dvar].mean(dim="time")
+    data_test = data[test][dvar].mean(dim="time")
+
+    n_q = data_ref.shape[0] // 2
+
+    all_ = np.concatenate((data_test.values, data_ref.values))
+    min_all = np.min(all_)
+    max_all = np.max(all_)
+
+    norm_ref = (data_ref - min_all) / (max_all - min_all)
+    norm_test = (data_test - min_all) / (max_all - min_all)
+
+    bins = np.linspace(0, 1, n_q)
+
+    freq_ref, _ = np.histogram(norm_ref, bins=bins)
+    freq_test, _ = np.histogram(norm_test, bins=bins)
+
+    cdf_ref = freq_ref.cumsum()
+    cdf_test = freq_test.cumsum()
+
+    _xlabel = "Normalized annual global means"
+    _ylabel = "N Ensemble members"
+
+    df_all = pd.concat(
+        [
+            pd.DataFrame(
+                {
+                    _xlabel: bins,
+                    _ylabel: [0, *cdf_ref],
+                    "Simulation": rplt.fmt_case(ctl),
+                }
+            ),
+            pd.DataFrame(
+                {
+                    _xlabel: bins,
+                    _ylabel: [0, *cdf_test],
+                    "Simulation": rplt.fmt_case(test),
+                }
+            ),
+        ]
+    )
+
+    fig, axes = plt.subplots(
+        1,
+        1,
+        figsize=(_figwidth, _figheight),
+        dpi=300,
+    )
+    sns.set_theme(
+        context="paper",
+        style="whitegrid",
+        font="sans-serif",
+        rc={"lines.linewidth": 1.5, "font.size": 8},
+    )
+    sns.lineplot(x=_xlabel, y=_ylabel, hue="Simulation", data=df_all, ax=axes)
+    plt.tight_layout()
+    fig.savefig(f"plt_cdf_{dvar}.pdf")
+    fig.savefig(f"plt_cdf_{dvar}.png")
+
+
+if __name__ == "__main__":
+    main(dvar="CLDLIQ")

detclim/plot_ensemble_mean_map.py

@@ -0,0 +1,130 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+import xarray as xr
+from cartopy import crs as ccrs
+
+import detclim.results_plot as rplt
+
+
+def main():
+    scratch_dir = Path("/lustre/storm/nwp501/proj-shared/mkelleher/e3sm_scratch")
+    ctl = "ctl"
+    test = "clubb_c1-5p0pct"
+    params_hum = {
+        "effgw_oro": "GW Orog",
+        "clubb_c1": "CLUBB C1",
+        "zmconv_c0_ocn": "ZM Conv C0-Ocean",
+        "ctl": "Control",
+    }
+
+    cases = {
+        ctl: "2025-09-16.F2010.ne30pg2_r05_oECv3_control",
+        test: "2025-09-16.F2010.ne30pg2_r05_oECv3_clubb_c1_2p520000",
+    }
+
+    ens_dirs = {_ens: Path(scratch_dir, _name, "run") for _ens, _name in cases.items()}
+
+    file_pattern = "{case}.eam.h0.aavg.nc"
+    ens_file = {
+        _ens: Path(ens_dirs[_ens], file_pattern.format(case=_name))
+        for _ens, _name in cases.items()
+    }
+
+    files_2d = {
+        _ens: Path(
+            ens_file[_ens].parent,
+            f"{cases[_ens]}.eam.h0.0001-03_0002-02_climo_ANN_r05_ensavg.nc",
+        )
+        for _ens in cases
+    }
+
+    data_2d = {_ens: xr.open_dataset(_file) for _ens, _file in files_2d.items()}
+
+    proj = ccrs.Robinson()
+    src_crs = ccrs.PlateCarree()
+
+    dvar = "CLDLIQ"
+    title_props = {"fontsize": 8}
+    _figwidth = 12.5 / 2.54
+    plot_3panel = False
+
+    cbar_props = {"location": "right", "pad": 0.01, "shrink": 0.8}
+    _figheight = _figwidth * 1.1
+    fig = plt.figure(
+        figsize=(_figwidth, _figheight),
+        dpi=300,
+    )
+    axes = [
+        fig.add_subplot(
+            2,
+            1,
+            idx + 1,
+            projection=proj,
+        )
+        for idx in range(2)
+    ]
+
+    diff_data = {}
+    for _ens, _data in data_2d.items():
+        _pltdata = _data.isel(time=0)
+        if "lev" in _pltdata.coords:
+            _pltdata = _pltdata.isel(lev=-1)
+        else:
+            _pltdata = _pltdata.squeeze()
+        _pltdata = _pltdata.eval(dvar)
+        if dvar in _data:
+            _name = _data[dvar].long_name
+            _units = _data[dvar].units
+        diff_data[_ens] = _pltdata
+
+    _cf = axes[0].pcolormesh(
+        _data["lon"],
+        _data["lat"],
+        diff_data[ctl],
+        transform=src_crs,
+        cmap="plasma",
+        rasterized=True,
+    )
+    axes[0].set_title(f"{rplt.fmt_case(ctl)} mean", fontsize=title_props["fontsize"])
+    _cax = fig.colorbar(_cf, ax=axes[0], **cbar_props)
+
+    _cax.set_label(_units, fontsize=title_props["fontsize"])
+
+    diff_2d = diff_data[test] - diff_data[ctl]
+    _absmax = np.max(np.abs(diff_2d.quantile([0.01, 0.99])))
+    _cf = axes[1].pcolormesh(
+        data_2d[test]["lon"],
+        data_2d[test]["lat"],
+        diff_2d,
+        vmin=-_absmax,
+        vmax=_absmax,
+        cmap="BrBG",
+        transform=src_crs,
+        rasterized=True,
+    )
+    axes[1].set_title(
+        f"Difference {rplt.fmt_case(test)} - {rplt.fmt_case(ctl)}",
+        fontsize=title_props["fontsize"],
+    )
+    _cax = fig.colorbar(_cf, ax=axes[1], **cbar_props)
+    _cax.set_label(_units, fontsize=title_props["fontsize"])
+    for axis in axes:
+        axis.coastlines(linewidth=0.5)
+        axis.gridlines(linestyle="--", linewidth=0.5, color="k")
+
+    fig.suptitle(f"{_name} [{_units}]", fontsize=12)
+    plt.tight_layout()
+    plt.savefig(
+        f"test_diff_{_name.replace(' ', '_').replace('(', '').replace(')', '')}_"
+        f"{ctl}-{test}.png"
+    )
+    plt.savefig(
+        f"test_diff_{_name.replace(' ', '_').replace('(', '').replace(')', '')}_"
+        f"{ctl}-{test}.pdf"
+    )
+
+
+if __name__ == "__main__":
+    main()