allow numbering for small fields

Author: bobmyhill

contrib/perplex/create_phase_diagram.py

@@ -45,67 +45,79 @@ def run_perplex(
     run_pssect(perplex_bindir, project_name, convert_to_pdf=False, verbose=False)
 
 
+# The following compositions are from Xu et al. (2008)
+# (https://doi.org/10.1016/j.epsl.2008.08.012)
+compositions = {
+    "pyrolite": burnman.Composition(
+        {
+            "SiO2": 38.71,
+            "MgO": 49.85,
+            "FeO": 6.17,
+            "CaO": 2.94,
+            "Al2O3": 2.22,
+            "Na2O": 0.11,
+        },
+        "molar",
+    ),
+    "basalt": burnman.Composition(
+        {
+            "SiO2": 51.75,
+            "MgO": 14.94,
+            "FeO": 7.06,
+            "CaO": 13.88,
+            "Al2O3": 10.19,
+            "Na2O": 2.18,
+        },
+        "molar",
+    ),
+    "harzburgite": burnman.Composition(
+        {
+            "SiO2": 36.07,
+            "MgO": 56.51,
+            "FeO": 6.07,
+            "CaO": 0.81,
+            "Al2O3": 0.53,
+            "Na2O": 0.001,
+        },
+        "molar",
+    ),
+    "modified_harzburgite": burnman.Composition(
+        {
+            "SiO2": 36.04,
+            "MgO": 56.54,
+            "FeO": 5.97,
+            "CaO": 0.79,
+            "Al2O3": 0.65,
+            "Na2O": 0.001,
+        },
+        "molar",
+    ),
+}
+
+
 if __name__ == "__main__":
     # Define the project parameters
-    compositions = {
-        "pyrolite": burnman.Composition(
-            {
-                "SiO2": 38.71,
-                "MgO": 49.85,
-                "FeO": 6.17,
-                "CaO": 2.94,
-                "Al2O3": 2.22,
-                "Na2O": 0.11,
-            },
-            "molar",
-        ),
-        "basalt": burnman.Composition(
-            {
-                "SiO2": 51.75,
-                "MgO": 14.94,
-                "FeO": 7.06,
-                "CaO": 13.88,
-                "Al2O3": 10.19,
-                "Na2O": 2.18,
-            },
-            "molar",
-        ),
-        "harzburgite": burnman.Composition(
-            {
-                "SiO2": 36.07,
-                "MgO": 56.51,
-                "FeO": 6.07,
-                "CaO": 0.81,
-                "Al2O3": 0.53,
-                "Na2O": 0.001,
-            },
-            "molar",
-        ),
-        "modified_harzburgite": burnman.Composition(
-            {
-                "SiO2": 36.04,
-                "MgO": 56.54,
-                "FeO": 5.97,
-                "CaO": 0.79,
-                "Al2O3": 0.65,
-                "Na2O": 0.001,
-            },
-            "molar",
-        ),
-    }
     project_name = "basalt"
     database = databases["stx21"]
     composition = compositions[project_name]
     pressure_range_total = [1.0e5, 140.0e9]
     temperature_range_total = [200.0, 3000.0]
-    n_pressures_per_split = 241
-    n_temperatures_per_split = 201
-    n_splits_pressure = 7
-    n_splits_temperature = 3
 
-    create_plots = True
-    outfile_base = f"{project_name}_table"
+    # Split pressure and temperature so that PerpleX
+    # no temperature splits seems to make diagram with
+    # less prominent discontinuities
+    n_pressures_per_split = 121
+    n_temperatures_per_split = 601
+    n_splits_pressure = 14
+    n_splits_temperature = 1
+
+    # If this script has already been run, and you just want to
+    # tweak the figure, perplex should not be run again.
+    perplex_should_be_run = True
 
+    # End of project definitions
+
+    outfile_base = f"{project_name}_table"
     perplex_dir = os.path.join(os.getcwd(), "perplex-installer/Perple_X")
     perplex_bindir = os.path.join(os.getcwd(), "perplex-installer/Perple_X/bin")
 
@@ -114,7 +126,7 @@ def run_perplex(
         temperature_range_total[1] - temperature_range_total[0]
     ) / n_splits_temperature
 
-    if create_plots:
+    if perplex_should_be_run:
         # Check if the project directory already exists
         # If it does, delete it to start fresh.
         try:
@@ -249,19 +261,26 @@ def run_perplex(
     fig = plt.figure(figsize=(12, 8))
     ax = fig.add_subplot(1, 1, 1)
 
-    pretty_plot_phase_diagram(
+    small_fields = pretty_plot_phase_diagram(
         ax,
         modes_filenames,
         phase_name_replacements,
         bounding_colors,
+        n_phases_bounds=[2.0, 7.0],
         smoothing_window=8,
         smoothing_order=1,
         linewidth=0.5,
         label_scaling=3.0,
         label_clearance=0.01,
-        n_phases_bounds=[2.0, 7.0],
+        number_small_fields=True,
     )
 
+    with open(f"{project_name}_field_ids.txt", "w") as f:
+        for i, small_field in enumerate(small_fields):
+            line = f"{i+1}: {small_field}"
+            print(line)
+            f.write(f"{line}\n")
+
     fig.savefig(f"{project_name}.pdf")
     plt.show()
 

contrib/perplex/perplex_utils.py

@@ -1060,12 +1060,13 @@ def pretty_plot_phase_diagram(
     werami_mode_tab_filenames,
     phase_name_replacements,
     bounding_colors=["#44015a", "#ffffff"],
+    n_phases_bounds=[2, 10],
     smoothing_window=4,
     smoothing_order=1,
     linewidth=0.5,
     label_scaling=3.0,
     label_clearance=0.01,
-    n_phases_bounds=[2, 10],
+    number_small_fields=True,
 ):
     """
     Plot the Perple_X calculated phase diagram on a matplotlib axis.
@@ -1078,17 +1079,26 @@ def pretty_plot_phase_diagram(
     :type phase_name_replacements: dict[str, str]
     :param bounding_colors: List of two colors for bounding colormap.
     :type bounding_colors: [str, str]
+    :param n_phases_bounds: Minimum and maximum number of phases for the colormap.
+    :type n_phases_bounds: [int, int]
     :param smoothing_window: Savitzky-Golay smoothing window length.
     :type smoothing_window: int
     :param smoothing_order: Savitzky-Golay smoothing polynomial order.
     :type smoothing_order: int
+    :param linewidth: Linewidth for the edges of fields.
+    :type linewidth: float
     :param label_scaling: Scaling factor corresponding to an
     approximate value of the width:height ratio of the label.
     :type label_scaling: float
     :param label_clearance: Only plot the label if there is a certain
-    clearance to the edge of the field. Scaled to the height of the domain.
+    clearance between the center and the edge of the field.
+    Scaled to the height of the domain.
     :type label_clearance: float
-    :return: None
+    :param number_small_fields: Replace assemblage labels with numbers if
+    there is not enough clearance around the label.
+    :type number_small_fields: bool
+    :return: List of assemblages corresponding to numbered small fields.
+    :rtype: [str]
     """
 
     polygon_data = []
@@ -1131,6 +1141,8 @@ def pretty_plot_phase_diagram(
     # ScalarMappable for colorbar
     sm = ScalarMappable(cmap=cmap, norm=norm)
 
+    small_fields = []
+    large_fields = set()
     for polygon in polygon_data:
 
         color = cmap(norm(polygon["n_phases"]))
@@ -1151,8 +1163,27 @@ def pretty_plot_phase_diagram(
         label_pos, label_dist = get_label_position_from_polygon_contour(
             contour, P_range, T_range, label_scaling
         )
+
         if label_dist > label_clearance:
+            large_fields.add(polygon["label"])
             ax.text(*label_pos, polygon["label"], fontsize=6, ha="center", va="center")
+        elif number_small_fields:
+            if len(polygon["edges"]) > 0:
+                small_fields.append([label_pos, polygon["label"]])
+
+    unique_ids = {}
+    unique_small_fields = []
+    i = 1
+    for _, label in small_fields:
+        if label not in unique_ids.keys() and label not in large_fields:
+            unique_ids[label] = i
+            unique_small_fields.append(label)
+            i = i + 1
+
+    for label_pos, label in small_fields:
+        ax.text(
+            *label_pos, f"{unique_ids[label]:02d}", fontsize=6, ha="center", va="center"
+        )
 
     cbar = plt.colorbar(sm, ax=ax, boundaries=boundaries, ticks=n_phases_range)
     cbar.set_label("n phases")
@@ -1164,3 +1195,5 @@ def pretty_plot_phase_diagram(
     ax.set_ylim(bounds[1])
     ax.grid(True, linestyle="--", color="gray", alpha=0.7)
     ax.set_axisbelow(False)
+
+    return unique_small_fields