Changeset 3819 for trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py

Timestamp:

Jul 1, 2025, 11:39:58 AM (5 days ago)

Author:

jbclement

Message:

PEM:
Addition of 2 figures in the outputs of "visu_evol_layering.py": (i) the stratification over time colored using RGB ternary mix of H2O ice, CO2 ice and dust; (ii) the dust-to-ice ratio in the stratification over time.
JBC

File:

• trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py (modified) (6 diffs)

trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py

@@ -11 +11 @@
 from netCDF4 import Dataset
 import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1.inset_locator import inset_axes
+from matplotlib.colors import LinearSegmentedColormap, LogNorm
 from scipy.interpolate import interp1d
 
@@ -414 +416 @@
             fig.suptitle(
                 f"Content variation over time for (Grid point {ig+1}, Slope {isl+1})",
-                fontsize=14
+                fontsize=14,
+                fontweight='bold'
             )
 
@@ -456 +459 @@
 
             fig.subplots_adjust(right=0.88)
+            fig.tight_layout(rect=[0, 0, 0.88, 1.0])
             cbar_ax = fig.add_axes([0.90, 0.15, 0.02, 0.7])
             fig.colorbar(im, cax=cbar_ax, orientation='vertical', label="Content")
-            fig.tight_layout(rect=[0, 0, 0.88, 1.0])
 
             fname = os.path.join(
                 output_folder, f"layering_evolution_ig{ig+1}_is{isl+1}.png"
+            )
+            fig.savefig(fname, dpi=150)
+
+
+def plot_stratification_rgb_over_time(
+    gridded_data,
+    ref_grid,
+    top_index,
+    heights_data,
+    date_time,
+    exclude_sub=False,
+    output_folder="."
+):
+    """
+    Plot stratification over time colored using RGB ternary mix of H2O ice (blue), CO2 ice (violet), and dust (orange).
+    Includes a triangular legend showing the mix proportions.
+    """
+
+    # Define constant colors
+    violet = np.array([255, 0, 255], dtype=float) / 255
+    blue   = np.array([0, 0, 255], dtype=float) / 255
+    orange = np.array([255, 165, 0], dtype=float) / 255
+
+    # Prepare elevation mask
+    mask_elev = (ref_grid >= 0.0) if exclude_sub else np.ones_like(ref_grid, dtype=bool)
+    elev = ref_grid[mask_elev]
+
+    # Generate legend image once
+    res = 300
+    u = np.linspace(0, 1, res)
+    v = np.linspace(0, np.sqrt(3)/2, res)
+    X, Y = np.meshgrid(u, v)
+    V_bary = 2 * Y / np.sqrt(3)
+    U_bary = X - 0.5 * V_bary
+    W_bary = 1 - U_bary - V_bary
+    mask_triangle = (U_bary >= 0) & (V_bary >= 0) & (W_bary >= 0)
+
+    legend_rgb = (
+        U_bary[..., None] * violet
+        + V_bary[..., None] * orange
+        + W_bary[..., None] * blue
+    )
+    legend_rgb = np.clip(legend_rgb, 0.0, 1.0)
+    legend_rgba = np.zeros((res, res, 4))
+    legend_rgba[..., :3] = legend_rgb
+    legend_rgba[..., 3] = mask_triangle.astype(float)
+
+    # Loop over grid and slope
+    h2o = gridded_data['h2o_ice']
+    co2 = gridded_data['co2_ice']
+    dust = gridded_data['dust']
+    ngrid, ntime, nslope, nz = h2o.shape
+
+    for ig in range(ngrid):
+        for isl in range(nslope):
+            # Compute RGB stratification over time
+            rgb = np.ones((nz, ntime, 3), dtype=float)
+            for t in range(ntime):
+                mask_z = np.arange(nz) < top_index[ig, t, isl]
+                if not mask_z.any():
+                    continue
+                cH2O = np.clip(h2o[ig, t, isl, mask_z], 0, None)
+                cCO2 = np.clip(co2[ig, t, isl, mask_z], 0, None)
+                cDust = np.clip(dust[ig, t, isl, mask_z], 0, None)
+                total = cH2O + cCO2 + cDust
+                total[total == 0] = 1.0
+                fH2O = cH2O / total
+                fCO2 = cCO2 / total
+                fDust = cDust / total
+                mix = (
+                    np.outer(fH2O, blue)
+                    + np.outer(fCO2, violet)
+                    + np.outer(fDust, orange)
+                )
+                mix = np.clip(mix, 0.0, 1.0)
+                rgb[mask_z, t, :] = mix
+
+            display_rgb = rgb[mask_elev, :, :]
+
+            # Create figure with legend
+            fig, (ax_main, ax_leg) = plt.subplots(
+                1, 2, figsize=(12, 5), dpi=200,
+                gridspec_kw={'width_ratios': [5, 1]}
+            )
+
+            # Main stratification panel
+            ax_main.imshow(
+                display_rgb,
+                aspect='auto',
+                extent=[date_time[0], date_time[-1], elev.min(), elev.max()],
+                interpolation='nearest',
+                origin='lower'
+            )
+            ax_main.set_facecolor('white')
+            ax_main.set_title(f"Ternary mix over time (Grid point {ig+1}, Slope {isl+1})", fontweight='bold')
+            ax_main.set_xlabel("Time (Mars years)")
+            ax_main.set_ylabel("Elevation (m)")
+
+            # Legend panel
+            ax_leg.imshow(
+                legend_rgba,
+                extent=[0, 1, 0, np.sqrt(3)/2],
+                origin='lower',
+                interpolation='nearest'
+            )
+
+            # Draw triangle border
+            triangle = np.array([[0, 0], [1, 0], [0.5, np.sqrt(3)/2], [0, 0]])
+            ax_leg.plot(triangle[:, 0], triangle[:, 1], 'k-', linewidth=1)
+
+            # Dashed gridlines
+            ticks = np.linspace(0.25, 0.75, 3)
+            for f in ticks:
+                ax_leg.plot([1 - f, 0.5 * (1 - f)], [0, (1 - f)*np.sqrt(3)/2], '--', color='k', linewidth=0.5)
+                ax_leg.plot([f, f + 0.5 * (1 - f)], [0, (1 - f)*np.sqrt(3)/2], '--', color='k', linewidth=0.5)
+                y = (np.sqrt(3)/2) * f
+                ax_leg.plot([0.5 * f, 1 - 0.5 * f], [y, y], '--', color='k', linewidth=0.5)
+
+            # Legend labels
+            ax_leg.text(0, -0.05, 'H2O ice', ha='center', va='top', fontsize=8)
+            ax_leg.text(1, -0.05, 'CO2 ice', ha='center', va='top', fontsize=8)
+            ax_leg.text(0.5, np.sqrt(3)/2 + 0.05, 'Dust', ha='center', va='bottom', fontsize=8)
+            ax_leg.axis('off')
+
+            plt.tight_layout()
+
+            # Save figure
+            fname = os.path.join(output_folder, f"layering_rgb_evolution_ig{ig+1}_is{isl+1}.png")
+            fig.savefig(fname, dpi=150, bbox_inches='tight')
+
+
+def plot_dust_to_ice_ratio_over_time(
+    gridded_data,
+    ref_grid,
+    top_index,
+    heights_data,
+    date_time,
+    exclude_sub=False,
+    output_folder="."
+):
+    """
+    Plot the dust-to-ice ratio in the stratification over time,
+    using a blue-to-orange colormap:
+    - blue: ice-dominated (low dust-to-ice ratio)
+    - orange: dust-dominated (high dust-to-ice ratio)
+    """
+    h2o = gridded_data['h2o_ice']
+    dust = gridded_data['dust']
+    ngrid, ntime, nslope, nz = h2o.shape
+
+    # Elevation mask
+    if exclude_sub:
+        elevation_mask = (ref_grid >= 0.0)
+        elev = ref_grid[elevation_mask]
+    else:
+        elevation_mask = np.ones_like(ref_grid, dtype=bool)
+        elev = ref_grid
+
+    # Define custom blue-to-orange colormap
+    blue = np.array([0, 0, 255], dtype=float) / 255
+    orange = np.array([255, 165, 0], dtype=float) / 255
+    custom_cmap = LinearSegmentedColormap.from_list('BlueOrange', [blue, orange], N=256)
+
+    # Log‑ratio bounds and small epsilon to avoid log(0)
+    vmin, vmax = -2, 1
+    epsilon = 1e-6
+
+    # Loop over grids and slopes
+    for ig in range(ngrid):
+        for isl in range(nslope):
+            log_ratio_array = np.full((nz, ntime), np.nan, dtype=np.float32)
+
+            # Compute log10(dust/ice) profile at each time step
+            for t in range(ntime):
+                zmax = top_index[ig, t, isl]
+                if zmax <= 0:
+                    continue
+
+                h2o_profile = np.clip(h2o[ig, t, isl, :zmax], 0, None)
+                dust_profile = np.clip(dust[ig, t, isl, :zmax], 0, None)
+
+                with np.errstate(divide='ignore', invalid='ignore'):
+                    ratio_profile = np.where(
+                        h2o_profile > 0,
+                        dust_profile / h2o_profile,
+                        10**(vmax + 1)
+                    )
+                    log_ratio = np.log10(ratio_profile + epsilon)
+                    log_ratio = np.clip(log_ratio, vmin, vmax)
+
+                log_ratio_array[:zmax, t] = log_ratio
+
+            # Convert back to linear ratio and apply elevation mask
+            ratio_array = 10**log_ratio_array
+            ratio_display = ratio_array[elevation_mask, :]
+
+            # Plot
+            fig, ax = plt.subplots(figsize=(8, 6), dpi=150)
+            im = ax.imshow(
+                ratio_display,
+                aspect='auto',
+                extent=[date_time[0], date_time[-1], elev.min(), elev.max()],
+                origin='lower',
+                interpolation='nearest',
+                cmap='managua_r',
+                norm=LogNorm(vmin=10**vmin, vmax=10**vmax)
+            )
+
+            # Add colorbar with simplified ratio labels
+            cbar = fig.colorbar(im, ax=ax, orientation='vertical')
+            cbar.set_label('Dust / H₂O ice (ratio)')
+
+            # Define custom ticks and labels
+            ticks = [1e-2, 1e-1, 1, 1e1]
+            labels = ['1:100', '1:10', '1:1', '10:1']
+            cbar.set_ticks(ticks)
+            cbar.set_ticklabels(labels)
+
+            # Save figure
+            plt.tight_layout()
+            fname = os.path.join(
+                output_folder,
+                f"dust_to_ice_ratio_grid{ig+1}_slope{isl+1}.png"
             )
             fig.savefig(fname, dpi=150)
@@ -493 +719 @@
             fig.suptitle(
                 f"Strata count & total height over time for (Grid point {ig+1}, Slope {isl+1})",
-                fontsize=14
+                fontsize=14,
+                fontweight='bold'
             )
 
@@ -549 +776 @@
     # Plot
     fig, axes = plt.subplots(3, 1, figsize=(8, 10), sharex=True)
-    fig.suptitle("Orbital Parameters vs Simulated Time", fontsize=14)
+    fig.suptitle("Orbital Parameters vs Simulated Time", fontsize=14, fontweight='bold')
 
     axes[0].plot(date_time, obl_interp, 'r+', linestyle='-')
@@ -626 +853 @@
         print(f"Warning: {date_time.size} timestamps vs {ntime} NetCDF files.")
 
-    # 13) Plot stratification over time
+    # 13) Plot stratification data over time
     plot_stratification_over_time(
         gridded_data, ref_grid, top_index, heights_data, date_time,
         exclude_sub=exclude_sub, output_folder="."
     )
-
-    # 14) Plot strata count & total height vs time
+    plot_stratification_rgb_over_time(
+        gridded_data, ref_grid, top_index, heights_data, date_time,
+        exclude_sub=exclude_sub, output_folder="."
+    )
+    plot_dust_to_ice_ratio_over_time(
+        gridded_data, ref_grid, top_index, heights_data, date_time,
+        exclude_sub=exclude_sub, output_folder="."
+    )
     plot_strata_count_and_total_height(heights_data, date_time, output_folder=".")
 
-    # 15) Plot orbital parameters
+    # 14) Plot orbital parameters
     plot_orbital_parameters(infofile, orbfile, date_time, output_folder=".")
 
-    # 16) Show all figures
+    # 15) Show all figures
     plt.show()