fix(plotnine): address review feedback for stereonet-equal-area

Attempt 2/3 - fixes based on AI review

Author: github-actions[bot]

plots/stereonet-equal-area/implementations/plotnine.py

@@ -1,31 +1,33 @@
-""" pyplots.ai
+"""pyplots.ai
 stereonet-equal-area: Structural Geology Stereonet (Equal-Area Projection)
 Library: plotnine 0.15.3 | Python 3.14.3
 Quality: 79/100 | Created: 2026-03-15
 """
 
 import matplotlib
-import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
 from plotnine import (
     aes,
+    after_stat,
     coord_fixed,
     element_blank,
     element_text,
+    geom_density_2d,
     geom_path,
     geom_point,
     geom_segment,
     geom_text,
     ggplot,
+    guides,
     labs,
+    scale_alpha_continuous,
     scale_color_manual,
     scale_shape_manual,
     scale_x_continuous,
     scale_y_continuous,
     theme,
 )
-from scipy.ndimage import gaussian_filter
 
 
 matplotlib.use("Agg")
@@ -57,7 +59,7 @@
 pole_y = pole_r * np.cos(pole_trend)
 poles_df = pd.DataFrame({"x": pole_x, "y": pole_y, "feature_type": feature_types})
 
-# Compute great circles (sample representative planes per type)
+# Compute great circles for representative planes
 gc_rows = []
 gc_indices = {"Bedding": [0, 10, 20, 30], "Fault": [40, 48, 56], "Joint": [65, 75, 85]}
 gc_id = 0
@@ -72,17 +74,15 @@
         dip_vec = np.array(
             [np.sin(dip_dir_rad) * np.cos(dip_rad), np.cos(dip_dir_rad) * np.cos(dip_rad), -np.sin(dip_rad)]
         )
-        angles = np.linspace(-np.pi / 2, np.pi / 2, 181)
-        for a in angles:
+        for a in np.linspace(-np.pi / 2, np.pi / 2, 181):
             pt = np.cos(a) * dip_vec + np.sin(a) * strike_vec
             if pt[2] > 0:
                 pt = -pt
             horiz = np.sqrt(pt[0] ** 2 + pt[1] ** 2)
             plunge = np.arctan2(-pt[2], horiz)
             trend = np.arctan2(pt[0], pt[1])
             r = np.sqrt(2) * np.sin((np.pi / 2 - plunge) / 2)
-            x = r * np.sin(trend)
-            y = r * np.cos(trend)
+            x, y = r * np.sin(trend), r * np.cos(trend)
             if x**2 + y**2 <= r_prim**2 * 1.01:
                 gc_rows.append({"x": x, "y": y, "feature_type": ftype, "gc_id": f"{ftype}_{gc_id}"})
         gc_id += 1
@@ -93,6 +93,24 @@
 circle_angles = np.linspace(0, 2 * np.pi, 361)
 prim_df = pd.DataFrame({"x": r_prim * np.cos(circle_angles), "y": r_prim * np.sin(circle_angles), "group": "primitive"})
 
+# Equal-area net grid lines (small circles at 30° dip intervals)
+grid_rows = []
+for dip_interval in range(30, 90, 30):
+    plunge_rad = np.radians(90 - dip_interval)
+    r_circle = np.sqrt(2) * np.sin((np.pi / 2 - plunge_rad) / 2)
+    for angle in np.linspace(0, 2 * np.pi, 181):
+        gx = r_circle * np.cos(angle)
+        gy = r_circle * np.sin(angle)
+        grid_rows.append({"x": gx, "y": gy, "grid_id": f"dip_{dip_interval}"})
+
+# Radial lines at 30° azimuth intervals
+for az in range(0, 360, 30):
+    az_rad = np.radians(az)
+    for t in np.linspace(0, r_prim, 50):
+        grid_rows.append({"x": t * np.sin(az_rad), "y": t * np.cos(az_rad), "grid_id": f"az_{az}"})
+
+grid_df = pd.DataFrame(grid_rows)
+
 # Degree tick marks every 10 degrees
 tick_rows = []
 for deg in range(0, 360, 10):
@@ -118,68 +136,63 @@
     if deg in (0, 90, 180, 270):
         continue
     rad = np.radians(deg)
-    offset = r_prim * 1.07
+    offset = r_prim * 1.08
     deg_labels.append({"x": offset * np.sin(rad), "y": offset * np.cos(rad), "label": f"{deg}°"})
 deg_label_df = pd.DataFrame(deg_labels)
 
-# Density contours (Kamb-style kernel density on projected pole coordinates)
-grid_res = 200
-gx_lin = np.linspace(-r_prim, r_prim, grid_res)
-gy_lin = np.linspace(-r_prim, r_prim, grid_res)
-gx_grid, gy_grid = np.meshgrid(gx_lin, gy_lin)
-density = np.zeros_like(gx_grid)
-
-sigma = 0.15
-for i in range(len(pole_x)):
-    dist_sq = (gx_grid - pole_x[i]) ** 2 + (gy_grid - pole_y[i]) ** 2
-    density += np.exp(-dist_sq / (2 * sigma**2))
-
-density = gaussian_filter(density, sigma=3)
-circle_mask = gx_grid**2 + gy_grid**2 > r_prim**2
-density[circle_mask] = np.nan
-
-# Extract contour coordinates (using matplotlib for numerical contour extraction only)
-fig_tmp, ax_tmp = plt.subplots()
-valid_density = density[~np.isnan(density)]
-levels = np.linspace(valid_density.max() * 0.2, valid_density.max() * 0.8, 5)
-cs = ax_tmp.contour(gx_lin, gy_lin, density, levels=levels)
-plt.close(fig_tmp)
-
-contour_rows = []
-contour_id = 0
-for level_segs in cs.allsegs:
-    for seg in level_segs:
-        if len(seg) > 0:
-            for xi, yi in seg:
-                if xi**2 + yi**2 <= r_prim**2 * 1.01:
-                    contour_rows.append({"x": xi, "y": yi, "contour_id": contour_id})
-            contour_id += 1
-
-contour_df = pd.DataFrame(contour_rows) if contour_rows else pd.DataFrame({"x": [], "y": [], "contour_id": []})
-
-# Colors (colorblind-safe: blue, orange, purple)
+# Cluster centroid annotations for data storytelling
+annotations = []
+for ftype in ["Bedding", "Fault", "Joint"]:
+    mask = poles_df["feature_type"] == ftype
+    cx, cy = poles_df.loc[mask, "x"].mean(), poles_df.loc[mask, "y"].mean()
+    mean_strike = strikes[np.array(feature_types) == ftype].mean()
+    mean_dip = dips[np.array(feature_types) == ftype].mean()
+    annotations.append({"x": cx, "y": cy - 0.12, "feature_type": ftype, "label": f"{mean_strike:.0f}°/{mean_dip:.0f}°"})
+annot_df = pd.DataFrame(annotations)
+
+# Colors (colorblind-safe)
 colors = {"Bedding": "#306998", "Fault": "#E5A023", "Joint": "#7B68A0"}
 shapes = {"Bedding": "o", "Fault": "D", "Joint": "s"}
 
-# Plot
+# Plot - using plotnine's geom_density_2d (stat_density_2d) for Kamb-style contouring
 plot = (
     ggplot()
+    # Grid lines (subtle equal-area net)
+    + geom_path(aes(x="x", y="y", group="grid_id"), data=grid_df, color="#CCCCCC", size=0.3, alpha=0.5)
+    # Primitive circle
     + geom_path(aes(x="x", y="y"), data=prim_df, color="#333333", size=1.2)
+    # Tick marks
     + geom_segment(aes(x="x1", y="y1", xend="x2", yend="y2"), data=tick_df, color="#333333", size=0.5)
-    + geom_path(
-        aes(x="x", y="y", group="contour_id"), data=contour_df, color="#777777", size=0.7, alpha=0.7, linetype="dashed"
+    # Density contours using plotnine's native stat_density_2d via geom_density_2d
+    # after_stat maps computed density level to alpha for visual depth
+    + geom_density_2d(
+        aes(x="x", y="y", alpha=after_stat("level")), data=poles_df, color="#666666", size=0.6, linetype="dashed"
     )
+    + scale_alpha_continuous(range=(0.3, 0.8))
+    + guides(alpha=False)
+    # Great circles
     + geom_path(aes(x="x", y="y", color="feature_type", group="gc_id"), data=gc_df, size=0.9, alpha=0.7)
+    # Poles to planes
     + geom_point(
         aes(x="x", y="y", color="feature_type", shape="feature_type"), data=poles_df, size=3.5, alpha=0.85, stroke=0.5
     )
+    # Cardinal directions
     + geom_text(aes(x="x", y="y", label="label"), data=dir_df, size=18, fontweight="bold", color="#222222")
-    + geom_text(aes(x="x", y="y", label="label"), data=deg_label_df, size=10, color="#555555")
+    # Degree labels (increased size for readability)
+    + geom_text(aes(x="x", y="y", label="label"), data=deg_label_df, size=13, color="#444444")
+    # Cluster orientation annotations
+    + geom_text(
+        aes(x="x", y="y", label="label", color="feature_type"),
+        data=annot_df,
+        size=11,
+        fontstyle="italic",
+        show_legend=False,
+    )
     + scale_color_manual(name="Feature Type", values=colors)
     + scale_shape_manual(name="Feature Type", values=shapes)
     + coord_fixed(ratio=1)
-    + scale_x_continuous(limits=(-1.8, 1.8))
-    + scale_y_continuous(limits=(-1.8, 1.8))
+    + scale_x_continuous(limits=(-1.85, 1.85))
+    + scale_y_continuous(limits=(-1.85, 1.85))
     + labs(title="stereonet-equal-area · plotnine · pyplots.ai")
     + theme(
         figure_size=(12, 12),